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

liwenhui

0 B# X* A2 K; ?, a8 q2 AEViews除了能解决计量经济学的估计问题以外，还提供一个编程环境用以解决复杂的问题。经过调试，我在EViews中实现了用龙格库塔方法求解常微分方程的数值解，供大家交流。
  ]! R  B" Q; H1 D/ u. v1 R演示中，我使用了如下常微分方程作为测试：
* Y1 j% @7 }4 o& ]* y! P

% X6 Y" d1 B5 g$ z
3 W) S9 C8 Q$ Z3 e0 Q) e) l这个方程的解析通解是：
' P# ^9 t5 h4 e' r& Z

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

1. '用龙格库塔法求解常微分方程dydx=-y*cos(x)+exp(-sin(x)),y(0)=0在区间[0,10*pi]的数值解4 A  A. N5 W6 c
   
2. '已知这个微分方程的解析解 y=exp(-sin(x))*x
   & H% m2 @+ Z) l+ D* H* k
3. 
   ! f0 \. W4 F* t( i
4. '生成一个workfile作为基本的数据容器5 @- {7 ]5 R, B# [. G; i
   
5. wfcreate (wf=temp) u 1000
   / J! E+ Z  a3 d  n' z! g) G6 q7 e
6. 
   + b3 J. R& G% E: D: g
7. '定义常量
   1 @2 L: t\" |4 U, j% G
8. scalar pi=3.14159
   % [\" W) U+ E4 P, C8 J
9. scalar a=0        '定义自变量下限\" T  Q; U$ _: E2 ?4 v# Z! r% p# I
   
10. scalar b=10*pi     '定义自变量上限
    6 U  x+ A& B0 v# V5 X
11. scalar  M=500       '定义步数( Y# h4 g4 s- d) M7 L9 E
    
12. scalar h=(b-a)/M   '计算每步之间的间隔
    / S$ ]* n4 \- c5 b# K. i% T
13. 
    \" g6 G$ a/ y  ]3 m
14. '定义一个矩阵来储存计算数据，其中第一列储存自变量数据，第二列储存因变量数据,第三列储存解析解的值用以作为比较
    6 K4 c- f' Q' `* M
15. matrix(M+1,3) F
    6 X) }, T# ^5 l( }- o/ C+ c6 a
16. ; x. |( Q& h& N5 Q8 l
    
17. '矩阵的第一行储存初值问题的初始条件
      }4 g\" ~/ B( E3 S6 F! ^
18. F(1,1)=0/ P2 O+ R; I2 n
    
19. F(1,2)=0! z4 I\" B\" t$ X7 K
    
20. F(1,3)=@exp(-sin(F(1,1)))*F(1,1)' K\" i- X7 I' y1 h8 p
    
21. + Z* ?: i/ h9 l+ D2 ]) S
    
22. '定义龙格库塔法的权重参数
    8 T1 ]) ~( G7 }9 Q# W! M! X
23. scalar k1
    6 y% G2 v5 \) i/ d\" i# T
24. scalar k2
      d, R7 Q3 T4 R
25. scalar k3
    / b( B! |5 o6 \3 p- @: \$ E8 S: B4 l
26. scalar k4
    1 w* z+ o' L6 \6 ]( e) P
27. 
      r; W: q- e: D. C
28. '定义权重的过程量' I7 {\" m) v( }$ @* U& t# J  _
    
29. scalar w1
    ( J2 M( z* V6 @; `+ D
30. scalar w2
    % H+ p4 }: Z; o1 ?1 J7 L1 g
31. scalar w3- v3 g+ U2 W  m4 Y; O' ~' ~
    
32. scalar w4
    , e* i9 J) P- _/ X2 Z4 S0 m9 |\" u2 i' V
33. 
    . [' H, h+ i3 _1 J' b: K
34. '程序主体  P7 k. s7 J6 ~# N  f! B
    
35. for !k=1 to M step 15 i: ?) {5 {5 e' n
    
36.   F(!k+1,1)=F(!k,1)+h' s5 t, Q( J7 _, ~
    
37.   '调用常微分方程计算权重9 ?6 I\" f7 y\" y, h  ]; a$ G
    
38.   call obj(w1,F(!k,1),F(!k,2))
    # u+ s5 j2 u# `# A5 n
39.     k1=w1*h
    9 \# E7 ^% e& d6 o& |8 C7 m! w! ^: I
40.   call obj(w2,F(!k,1)+h/2,F(!k,2)+k1/2)* r, n& g9 u4 J6 p\" ~
    
41.     k2=w2*h# ^2 g) B\" s& U6 V% Q
    
42.   call obj(w3,F(!k,1)+h/2,F(!k,2)+k2/2)$ m: {6 [9 `! s+ [
    
43.     k3=w3*h9 B9 R% P' i6 }$ g
    
44.   call obj(w4,F(!k,1)+h,F(!k,2)+k3)5 p; x/ p# g) n3 {5 ^5 P: z( w
    
45.     k4=w4*h
    0 m) F6 h% z5 w% U# n
46.   '计算函数估计值
    8 Z6 F; ]6 r. _; b; S
47.   F(!k+1,2)=F(!k,2)+(k1+2*k2+2*k3+k4)/6% b% N) Z3 p$ A\" \- c6 t% h! m
    
48.   '计算函数解析值
    1 H/ |6 l, L5 u5 e
49.   F(!k+1,3)=@exp(-sin(F(!k+1,1)))*F(!k+1,1)
    , k0 C: X' X$ A7 H$ ~
50. next4 X+ ?( P7 E8 v
    
51. 
    ! @3 v! s( v8 p% e
52. '显示最终结果' D+ q8 K3 W/ v9 \& o
    
53. freeze F.xyline, I$ L' G: L& G0 ?% M\" W
    
54. freeze F
    $ u* C( j% }. P3 {5 `
55. 
    ! D' B& @8 }2 i0 W  ~# S7 v
56. '定义常微分方程
    * q) ?' L& ^) q
57. subroutine obj(scalar dydx,scalar x,scalar y)
    , {7 m7 p\" g+ W' [
58.   dydx=-y*@cos(x)+@exp(-@sin(x))
    ' H4 {) u! N! a3 N& Q8 l) [) C
59. endsub
    3 F2 X' \1 k1 F

复制代码

运行后求得结果如下：
7 `  F: m$ f! f7 a% A( k

$ U1 u* A! c9 ?; I6 `/ w0 _
" u1 I  ?" c( M' n$ ^其中C2列是数值解，C3列是解析解，比较之下，这二者之间无明显差异。

0 U( W1 C' T3 S1 }! S

, g* }6 l$ i( U, E