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

liwenhui

  u( z! W; D; k: m: N0 N
EViews除了能解决计量经济学的估计问题以外，还提供一个编程环境用以解决复杂的问题。经过调试，我在EViews中实现了用龙格库塔方法求解常微分方程的数值解，供大家交流。
演示中，我使用了如下常微分方程作为测试：

这个方程的解析通解是：
, d4 \9 R- Q% e/ m* r8 i2 V% v; I0 ]

. _6 t( m" N2 ?- t- }  s
9 ^9 R( w5 o- y& ]/ I使用“龙格库塔方法”，编制的EViews程序如下：

1. '用龙格库塔法求解常微分方程dydx=-y*cos(x)+exp(-sin(x)),y(0)=0在区间[0,10*pi]的数值解
   3 e! M5 `9 }\" j- R2 ?
2. '已知这个微分方程的解析解 y=exp(-sin(x))*x6 i- Y& ^$ }7 }# {) A
   
3. 
   1 h3 O5 M% Y! Q% Z! h5 O3 {: P8 o, x
4. '生成一个workfile作为基本的数据容器6 R0 p; I7 f- t, f& e
   
5. wfcreate (wf=temp) u 1000
   $ p. n\" r* E! p& w5 s1 n6 ?\" t
6. 
   ! j7 p6 L! O& Q  }
7. '定义常量
   ( s- O. s& @/ F2 G7 z
8. scalar pi=3.14159
   ! M! M( d4 X6 I  t. ^8 F, Y
9. scalar a=0        '定义自变量下限
   3 j\" K0 W: c% N% C
10. scalar b=10*pi     '定义自变量上限
    3 {7 E( W% n2 y
11. scalar  M=500       '定义步数
    ' L4 Q0 C3 i# h6 E
12. scalar h=(b-a)/M   '计算每步之间的间隔
    + j5 q$ I0 F\" V8 Y# @6 N  F, d% B
13. 
    $ b6 K, @* z' l, E\" X- E
14. '定义一个矩阵来储存计算数据，其中第一列储存自变量数据，第二列储存因变量数据,第三列储存解析解的值用以作为比较
    ( p0 O\" r& W7 i  G  d
15. matrix(M+1,3) F
    3 O6 S+ O, m3 Q$ U0 @  c& J0 L2 f
16. 9 \) Q$ A( @; ^0 P2 }! c
    
17. '矩阵的第一行储存初值问题的初始条件. f4 r, t  E6 w0 s. Y
    
18. F(1,1)=07 E% N! z\" ?% p7 h, t
    
19. F(1,2)=00 i' U6 n& J1 P! N6 k
    
20. F(1,3)=@exp(-sin(F(1,1)))*F(1,1)
    1 w: j% Q# F, `8 w% P\" g
21. 3 y2 r1 M: U2 s/ f; j4 A
    
22. '定义龙格库塔法的权重参数
    ) ^0 X9 ^( |$ b% E& [% I/ `
23. scalar k12 w8 H; a! @) F) C# v
    
24. scalar k2& \; o' ~. i3 l\" |+ o
    
25. scalar k35 J4 j# ?# W  c2 }! q6 Z
    
26. scalar k4
    5 C3 v! a\" q7 j, N' p
27. 
    / |: Y) i) H, P5 i3 X
28. '定义权重的过程量
    4 S0 N2 U$ ^; ]! O
29. scalar w1' \9 a8 [0 l# k8 O2 B: c\" d
    
30. scalar w2
    ; I; P; _2 C2 ?% m/ M4 y
31. scalar w3
    : a& k8 X# c# @# G  S: s
32. scalar w43 b) g- U* o) z2 ~% E4 \
    
33. ' l7 j\" A3 f5 G/ }* _
    
34. '程序主体: U2 s9 d$ v( I% Y% |
    
35. for !k=1 to M step 1
    + w/ \$ F* F. }- H\" ]
36.   F(!k+1,1)=F(!k,1)+h7 g4 i+ y* {1 t1 U, C! K1 p
    
37.   '调用常微分方程计算权重' A7 D3 a: O( H' Q
    
38.   call obj(w1,F(!k,1),F(!k,2)) 9 U) {1 k  ~6 \1 g- Z
    
39.     k1=w1*h% y# X  m! ~# Z
    
40.   call obj(w2,F(!k,1)+h/2,F(!k,2)+k1/2)+ h+ Y; m& J( d\" T. n9 t7 e- K
    
41.     k2=w2*h1 L( F0 s/ r( S9 Y5 j4 @9 F\" ?! d
    
42.   call obj(w3,F(!k,1)+h/2,F(!k,2)+k2/2)
    ' A0 R) [2 ~8 ~9 n6 ]
43.     k3=w3*h7 ?5 K7 W- u$ E+ |  R
    
44.   call obj(w4,F(!k,1)+h,F(!k,2)+k3)
    4 A) B1 B' ]/ k9 z$ z! t5 q- t( @
45.     k4=w4*h
    1 q% K/ h1 o8 a$ j4 [/ m
46.   '计算函数估计值
    ' u2 l5 T5 x1 m: {& K% C2 n
47.   F(!k+1,2)=F(!k,2)+(k1+2*k2+2*k3+k4)/6
    7 i5 F1 M( N% U6 z8 T
48.   '计算函数解析值  W+ X' |+ \5 o
    
49.   F(!k+1,3)=@exp(-sin(F(!k+1,1)))*F(!k+1,1)
    3 J( z% _* `: M. J
50. next& W3 f2 v+ m2 l- }6 [. w6 l& _3 G( l
    
51. 
    6 x3 v# I0 n: G- h
52. '显示最终结果
    7 m! n  ]/ [! p: d# K9 v
53. freeze F.xyline
    1 u0 F3 E; \& \$ h- e+ w
54. freeze F8 {: S% s2 X+ G! K0 @1 h: {! G
    
55. . T9 e9 r. ~2 Z. ^$ @  Q# q
    
56. '定义常微分方程
    6 Y. [, \5 F% x0 R, Z( j
57. subroutine obj(scalar dydx,scalar x,scalar y)/ k; A  ^9 L9 s4 E
    
58.   dydx=-y*@cos(x)+@exp(-@sin(x))' W1 \, A4 a: N% Y
    
59. endsub
    * M) A- W- ^* J& t: V

复制代码

运行后求得结果如下：

$ x) u1 f( K& O6 m7 s) T) i

其中C2列是数值解，C3列是解析解，比较之下，这二者之间无明显差异。
1 R; R* E  l4 n4 Y( G5 c, U8 o8 e
( D& m' p  Y* Y( n+ U! s% a; w( h
' e8 s8 @; K! r/ Y
+ Y# o2 P, M- v6 e$ W- p( b
" J  P8 O0 x; t7 `. q