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

liwenhui

4 `/ @; Y( ]7 ~" B! y
EViews除了能解决计量经济学的估计问题以外，还提供一个编程环境用以解决复杂的问题。经过调试，我在EViews中实现了用龙格库塔方法求解常微分方程的数值解，供大家交流。
演示中，我使用了如下常微分方程作为测试：
/ v+ Q# \$ k2 I8 k! ^

* r8 E4 E" P1 v2 A3 K7 C这个方程的解析通解是：

使用“龙格库塔方法”，编制的EViews程序如下：

1. '用龙格库塔法求解常微分方程dydx=-y*cos(x)+exp(-sin(x)),y(0)=0在区间[0,10*pi]的数值解& M) i7 `# ]! H; D
   
2. '已知这个微分方程的解析解 y=exp(-sin(x))*x& Y% v/ S* g1 E0 ?! ^* f4 L$ H5 B) k
   
3. 
   5 K4 [! D. N, e7 W  \4 R
4. '生成一个workfile作为基本的数据容器9 G/ R  R( w1 K. h0 c+ B
   
5. wfcreate (wf=temp) u 1000& q: ?' q\" U: s
   
6. 
   ! t% K2 y( ?: _8 F' ~
7. '定义常量
   0 j& v6 `( W  M, u% _1 D
8. scalar pi=3.14159
   6 f/ O8 ?\" h9 Y7 M4 p
9. scalar a=0        '定义自变量下限
   6 i\" A. m& n) E+ r, Y+ L
10. scalar b=10*pi     '定义自变量上限0 j* c0 g. H0 P: ?0 F* U+ w
    
11. scalar  M=500       '定义步数
    3 |  l' f) ]* X  Q/ a
12. scalar h=(b-a)/M   '计算每步之间的间隔
    ! h$ d0 {( h& a' v, }0 l
13. # R' @, u, W% `4 w
    
14. '定义一个矩阵来储存计算数据，其中第一列储存自变量数据，第二列储存因变量数据,第三列储存解析解的值用以作为比较
    8 H8 w2 h% ]\" s
15. matrix(M+1,3) F
    # H# F, ^% _9 W6 X9 d1 W
16. / n+ l7 [3 ]! f
    
17. '矩阵的第一行储存初值问题的初始条件' Q6 O* Q/ M, h- u
    
18. F(1,1)=05 r1 K) h6 j+ _1 C4 ]% H4 ^: \
    
19. F(1,2)=09 R# b7 t8 t) D/ ]
    
20. F(1,3)=@exp(-sin(F(1,1)))*F(1,1)
    8 f6 ]% k# i2 I, D5 h
21.   z/ i\" P& ^* Y: P. Z$ Y' [
    
22. '定义龙格库塔法的权重参数
    $ g$ H1 x5 N8 g
23. scalar k1
    % y3 N$ Q' c( Y) Q' Z\" U2 G
24. scalar k28 Y2 ~\" H* d' u2 i/ R/ C
    
25. scalar k31 m7 B: a7 l. M1 \4 }
    
26. scalar k4( A\" g. J( k6 D/ n* ?7 p
    
27. 
    8 t\" m) w2 l/ @
28. '定义权重的过程量
    2 H7 L3 z; B2 U! S5 [! g
29. scalar w15 t! l' Y, g3 M\" k4 ~\" N9 p
    
30. scalar w2
    ; r+ G, N$ z  \. A- z1 j# ^
31. scalar w3- z1 r* b2 G* K% @
    
32. scalar w4
      f  G3 k. X* e: `1 ]- T9 k% S+ a) l' V, Z
33. 1 f% O3 c\" M7 j+ `+ K; j# z/ }
    
34. '程序主体
    : X8 u, B: ^# y  z7 [2 A
35. for !k=1 to M step 1\" u/ e9 y; T$ d8 n0 e/ T# g
    
36.   F(!k+1,1)=F(!k,1)+h\" _& f4 Z% Q& i+ D4 i* t
    
37.   '调用常微分方程计算权重6 f' j6 C$ N0 a
    
38.   call obj(w1,F(!k,1),F(!k,2))
    $ h: f5 _$ }& u! t- t! I7 K
39.     k1=w1*h
    / k1 A/ m( I: o& X' B+ s& }4 P2 \
40.   call obj(w2,F(!k,1)+h/2,F(!k,2)+k1/2)
    1 k# f\" n( J1 Z$ a# f0 ?' k( B
41.     k2=w2*h! Z% [0 C: @$ k
    
42.   call obj(w3,F(!k,1)+h/2,F(!k,2)+k2/2)6 X, V9 I3 q5 W% A( p- u0 f
    
43.     k3=w3*h
    ' w2 ^7 X. Q# n1 `
44.   call obj(w4,F(!k,1)+h,F(!k,2)+k3)3 b# ]: X\" O% J( Z( J8 e2 r' @3 {3 ?
    
45.     k4=w4*h
    ! K* K$ O* ]+ \
46.   '计算函数估计值, c( S, i0 s/ E: f8 J# @( f# N
    
47.   F(!k+1,2)=F(!k,2)+(k1+2*k2+2*k3+k4)/6. s\" j9 q. o' K
    
48.   '计算函数解析值' D% @- ~- g6 c* U1 X% Y/ P1 C, M
    
49.   F(!k+1,3)=@exp(-sin(F(!k+1,1)))*F(!k+1,1)
    7 b. {5 L6 }  ~) P\" O! i, v! S
50. next
    9 C6 h2 `' E! ^( f% Q* w
51. 
    ( q; A/ _\" r. Y/ @  m$ t
52. '显示最终结果4 W: S* R0 R, e
    
53. freeze F.xyline
    + F3 c% i6 w( E% Z\" x
54. freeze F
    % Z\" u2 u& J* u. p/ e1 W
55. % e+ E7 Y. t' i\" o% x6 A; d- I
    
56. '定义常微分方程
    1 T* Z5 c7 Z2 Z  U\" @. U
57. subroutine obj(scalar dydx,scalar x,scalar y)8 l7 M) C2 b3 Y* k3 \
    
58.   dydx=-y*@cos(x)+@exp(-@sin(x))
    / [! ^7 P, X; z% W/ k6 R/ n
59. endsub
    8 G* H1 O$ F7 D\" f

复制代码

运行后求得结果如下：

  r% p1 ~0 b; ~7 c5 ?( K  p& r/ O

其中C2列是数值解，C3列是解析解，比较之下，这二者之间无明显差异。
2 ^3 ~) f* E2 P, E
8 P# P) w% m. U0 C7 V


" I: ~: c2 @0 K, D* P3 b/ f