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

liwenhui

4 W4 e9 u6 X$ \0 o9 J: b6 w
EViews除了能解决计量经济学的估计问题以外，还提供一个编程环境用以解决复杂的问题。经过调试，我在EViews中实现了用龙格库塔方法求解常微分方程的数值解，供大家交流。
演示中，我使用了如下常微分方程作为测试：
, }$ f2 B% ^/ w  D1 r! L6 n

这个方程的解析通解是：
5 {& ]4 _0 w* I2 g6 w9 W  _' x% l+ {

9 {  ~) V, Q, y( S: r
' c+ n+ z" f1 h/ W3 V. v6 D% _, G' G5 e使用“龙格库塔方法”，编制的EViews程序如下：

1. '用龙格库塔法求解常微分方程dydx=-y*cos(x)+exp(-sin(x)),y(0)=0在区间[0,10*pi]的数值解
   5 K- q$ f0 n8 m7 L
2. '已知这个微分方程的解析解 y=exp(-sin(x))*x# Y+ ?0 z$ R) A$ f, n9 q
   
3. 
   4 p& c- N8 n7 w! \6 u0 E
4. '生成一个workfile作为基本的数据容器
   1 R0 q1 z- E+ p+ r
5. wfcreate (wf=temp) u 1000
   : b8 ?5 ?% z' ]- W; u  r+ l
6. 
   ) |, r& r6 f( \: t, X
7. '定义常量5 `) `, ?+ q0 W  D/ ~' C. J
   
8. scalar pi=3.14159
   : G4 u0 P7 S* g3 L: ~
9. scalar a=0        '定义自变量下限2 [+ J6 p& p* Y) u
   
10. scalar b=10*pi     '定义自变量上限
    7 v& A5 \. l- O& h! _6 S
11. scalar  M=500       '定义步数
    8 L5 ]& x  |/ \) s6 X) V- ?9 o
12. scalar h=(b-a)/M   '计算每步之间的间隔
    9 U4 O$ Q: ^- f; I) U' p
13. 
    ) m0 u+ }+ o) J  i
14. '定义一个矩阵来储存计算数据，其中第一列储存自变量数据，第二列储存因变量数据,第三列储存解析解的值用以作为比较8 l- t# ?\" |4 l8 K0 ^1 h8 o
    
15. matrix(M+1,3) F
    9 }( J$ C: F& e8 f; ?, h9 }
16. 
    ; i# D8 t$ n/ z8 r7 U* q
17. '矩阵的第一行储存初值问题的初始条件
    5 T+ z\" W# s6 c( K; |. E( O, `
18. F(1,1)=0
    ; w$ ~! A$ ^5 S* \
19. F(1,2)=0* _$ U5 g) w& \, B\" `
    
20. F(1,3)=@exp(-sin(F(1,1)))*F(1,1)
    & P3 G. N7 F: E
21. ( J/ \; q( P! l
    
22. '定义龙格库塔法的权重参数\" j, r. G# \1 m! H# q: `$ O
    
23. scalar k1% B, |) `* Z  [7 c  n1 l
    
24. scalar k2
    0 n$ b7 O- @& o
25. scalar k3
    & R7 U3 A* Z0 h\" w- B
26. scalar k42 d/ H6 w, z\" t0 `. t1 m/ q
    
27. 6 y3 i- [, B- Q0 b& j5 y2 \
    
28. '定义权重的过程量
    8 v9 z) X, p5 z4 Z' y
29. scalar w1
    $ S, M- j  N. c# D
30. scalar w24 P, [7 q7 J9 R5 v' V# G
    
31. scalar w3
    & @8 Z; R3 [4 ?3 @+ W2 \
32. scalar w4
      `' P2 |$ E$ x' l7 z- B
33. 
    # R: D+ `) a: r, ~, c
34. '程序主体/ ?  I+ q$ ?$ ~* C1 r0 m
    
35. for !k=1 to M step 1' t6 r- \; G/ Z- ^
    
36.   F(!k+1,1)=F(!k,1)+h
      S- t1 o6 _; M+ [3 o# i3 n: E$ T
37.   '调用常微分方程计算权重# m\" n6 X' e+ B
    
38.   call obj(w1,F(!k,1),F(!k,2))
    ! t% u1 i4 q4 B5 B' x! w! z- c
39.     k1=w1*h
    . X1 J$ n! k! S6 o& g) [: [& f\" ~1 [
40.   call obj(w2,F(!k,1)+h/2,F(!k,2)+k1/2)0 f& e+ w% B( b7 C  o+ R& W
    
41.     k2=w2*h\" e  D0 A% J1 K
    
42.   call obj(w3,F(!k,1)+h/2,F(!k,2)+k2/2)% L+ t7 y# n& M
    
43.     k3=w3*h
    . W/ D. }, \- f% N
44.   call obj(w4,F(!k,1)+h,F(!k,2)+k3)  D3 I/ b) B) E2 C9 P
    
45.     k4=w4*h5 A  X2 V5 [+ s% l' a+ m- E
    
46.   '计算函数估计值
    & M5 \+ s) q2 f  d9 w
47.   F(!k+1,2)=F(!k,2)+(k1+2*k2+2*k3+k4)/6
      P9 K% F1 i2 I+ q# T! e
48.   '计算函数解析值+ y& I' w6 W6 S2 z% Z- t# q; ~\" V
    
49.   F(!k+1,3)=@exp(-sin(F(!k+1,1)))*F(!k+1,1)) C! Q/ g0 N% w. q\" G
    
50. next
    # _/ X* k/ M) a4 B$ `) ]
51. 
    6 U) c& @6 A; j
52. '显示最终结果
      E( U- T- z; Z, {6 J( y7 ?
53. freeze F.xyline
    6 B- N; {: |! \4 r( F
54. freeze F9 Y8 B; l  m4 F7 E2 @  a
    
55. / w. H6 I/ z6 q2 v7 n
    
56. '定义常微分方程$ M2 z. `0 O\" b  r/ Y3 @# d. t
    
57. subroutine obj(scalar dydx,scalar x,scalar y)0 v0 e6 I. C; _  O* p' K( k; B$ e
    
58.   dydx=-y*@cos(x)+@exp(-@sin(x))\" b* e. f0 f- b
    
59. endsub) l- d. [4 d7 O\" K/ S2 H
    

复制代码

运行后求得结果如下：
1 c( q. |( I- N8 C6 l
1 B8 {. [1 }7 n8 c- y; J

其中C2列是数值解，C3列是解析解，比较之下，这二者之间无明显差异。


5 g+ j, U4 X( X6 s1 r