在EViews中实现数值求解常微分方程（ODE）

liwenhui

# S% j' g, d" h3 O) P* s
2 M- a+ d  f1 T- KEViews除了能解决计量经济学的估计问题以外，还提供一个编程环境用以解决复杂的问题。经过调试，我在EViews中实现了用龙格库塔方法求解常微分方程的数值解，供大家交流。
演示中，我使用了如下常微分方程作为测试：

这个方程的解析通解是：

4 l: v# ?1 |. q
使用“龙格库塔方法”，编制的EViews程序如下：

1. '用龙格库塔法求解常微分方程dydx=-y*cos(x)+exp(-sin(x)),y(0)=0在区间[0,10*pi]的数值解
   2 o; ]8 _! B. x5 [: v
2. '已知这个微分方程的解析解 y=exp(-sin(x))*x$ o, r\" h0 A/ v( d
   
3. 
   8 I+ R2 X. a8 S6 l- `4 K) V3 s; w  H
4. '生成一个workfile作为基本的数据容器
   1 O! z& z, \( b/ X8 |1 g
5. wfcreate (wf=temp) u 1000+ a' Z- }( {' S6 }7 u
   
6. 
   9 W: y! V$ L7 w$ K6 F# H\" h# f4 c, Z* `
7. '定义常量% [0 R+ J% ]6 E
   
8. scalar pi=3.141593 [) i' q3 c\" p6 D\" s; k* P
   
9. scalar a=0        '定义自变量下限8 W( `5 f# Q% C* `* o
   
10. scalar b=10*pi     '定义自变量上限
    & j, v( v6 u# W$ D. E2 w3 T
11. scalar  M=500       '定义步数
    8 l7 D8 f\" B) O* j: v( t
12. scalar h=(b-a)/M   '计算每步之间的间隔2 [& l2 O' ^\" q! |3 V
    
13. , g% p3 @: X% G% t$ }3 L
    
14. '定义一个矩阵来储存计算数据，其中第一列储存自变量数据，第二列储存因变量数据,第三列储存解析解的值用以作为比较
    ( d) W3 [3 S. W\" |/ n! {% f, u
15. matrix(M+1,3) F
    6 Y7 y; X+ y  A9 |
16. ; t  \/ Y: U' y, W9 w) l4 b; b4 m
    
17. '矩阵的第一行储存初值问题的初始条件# v. [' w: I5 W0 y2 n& [0 Q
    
18. F(1,1)=00 M  C2 [6 S: g# Q: C1 Q; Q& _4 a; w
    
19. F(1,2)=07 l; R  d9 c, c0 v8 R4 c
    
20. F(1,3)=@exp(-sin(F(1,1)))*F(1,1)
    6 ^1 Y/ @\" o- }
21. 
    4 D1 n' U( u  {( A
22. '定义龙格库塔法的权重参数
    4 m5 A( T# ]( e& z& e: W2 k+ Y
23. scalar k1
    9 H/ f7 y* w( F1 B
24. scalar k2
    ( n' \4 ?- M0 l( R) n
25. scalar k3# m% @: ?) P! j2 M\" b
    
26. scalar k4# d: y& l* U) a7 U/ z9 p
    
27. 
    1 n2 ]8 _$ G! z* b
28. '定义权重的过程量: p& i8 f* I% c2 v5 _- e
    
29. scalar w1
    2 q: s' D/ e% l  F! h; ^: r
30. scalar w22 [/ g* }' ~4 y2 n0 Y8 ^+ ]
    
31. scalar w3
    ! k  \0 M# i( i& W; Z( X
32. scalar w4
    4 m) B7 t; }* z0 f
33. # k9 L4 S4 J: M7 {
    
34. '程序主体
    ( q3 R\" A. v\" c  s
35. for !k=1 to M step 14 o5 S8 }& D+ H% f+ S: M
    
36.   F(!k+1,1)=F(!k,1)+h9 W! e& U5 W( z3 Q2 E$ T
    
37.   '调用常微分方程计算权重
    $ ~7 R0 g2 i% c\" i0 V0 i* G
38.   call obj(w1,F(!k,1),F(!k,2))
    % x\" G1 d9 E# u  _5 m! y# T7 h6 p
39.     k1=w1*h
    # L/ c! k, N$ I
40.   call obj(w2,F(!k,1)+h/2,F(!k,2)+k1/2)
    & Q( W. D# p0 X2 W
41.     k2=w2*h
    ( U3 P& b\" j6 w! E+ Q
42.   call obj(w3,F(!k,1)+h/2,F(!k,2)+k2/2)
    * d\" R% q4 ?! H- I- _
43.     k3=w3*h$ E1 V2 [8 W* }: }% ~/ W! d+ B
    
44.   call obj(w4,F(!k,1)+h,F(!k,2)+k3)
    & m; w% f7 Y9 v
45.     k4=w4*h: b2 Y* h. ~. a7 s3 `, j( {
    
46.   '计算函数估计值
    - X+ A  N. I/ E, a6 y# F\" p- ~
47.   F(!k+1,2)=F(!k,2)+(k1+2*k2+2*k3+k4)/69 Z5 c3 v) U/ p* z5 K* I$ k$ x4 x
    
48.   '计算函数解析值\" f* v( b! g. V) b
    
49.   F(!k+1,3)=@exp(-sin(F(!k+1,1)))*F(!k+1,1)
    + n9 w# A$ h( E3 o$ u  k2 B
50. next
    $ y' \1 k4 P- g7 ^, k# u' K* J
51. 3 B7 C9 w# F4 \5 `' I# _/ r! \4 l
    
52. '显示最终结果
    2 x# j1 a) p; n* I- p8 @# @
53. freeze F.xyline. M$ C; s! \7 }+ v: T+ g; q8 f
    
54. freeze F
    ; E1 C! T6 d! N( h5 H
55. 6 Z# W3 t8 A8 T2 l
    
56. '定义常微分方程
      g% e4 I! J. Y5 Y
57. subroutine obj(scalar dydx,scalar x,scalar y)' b  [% `5 N; D' G+ K5 U- G
    
58.   dydx=-y*@cos(x)+@exp(-@sin(x))
    - R$ A# e3 V$ {1 w2 J\" y& }
59. endsub7 K  ~. D. l* e7 t/ p
    

复制代码

运行后求得结果如下：

  o* [9 q5 ~0 b; j- N

) z0 E( E7 ~0 [
其中C2列是数值解，C3列是解析解，比较之下，这二者之间无明显差异。
% ?. T& j* N$ Z: m7 R" K

# V; [6 G, T: a& _% N
) D* u$ \. C& d4 r8 W8 V" F0 Y
/ j6 }( \2 c7 l( b' B& @6 [