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

liwenhui

EViews除了能解决计量经济学的估计问题以外，还提供一个编程环境用以解决复杂的问题。经过调试，我在EViews中实现了用龙格库塔方法求解常微分方程的数值解，供大家交流。
% i+ N# R; F* {& S, e' J, Y# |演示中，我使用了如下常微分方程作为测试：

这个方程的解析通解是：

5 X) _9 b/ f3 s3 ?& x7 S( g
使用“龙格库塔方法”，编制的EViews程序如下：

1. '用龙格库塔法求解常微分方程dydx=-y*cos(x)+exp(-sin(x)),y(0)=0在区间[0,10*pi]的数值解
   ' q3 u4 j, k5 s% w8 L
2. '已知这个微分方程的解析解 y=exp(-sin(x))*x
   2 J$ n& {( e( y3 s8 P
3. 7 W3 J. Y* z' a! U; s) N
   
4. '生成一个workfile作为基本的数据容器
   7 P# d9 H: |) }- R
5. wfcreate (wf=temp) u 1000
   % N+ A$ Q2 p) g4 J1 h7 N
6. 
   - Q5 g& m. |4 G2 F% \  t
7. '定义常量' m) e' H0 `5 v7 U  W
   
8. scalar pi=3.141595 m/ E) R# w4 u# p0 t; m
   
9. scalar a=0        '定义自变量下限\" L3 s' c3 e& q2 T# X* v
   
10. scalar b=10*pi     '定义自变量上限
    7 D3 c: D! F1 t9 i! C# p\" p3 q
11. scalar  M=500       '定义步数
    ' E* M# C) N7 m1 m9 R% V
12. scalar h=(b-a)/M   '计算每步之间的间隔
    4 c9 n  b4 Q( l: o( M0 N
13. $ W4 Y7 |( a/ n
    
14. '定义一个矩阵来储存计算数据，其中第一列储存自变量数据，第二列储存因变量数据,第三列储存解析解的值用以作为比较
    ' ^4 W, U: U: [\" u0 n
15. matrix(M+1,3) F
    ) B# H4 k) ?( S( v+ t  _' s
16. 
    ! g& P3 J) q/ x% h' E
17. '矩阵的第一行储存初值问题的初始条件/ ~9 I* z! Z6 ?6 h8 Z8 F  ?; k
    
18. F(1,1)=09 _' P& f1 b9 h0 C9 u; \7 |
    
19. F(1,2)=0
    ' b1 e& q& l+ g  [
20. F(1,3)=@exp(-sin(F(1,1)))*F(1,1)- x, p, E\" m! m0 O2 ]
    
21. 
    . Y  w4 `/ L% g\" h; a/ A& c
22. '定义龙格库塔法的权重参数
    2 h6 b6 s7 y5 W% Z6 R
23. scalar k1! r8 \3 l# U3 [5 B/ }9 M, M( @
    
24. scalar k2: T6 a3 O1 r* z\" {* D
    
25. scalar k31 Q  ?7 _7 O\" x) I
    
26. scalar k4
    & a+ v0 F- N# Y. U
27. 
    9 h5 m+ V: B* i- A$ J( I1 o
28. '定义权重的过程量# k, i: u: l+ M( Z2 V
    
29. scalar w1
    \" [1 z\" N, q* u; g/ a0 s% C  I
30. scalar w2* B) |- d* Q3 `$ x6 Y6 m
    
31. scalar w3
    . }  v# E4 N: o7 A\" v
32. scalar w4
    : w3 Q6 h) z* x/ u
33. ! I3 ]; `; G( V, q
    
34. '程序主体& p$ c5 ~: ?( d% M- ]
    
35. for !k=1 to M step 1
    \" x' ?2 a; x\" j( D
36.   F(!k+1,1)=F(!k,1)+h, ^6 t# A+ v& o9 {\" F
    
37.   '调用常微分方程计算权重
    ( N' V; N4 L5 `1 i( Q: P/ [
38.   call obj(w1,F(!k,1),F(!k,2))
    1 F5 U2 S0 u\" [
39.     k1=w1*h
    ' a: y1 Y' Q! H8 ?
40.   call obj(w2,F(!k,1)+h/2,F(!k,2)+k1/2)
    $ t! d7 g9 Z0 `
41.     k2=w2*h
    , ?6 Z6 V! W\" O8 l
42.   call obj(w3,F(!k,1)+h/2,F(!k,2)+k2/2)
    0 `: D/ V# m5 b\" ]$ e
43.     k3=w3*h( X' e1 ~! L0 }' a* `# w
    
44.   call obj(w4,F(!k,1)+h,F(!k,2)+k3)
    $ e3 v* D2 d. J8 L  A: \/ e& k
45.     k4=w4*h0 R4 F2 G9 \( j5 |6 k
    
46.   '计算函数估计值
    7 r8 N3 v1 e( h6 t$ ^& p
47.   F(!k+1,2)=F(!k,2)+(k1+2*k2+2*k3+k4)/6
    ! x7 p# t1 K- z6 e% U/ e
48.   '计算函数解析值# u5 o2 `: ~: l, g1 x! A' {2 C
    
49.   F(!k+1,3)=@exp(-sin(F(!k+1,1)))*F(!k+1,1)
    - f) A- |1 ~& y9 ^/ [+ ~
50. next/ S  p( \4 W\" T# {) K
    
51. 
    8 z7 B! m; d* G/ B  C
52. '显示最终结果
    7 X1 @( R* z6 @3 s- @
53. freeze F.xyline
    7 ^& e' \\" Z$ {9 Q) h\" C
54. freeze F- i' B4 L\" }1 e# s/ ?# J. ]
    
55. . F; g# n5 d0 ?# Y\" I
    
56. '定义常微分方程
    4 o1 N7 f; a\" `# K9 Z
57. subroutine obj(scalar dydx,scalar x,scalar y)( Y; g+ s1 |% _& I
    
58.   dydx=-y*@cos(x)+@exp(-@sin(x))
    + T% k% m  v* R) ~! y
59. endsub9 O& R- U1 S: N
    

复制代码

运行后求得结果如下：

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

+ `( X8 o7 g3 C' W% t& x1 P0 k8 I' ~" d