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

liwenhui

7 k( I' w7 [6 A, s# x: D
" T* e& ~$ T8 B8 nEViews除了能解决计量经济学的估计问题以外，还提供一个编程环境用以解决复杂的问题。经过调试，我在EViews中实现了用龙格库塔方法求解常微分方程的数值解，供大家交流。
: n  K& U0 q8 K& M演示中，我使用了如下常微分方程作为测试：

2 O; n  L8 _2 s" y7 E这个方程的解析通解是：
9 G. d2 g% ^  k1 S7 Y

; O) U3 m9 I( }/ D+ }
' z2 T3 D9 V. a使用“龙格库塔方法”，编制的EViews程序如下：

1. '用龙格库塔法求解常微分方程dydx=-y*cos(x)+exp(-sin(x)),y(0)=0在区间[0,10*pi]的数值解2 _# r+ p- B: e2 G2 p
   
2. '已知这个微分方程的解析解 y=exp(-sin(x))*x
   , Z  Q- |4 H; z3 i) x3 E0 ]: H1 ~
3. 
   ( ~+ d( ^  t4 J4 [- @+ q/ ?  O, V
4. '生成一个workfile作为基本的数据容器
   0 A- i' Z/ r: `
5. wfcreate (wf=temp) u 10004 P7 H# C: C' e5 J7 N- ~, ~0 ^
   
6. 
   , G( o: k4 X4 W* O6 z6 D! l
7. '定义常量
   + a& R4 P' _. S0 z\" l' K
8. scalar pi=3.14159
   ; F# T; [  F8 s& z/ F# e/ V0 v* v7 W
9. scalar a=0        '定义自变量下限2 Z; g9 l3 _8 `\" w. n
   
10. scalar b=10*pi     '定义自变量上限$ }4 k: Z; _) U; h' B3 P
    
11. scalar  M=500       '定义步数7 k! |\" P) Z; U7 `; b! x# I
    
12. scalar h=(b-a)/M   '计算每步之间的间隔
    , ^* n. C/ Y, V6 I$ ?
13. 
    & a! L) }* U: Q0 P9 Q$ `$ m* e& R! R
14. '定义一个矩阵来储存计算数据，其中第一列储存自变量数据，第二列储存因变量数据,第三列储存解析解的值用以作为比较* K9 l- G- v! G+ \* k0 O$ C8 o! U3 r  L
    
15. matrix(M+1,3) F % d: J  h. R. Z
    
16. 
    5 d' e\" ~; X/ N7 a
17. '矩阵的第一行储存初值问题的初始条件
    & p  m& R! ^' K1 O\" A. P; ~
18. F(1,1)=0- P% p; S0 C3 h6 o. W6 G6 w8 x* X
    
19. F(1,2)=0$ U7 v  D+ t9 ^- ^3 _
    
20. F(1,3)=@exp(-sin(F(1,1)))*F(1,1)- q5 L& K4 L; r7 d6 j+ p& b6 G
    
21. 4 h; D; o+ A' i% q3 y
    
22. '定义龙格库塔法的权重参数
    - K% L, _, U! N7 M5 h+ R) U7 Q
23. scalar k1$ H- z, T& L& O3 l4 h
    
24. scalar k2* n: Q3 ]5 q% h  g
    
25. scalar k3
    7 L9 x7 Z; \/ _, l
26. scalar k4
    7 j1 M' u) e# o$ X2 g8 J- j
27. 0 z0 `) h' o, b9 G. i1 v
    
28. '定义权重的过程量. C# e+ K& ]+ {6 c, D7 j
    
29. scalar w1+ i, d$ u$ M7 Q6 S4 @# ]! |* M
    
30. scalar w2
    8 G/ z, Q; u+ i6 T( D1 a8 t; T
31. scalar w31 m' I# Z( V- K6 Y) ?
    
32. scalar w4, }$ W2 \, G2 i\" h
    
33. 
    6 H! V7 O3 E1 Q1 r/ v; [4 U\" F
34. '程序主体2 W8 G7 i0 O2 P7 k
    
35. for !k=1 to M step 1
    3 d. ?: ]\" U+ B: ~% `! G2 \. r- l
36.   F(!k+1,1)=F(!k,1)+h
    1 |- B) P+ S( J. J$ G2 U4 C, r! a) k
37.   '调用常微分方程计算权重\" w7 b\" a) f! r4 E3 r
    
38.   call obj(w1,F(!k,1),F(!k,2))
    \" c! R; c1 R) k
39.     k1=w1*h
    # a6 Y: D( o/ J# i2 I1 D. u
40.   call obj(w2,F(!k,1)+h/2,F(!k,2)+k1/2)
    + l- k, r. ^. V
41.     k2=w2*h
    ; N  X9 U9 q! w; x; w+ {& h
42.   call obj(w3,F(!k,1)+h/2,F(!k,2)+k2/2)( h1 O$ G) Y2 J; q
    
43.     k3=w3*h9 S$ R- _5 c' ?
    
44.   call obj(w4,F(!k,1)+h,F(!k,2)+k3)2 Q& q3 V0 g9 v5 U! [
    
45.     k4=w4*h
    7 K2 b. F; k1 r6 S4 a0 O
46.   '计算函数估计值; O6 z: l0 D5 p1 |8 S
    
47.   F(!k+1,2)=F(!k,2)+(k1+2*k2+2*k3+k4)/6! \# `6 J) k3 K/ \- B
    
48.   '计算函数解析值  W- q1 `1 N- c( w/ Q' B+ l/ I
    
49.   F(!k+1,3)=@exp(-sin(F(!k+1,1)))*F(!k+1,1)
    ! Z4 Q' l6 D3 r3 z# b
50. next& Y8 g6 @$ X# y0 {' \4 l7 s
    
51. - C, L! r1 A; f7 L/ E
    
52. '显示最终结果1 L; }: q) e8 e$ L- ^. M! x. n9 J
    
53. freeze F.xyline
    ( X4 i2 \4 N' v
54. freeze F
    / e* F; e, e9 m5 s) e1 l' Z0 T0 ]7 ^\" B
55. % K' \, n% a6 {9 y* [
    
56. '定义常微分方程
    8 y/ w4 h  j9 m) m  E
57. subroutine obj(scalar dydx,scalar x,scalar y)
    - Q& P; W\" w& u
58.   dydx=-y*@cos(x)+@exp(-@sin(x))/ p! S+ }8 ~4 B! C2 d1 ]
    
59. endsub
    \" J) Z/ d- v% S. {' e8 U# x

复制代码

运行后求得结果如下：
. Y2 D- k0 n0 x( k& k

  W: T7 ~- f" q& m' x9 G
9 M: T& Z  r- c& U3 `其中C2列是数值解，C3列是解析解，比较之下，这二者之间无明显差异。
! E9 c0 M2 m4 U: J
6 R) l! {4 |1 E. u! K' F
' q6 U; a5 c- U. c4 b) K7 _8 T