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

liwenhui

EViews除了能解决计量经济学的估计问题以外，还提供一个编程环境用以解决复杂的问题。经过调试，我在EViews中实现了用龙格库塔方法求解常微分方程的数值解，供大家交流。
+ ?% I; x% n! v$ z/ B% m1 {; m& A演示中，我使用了如下常微分方程作为测试：
  F3 x$ B3 e0 S8 p

- g$ K9 S* r% I) f& D5 r这个方程的解析通解是：
* J: q, n  i3 o* q, k6 u

, {2 s8 \3 Q2 U/ s* i: H  U, m& C3 ?
- E7 N1 O. ^5 a2 q4 [使用“龙格库塔方法”，编制的EViews程序如下：

1. '用龙格库塔法求解常微分方程dydx=-y*cos(x)+exp(-sin(x)),y(0)=0在区间[0,10*pi]的数值解
   ! l! U5 u, a8 l2 z
2. '已知这个微分方程的解析解 y=exp(-sin(x))*x
   * |8 [# M0 L\" S7 o
3. ; F6 u4 ^9 T6 P2 p1 ]8 Z\" y
   
4. '生成一个workfile作为基本的数据容器
   ) L* |0 M& F% w0 e& X/ n4 v
5. wfcreate (wf=temp) u 1000
   $ u  a' g$ \3 F
6. 
   , d- V1 L3 U% p! O2 f7 s2 k
7. '定义常量
   * F4 b# [7 [\" V$ G& q, A6 @0 G& `
8. scalar pi=3.141595 N/ m7 T' \3 p
   
9. scalar a=0        '定义自变量下限
   - L- {) q  `, w8 g1 O
10. scalar b=10*pi     '定义自变量上限
    4 p9 ]7 g3 t7 ]: M
11. scalar  M=500       '定义步数/ R2 `3 u) y3 Q3 X\" M% |, ~
    
12. scalar h=(b-a)/M   '计算每步之间的间隔
    # T$ L8 R( F4 s# v
13. ' D- k% Q7 ~2 c* Q* U
    
14. '定义一个矩阵来储存计算数据，其中第一列储存自变量数据，第二列储存因变量数据,第三列储存解析解的值用以作为比较& F/ r9 `\" m& L5 o0 H
    
15. matrix(M+1,3) F
    8 [& y! V% ?  ]# b9 B: K' [2 J
16. % O6 J+ X' b. R5 j
    
17. '矩阵的第一行储存初值问题的初始条件; ?! o% a$ s0 W/ x( V
    
18. F(1,1)=04 Y& }( W0 W9 m\" D
    
19. F(1,2)=09 \2 {: P% @9 _9 L1 S; n3 ~
    
20. F(1,3)=@exp(-sin(F(1,1)))*F(1,1)3 @* T$ n3 e# n4 V' q8 d' M9 G
    
21. 
    1 X- ~2 e+ E+ s' n1 e0 O3 m% D
22. '定义龙格库塔法的权重参数
    / X% @1 Z! |9 H3 [- G. |
23. scalar k1; c3 P7 D8 h, h; p
    
24. scalar k2
    $ d$ X3 \4 i. w  q
25. scalar k3# p. t, v* d1 d- T
    
26. scalar k4# d& c) n: g6 `1 m* `- o  `4 t
    
27. 
    7 U& _4 w; H$ h, \( k
28. '定义权重的过程量7 |4 s2 I& K2 h\" X! w9 \. |$ C# N
    
29. scalar w1- |3 y3 H: ~: @, J$ O
    
30. scalar w29 e! S$ l* P* e) k- k6 c
    
31. scalar w3
    ) |8 x) Z* n' b, F2 e. ^1 W
32. scalar w4
    ) n2 ?/ ^2 F: v$ {5 K8 J
33. 
    ' n* X9 ^2 s* i! A& m
34. '程序主体; `: o' {/ S6 P4 h, R1 k
    
35. for !k=1 to M step 14 ~) B' [1 R, W& X\" [1 V
    
36.   F(!k+1,1)=F(!k,1)+h$ U/ e( U6 J- D& w
    
37.   '调用常微分方程计算权重( y5 l) t6 o9 L! G
    
38.   call obj(w1,F(!k,1),F(!k,2))
      A+ `, q2 q7 E0 H  {# {& f: s4 V
39.     k1=w1*h
    1 Q\" }3 I5 s\" n9 D! }
40.   call obj(w2,F(!k,1)+h/2,F(!k,2)+k1/2)8 p4 ]3 L* R) a* c+ c
    
41.     k2=w2*h
    1 u# ~; r8 [: H! _
42.   call obj(w3,F(!k,1)+h/2,F(!k,2)+k2/2)
      l: A9 B. H% v% S
43.     k3=w3*h0 k  `7 R% h' p
    
44.   call obj(w4,F(!k,1)+h,F(!k,2)+k3)3 B. X- S3 G: W# P; t9 [
    
45.     k4=w4*h
    & o: E2 E/ P8 ~: @
46.   '计算函数估计值
    1 y4 p0 R- j% n8 }1 e% f' y
47.   F(!k+1,2)=F(!k,2)+(k1+2*k2+2*k3+k4)/60 \2 ?* ?1 E6 {; `0 L1 \+ a5 o
    
48.   '计算函数解析值: t( ^- N* W  j* O, o& W2 v
    
49.   F(!k+1,3)=@exp(-sin(F(!k+1,1)))*F(!k+1,1)
    ! o2 ]& c2 F; p+ U' I
50. next' d: D$ H7 ~8 y' s2 n& Z
    
51. / f! S6 s  o' g3 J3 `, b
    
52. '显示最终结果3 a5 g; r0 ]) E\" i: w4 G; i
    
53. freeze F.xyline9 m5 m4 e$ ^2 R; a( @# z+ R
    
54. freeze F
    ( A3 Q+ p9 ^2 f& N, d
55. - ~) m; i  b1 y; i9 }8 F0 u$ r
    
56. '定义常微分方程& s4 \! Y3 u: y  _# d7 q1 K3 k
    
57. subroutine obj(scalar dydx,scalar x,scalar y)- G8 r. |$ F' e
    
58.   dydx=-y*@cos(x)+@exp(-@sin(x)). }. Y' E- v4 W0 }\" _
    
59. endsub
    ! h8 \0 r3 I  R4 S3 @( ]8 A+ I

复制代码

运行后求得结果如下：
" `! X1 t9 o; e
* h6 \" l( L% m7 X) i! r( w

4 a. \/ n' K# S' j% x其中C2列是数值解，C3列是解析解，比较之下，这二者之间无明显差异。