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

liwenhui

9 F4 C: G4 I; w, H' i. t! _! `" D4 X
EViews除了能解决计量经济学的估计问题以外，还提供一个编程环境用以解决复杂的问题。经过调试，我在EViews中实现了用龙格库塔方法求解常微分方程的数值解，供大家交流。
演示中，我使用了如下常微分方程作为测试：

这个方程的解析通解是：

$ w  \6 L7 M3 O8 H4 Q1 u( P
# r0 v( o: x$ K; |" C使用“龙格库塔方法”，编制的EViews程序如下：

1. '用龙格库塔法求解常微分方程dydx=-y*cos(x)+exp(-sin(x)),y(0)=0在区间[0,10*pi]的数值解& d7 N) J& Q' R# M6 `3 ~! {3 y
   
2. '已知这个微分方程的解析解 y=exp(-sin(x))*x4 U. A% p3 M3 E. Q: ]: D2 w
   
3. 
   6 Z+ v2 ?  w, Y# R% m1 J\" N
4. '生成一个workfile作为基本的数据容器
   . S; O4 Q, o& y' L
5. wfcreate (wf=temp) u 1000( N8 ^* h% p, o3 ~1 N, k
   
6. 
   / p' X: x6 ^4 C; h' k8 p\" I3 ?
7. '定义常量
   3 d& [3 i/ c) ~) R1 m
8. scalar pi=3.14159
   ) B1 ~# G' a! q( a
9. scalar a=0        '定义自变量下限6 Q5 l8 }' Y3 x$ Z  G! R
   
10. scalar b=10*pi     '定义自变量上限
    6 s8 u: k. ]# ~\" V
11. scalar  M=500       '定义步数
    * f& S! n1 F  c7 U; N2 X2 {. j
12. scalar h=(b-a)/M   '计算每步之间的间隔5 r6 n: n; D7 Y$ b
    
13. 8 }' J\" {\" {- _& v4 P2 g6 U
    
14. '定义一个矩阵来储存计算数据，其中第一列储存自变量数据，第二列储存因变量数据,第三列储存解析解的值用以作为比较% n6 J* L6 T( I+ e
    
15. matrix(M+1,3) F . C: i$ D. Q6 u8 A( _5 r! V
    
16. 
    4 g1 c3 B* m6 f) E
17. '矩阵的第一行储存初值问题的初始条件
    ) V7 ]2 u. a4 s( M
18. F(1,1)=0% Y0 c9 o5 D( g& B\" V0 N
    
19. F(1,2)=06 c- @/ A* i) m. O; X8 c. k
    
20. F(1,3)=@exp(-sin(F(1,1)))*F(1,1)
    # [- C1 G( c- F' q6 N- i3 _  B
21. / X7 h9 T7 P6 E
    
22. '定义龙格库塔法的权重参数
    2 }% j7 Z+ e+ G+ I' z' S5 b, }
23. scalar k1# [5 B* A% x9 p6 b! D: b5 R) a
    
24. scalar k2\" k: E  V! k9 w: r9 g& l
    
25. scalar k3  z; B' r+ `& B8 Q+ e
    
26. scalar k4
    4 ?2 N- Q; [5 a
27. * g& @6 I( {0 j# Z: Y! s3 C
    
28. '定义权重的过程量
    \" W6 ?, h$ J6 x6 n% [4 C
29. scalar w1, D\" ~\" r* V9 c5 K\" y4 ]
    
30. scalar w2
    1 e& J! Y. l: M1 d8 {1 V( S8 ?( P
31. scalar w3
    * t% x: S$ [4 a# V8 {. ~
32. scalar w4
    2 K/ b0 Y% ]( D! D% p( d4 A6 Q$ e
33. 
    # W0 ^- A! Y5 ?( u3 F
34. '程序主体
    0 ]- X$ x! r3 b7 R' t2 E8 K: }: b
35. for !k=1 to M step 17 v: r$ m$ v+ H8 f1 c\" n
    
36.   F(!k+1,1)=F(!k,1)+h
    . V% X4 P1 Z' `2 u
37.   '调用常微分方程计算权重$ l' x, H; n* z0 Y9 r\" A2 C
    
38.   call obj(w1,F(!k,1),F(!k,2))   }$ L; z& |! d( V& F- j( |
    
39.     k1=w1*h! T. V) r) Z1 P6 [' ^5 e
    
40.   call obj(w2,F(!k,1)+h/2,F(!k,2)+k1/2)+ C! W4 T% G: v# {& C0 @$ `6 d
    
41.     k2=w2*h
    + J( Z& L7 t3 G. e' E$ f; B+ h  V
42.   call obj(w3,F(!k,1)+h/2,F(!k,2)+k2/2)
    1 ^\" _* {% K\" Z* N. S) {4 A; `
43.     k3=w3*h, {: g+ d4 B2 {$ G& r! j
    
44.   call obj(w4,F(!k,1)+h,F(!k,2)+k3)
    * Z9 m' g! a# L3 p\" n+ [
45.     k4=w4*h# c$ o( H; T, g; i0 S: _7 M# x
    
46.   '计算函数估计值
    6 u# X! L! @0 ^  D* E
47.   F(!k+1,2)=F(!k,2)+(k1+2*k2+2*k3+k4)/62 v' \0 k/ V4 C) E
    
48.   '计算函数解析值5 i7 u8 A1 _( b1 k
    
49.   F(!k+1,3)=@exp(-sin(F(!k+1,1)))*F(!k+1,1)8 m! j8 T2 `\" J$ {+ I
    
50. next6 z1 D( I8 j, [' i, r+ i
    
51. 5 s# a( n3 [. y. i; P8 w
    
52. '显示最终结果# _2 i) u% {* ~$ c
    
53. freeze F.xyline
    ; q. ?' [8 x9 ^- X
54. freeze F$ m8 p  \2 V4 U; Y' H
    
55. 
    & p$ w\" m9 t4 U. K\" F- r1 O
56. '定义常微分方程5 _+ |' H3 c2 w3 K  H9 I
    
57. subroutine obj(scalar dydx,scalar x,scalar y)
    ! a4 D8 f9 Z' H! C
58.   dydx=-y*@cos(x)+@exp(-@sin(x))! D5 M2 t  S8 B/ \/ ^! _8 K1 r
    
59. endsub
    ( q\" j! _: r( d8 M$ ]

复制代码

运行后求得结果如下：

$ i$ C" l# t- h4 O; `' q6 a) _0 u
其中C2列是数值解，C3列是解析解，比较之下，这二者之间无明显差异。


  L4 N7 p+ b" ~/ x  N/ O
" D, R7 P5 l0 C7 e  w
& K% g" q+ n9 ~& ~  ~0 ^; e