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

liwenhui

( U' h: U# [% Q9 Q
- ]+ {; o) g9 h% VEViews除了能解决计量经济学的估计问题以外，还提供一个编程环境用以解决复杂的问题。经过调试，我在EViews中实现了用龙格库塔方法求解常微分方程的数值解，供大家交流。
演示中，我使用了如下常微分方程作为测试：

+ }% w3 j6 B/ Z, s0 B1 V- A这个方程的解析通解是：
* c) p" N# f9 g

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

1. '用龙格库塔法求解常微分方程dydx=-y*cos(x)+exp(-sin(x)),y(0)=0在区间[0,10*pi]的数值解
   ( F6 ~' R0 l( e9 U+ h  Y& _
2. '已知这个微分方程的解析解 y=exp(-sin(x))*x
   / @4 l) z\" e+ O\" X\" E
3. 
   ; Y) p5 S2 o, I
4. '生成一个workfile作为基本的数据容器
   6 o* q% u# ~\" ~: t6 _# Q
5. wfcreate (wf=temp) u 1000
   \" e\" |7 D4 b) ]6 F  r7 J( [7 ^2 B
6. . e& X; R! w0 d2 ]
   
7. '定义常量
   , U2 `: C( ]7 u1 y\" K) j8 j
8. scalar pi=3.14159- O+ y+ V+ T6 F
   
9. scalar a=0        '定义自变量下限* o- ~: V3 G\" b' K  [: Z
   
10. scalar b=10*pi     '定义自变量上限7 R6 X+ l6 X: b# ~0 }
    
11. scalar  M=500       '定义步数
    ' o1 G/ ?) G2 @& X3 B
12. scalar h=(b-a)/M   '计算每步之间的间隔/ z5 P- m! S! Y) e' J5 H5 V4 d\" d
    
13. 
    9 v9 H  U% D* \/ q
14. '定义一个矩阵来储存计算数据，其中第一列储存自变量数据，第二列储存因变量数据,第三列储存解析解的值用以作为比较
    ' A# _0 u* k+ h\" }* N7 {4 F
15. matrix(M+1,3) F
    1 m9 d) ]* p, Y( V' k, j. y
16. 
    2 w0 [, D/ \7 d+ M, C' u' r! h\" u
17. '矩阵的第一行储存初值问题的初始条件
    ; J# g\" f' `5 \) ]2 ]\" @
18. F(1,1)=0, P4 }8 l0 w8 Q- R
    
19. F(1,2)=0
    5 D+ C* A2 {! h: j) c! Q  l# D
20. F(1,3)=@exp(-sin(F(1,1)))*F(1,1), H1 S2 D. R) T/ g
    
21. 
    8 N1 H& P! H+ g( a  C
22. '定义龙格库塔法的权重参数
    3 H7 |1 g7 L\" e  h6 ^4 l
23. scalar k17 O) L: @  E' L& U- q* N/ ]& y
    
24. scalar k2* b- x/ i/ u* g+ T\" B. ?& Z
    
25. scalar k3
    & r4 ^+ z. O  F& Y5 o0 P
26. scalar k4
    4 `/ z8 J& q# ]& B: u
27. . }/ J$ s6 X/ V9 j
    
28. '定义权重的过程量\" u* A; F: U& g9 V5 f
    
29. scalar w1
    / R0 e7 O- |3 v. Y) k# @1 m4 O1 h2 y
30. scalar w2; _4 z) n; m4 ]: a
    
31. scalar w3. O7 ~, i9 X# W. h
    
32. scalar w4
    . J7 J; c+ L! M  J0 D
33. 
    4 b  c7 D% h) h
34. '程序主体
      M1 t4 h: b: w6 ~) K
35. for !k=1 to M step 10 E: S  H  W) h* |) e. q3 `
    
36.   F(!k+1,1)=F(!k,1)+h' S* p/ m0 E7 y& z
    
37.   '调用常微分方程计算权重
    & B7 Z3 X1 Y1 ~9 d( w4 {
38.   call obj(w1,F(!k,1),F(!k,2))
    & p7 n2 r9 B, s8 A# L/ n' M) V5 g
39.     k1=w1*h6 C; {3 G7 }\" L8 Z
    
40.   call obj(w2,F(!k,1)+h/2,F(!k,2)+k1/2)) a* w3 l$ u- D9 _, W2 I
    
41.     k2=w2*h# U5 V$ _, s  H' |2 B3 O
    
42.   call obj(w3,F(!k,1)+h/2,F(!k,2)+k2/2)
    ; `2 o- ~9 n$ H) S' e; K% |2 Z  k+ l
43.     k3=w3*h
    0 R4 P7 L: o  \/ n; J
44.   call obj(w4,F(!k,1)+h,F(!k,2)+k3)- D1 [% G, U. @  t5 o- ]1 J$ O
    
45.     k4=w4*h
    ) Z0 J1 F$ m+ p\" H: X
46.   '计算函数估计值# e\" }& m$ o; V# N) i
    
47.   F(!k+1,2)=F(!k,2)+(k1+2*k2+2*k3+k4)/6
    0 A! w& @, Y1 R+ G
48.   '计算函数解析值9 V2 I7 f0 ?+ G
    
49.   F(!k+1,3)=@exp(-sin(F(!k+1,1)))*F(!k+1,1)+ `; |# Q+ t! Z  B\" B# O* \: ]) j
    
50. next0 K/ b& f' p% S! L
    
51. 
    $ S8 J, m$ p' C6 A+ `+ j
52. '显示最终结果* R' L7 M; u: {\" x; z0 U8 p
    
53. freeze F.xyline, h$ S  F1 Z6 S$ E. f; b
    
54. freeze F
    1 j7 t( J$ {0 D4 T' |* f
55. \" m& Q# w4 [7 Z\" }/ p
    
56. '定义常微分方程
    / i0 A. @1 _' |8 s( N
57. subroutine obj(scalar dydx,scalar x,scalar y)
    : O. U3 A+ H: g+ M/ b- h
58.   dydx=-y*@cos(x)+@exp(-@sin(x))  o& _& F\" A0 `\" b2 A
    
59. endsub
    ! j( u6 l  f1 K4 m4 X

复制代码

运行后求得结果如下：

; m# z" b- I- O" F

其中C2列是数值解，C3列是解析解，比较之下，这二者之间无明显差异。
0 j) H. [1 ?  W
  ?; |) z" u9 Y. Z% ^7 p# W

/ Z$ W9 U  W! r6 L/ n! t: T
" B# o+ T* a& o' @, m+ i0 F/ @& G6 O