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

liwenhui

- v, L5 E8 ]( F2 VEViews除了能解决计量经济学的估计问题以外，还提供一个编程环境用以解决复杂的问题。经过调试，我在EViews中实现了用龙格库塔方法求解常微分方程的数值解，供大家交流。
演示中，我使用了如下常微分方程作为测试：

这个方程的解析通解是：

4 \2 N! K: P  B3 H8 G2 D
3 @' O9 C% `* c9 T- U使用“龙格库塔方法”，编制的EViews程序如下：

1. '用龙格库塔法求解常微分方程dydx=-y*cos(x)+exp(-sin(x)),y(0)=0在区间[0,10*pi]的数值解
   0 F, s. ]& M2 O& s& T* ~; S5 `
2. '已知这个微分方程的解析解 y=exp(-sin(x))*x; |7 S* X\" ~* I, p
   
3. , \5 o! d8 @9 w7 K
   
4. '生成一个workfile作为基本的数据容器
   3 e+ B0 {, n9 K; }0 q2 |1 Q
5. wfcreate (wf=temp) u 1000
   \" n* q# |2 s. W# X, J
6. 2 I+ X. I. _' i0 S
   
7. '定义常量3 P! @; @! R1 \# r
   
8. scalar pi=3.14159
   ) k\" x1 G9 r5 P
9. scalar a=0        '定义自变量下限2 ~8 N& I) W% ?1 `
   
10. scalar b=10*pi     '定义自变量上限
    + t\" g+ u$ {+ S! B3 Y5 S1 `7 e
11. scalar  M=500       '定义步数2 Q  Q# O* Q# t( D1 k
    
12. scalar h=(b-a)/M   '计算每步之间的间隔
    ( x+ o, v, u, `
13. 
    + C$ A7 h: {% s\" w& K' ~; K2 s
14. '定义一个矩阵来储存计算数据，其中第一列储存自变量数据，第二列储存因变量数据,第三列储存解析解的值用以作为比较
    % J1 e0 I\" `# L2 v
15. matrix(M+1,3) F 2 G2 U9 @7 f) d2 v6 D
    
16. 0 l% X! Z; n! I+ n
    
17. '矩阵的第一行储存初值问题的初始条件
    : ^: `% J0 E1 t\" ^- h0 X2 }0 N
18. F(1,1)=0
    7 N9 g1 v# D% w9 j' H+ c
19. F(1,2)=0
    5 H3 s! Y$ _0 c6 F
20. F(1,3)=@exp(-sin(F(1,1)))*F(1,1)
    6 U+ o/ j- n* @7 Y: p
21. 
    1 b3 q/ i2 z$ B  _' N: r0 s
22. '定义龙格库塔法的权重参数
    - F! ^- X$ `3 t3 h( T
23. scalar k1
    # F3 V$ U\" w/ y  r# L7 Q
24. scalar k2
    : l9 g3 Y7 l8 W% U: Q0 q# B& b/ s
25. scalar k31 [; s/ G' Z: h6 e
    
26. scalar k4
      w\" I/ }8 C& W( o8 p
27. ! H2 [0 X  \5 l# G( v( n
    
28. '定义权重的过程量6 e! X5 \6 ~% v- \
    
29. scalar w1
    8 O, s+ {+ r  E/ j' f( n1 d
30. scalar w2: V0 P: }3 E. D\" z5 i+ x% V
    
31. scalar w3
    - W$ C: ^1 V6 m/ V% s
32. scalar w4
    9 {  L! x9 a6 o! v( K8 z
33. 
    \" f3 e) v. z\" H/ Z% C2 P3 ?/ m- q
34. '程序主体3 R6 L; n- K/ @1 r
    
35. for !k=1 to M step 1
    + d! ]: H7 [# m) ^
36.   F(!k+1,1)=F(!k,1)+h
    9 o3 y' s( P, Q! C5 l0 l  X
37.   '调用常微分方程计算权重( }1 c+ [2 r+ m' Z: o7 X6 j8 j4 ?
    
38.   call obj(w1,F(!k,1),F(!k,2))
    2 O7 W  T) z! b3 k
39.     k1=w1*h, ^5 A3 m  o: E# u) D- t) v; T
    
40.   call obj(w2,F(!k,1)+h/2,F(!k,2)+k1/2)
    : j9 j: A4 r& Y* Y. n9 X1 q
41.     k2=w2*h3 u! x2 \- W% K
    
42.   call obj(w3,F(!k,1)+h/2,F(!k,2)+k2/2)) @, F8 f; ?: P. _* r7 L7 A* u
    
43.     k3=w3*h
    ; A/ F: c6 j3 I' U# ?' l8 Z5 q1 n
44.   call obj(w4,F(!k,1)+h,F(!k,2)+k3)
    % C! X6 O% h\" k1 |5 g6 y\" n% K
45.     k4=w4*h5 X2 r, X# @! X: `& ^. s
    
46.   '计算函数估计值
    . G# n9 h7 n  e9 ^9 V: s) Z3 j
47.   F(!k+1,2)=F(!k,2)+(k1+2*k2+2*k3+k4)/6
    ) j- j& l+ G9 O, @6 n: U, S( Y- ]2 z
48.   '计算函数解析值9 t& |$ W: I\" a
    
49.   F(!k+1,3)=@exp(-sin(F(!k+1,1)))*F(!k+1,1)
    0 ~  j7 X% Q) Z9 {) i% I; R
50. next\" o; |* H# C' O% f
    
51. 1 [9 k  V; x\" ]+ y5 D0 s% \# c
    
52. '显示最终结果/ D0 |9 h- U, S/ ~1 p  V! _( c; G' h
    
53. freeze F.xyline' _2 h6 ]& W' w4 u9 R
    
54. freeze F
    3 |  p. M3 n+ f0 Z' g
55. 
    2 Z\" E\" H4 |+ g& k) z1 d% ?; m
56. '定义常微分方程
    ' o1 U$ d+ f; u& P+ N& P
57. subroutine obj(scalar dydx,scalar x,scalar y)
    + a* v# J  T\" D  P0 {( M* Q6 `0 y
58.   dydx=-y*@cos(x)+@exp(-@sin(x))
    : Y: n! O& n2 N% T8 ^
59. endsub
    8 h# Z\" Y! I9 N8 w3 Z5 [

复制代码

运行后求得结果如下：
. [; Z( |1 b8 G/ o

3 w* b. k# R+ c- U2 Z, \
; l. q! [7 u- i& a. s( X其中C2列是数值解，C3列是解析解，比较之下，这二者之间无明显差异。
" s( P. f7 D. i
5 X6 Q$ M$ e% [, g1 I9 N

  p, |+ o* ~8 ]5 U. P7 O
: [1 T, _: I) G, K& s# ^