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

liwenhui

( Y. U4 \: C/ y1 C6 J1 G
" l6 L' q4 P  B9 P2 @EViews除了能解决计量经济学的估计问题以外，还提供一个编程环境用以解决复杂的问题。经过调试，我在EViews中实现了用龙格库塔方法求解常微分方程的数值解，供大家交流。
演示中，我使用了如下常微分方程作为测试：
2 s8 D$ S. v5 @' d  s1 K/ _& f, s; l% L" V

& U1 M$ h7 n) U' K4 j; `4 w- `
这个方程的解析通解是：

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

1. '用龙格库塔法求解常微分方程dydx=-y*cos(x)+exp(-sin(x)),y(0)=0在区间[0,10*pi]的数值解: o2 [  g1 N  @7 c\" B7 ]9 t
   
2. '已知这个微分方程的解析解 y=exp(-sin(x))*x
   ( J1 d6 v! O: w- d. V8 o7 Z, \: a* |
3. 
   - ?; c2 ~1 A' w) m* f4 T# U
4. '生成一个workfile作为基本的数据容器8 T6 Q& `5 B7 h4 Z
   
5. wfcreate (wf=temp) u 1000- P# e( i3 a% W/ }- M
   
6. 
     S) A* d. B) P7 f
7. '定义常量
   8 p. E8 b- P' K- X# R
8. scalar pi=3.14159
   ; F6 D9 \1 q4 F1 r
9. scalar a=0        '定义自变量下限( c; U$ L- i. g% E/ t% n9 {
   
10. scalar b=10*pi     '定义自变量上限) U& Z7 y( s  w0 T
    
11. scalar  M=500       '定义步数; z7 F- J$ s% {0 Y: |
    
12. scalar h=(b-a)/M   '计算每步之间的间隔
    ) j8 e# J% `- [5 a+ @, i5 Z' o1 F
13. - h. V6 `: [# T, C* u8 t+ i
    
14. '定义一个矩阵来储存计算数据，其中第一列储存自变量数据，第二列储存因变量数据,第三列储存解析解的值用以作为比较% x  S2 J* |5 \/ j/ ~
    
15. matrix(M+1,3) F . `6 t4 [$ T# \9 W- N
    
16. 
    . E7 C4 F7 S: ?
17. '矩阵的第一行储存初值问题的初始条件2 F8 Q- C% `& h/ K  ?5 k. V( g+ O
    
18. F(1,1)=03 H\" O& Z. M& y  j0 ]
    
19. F(1,2)=0
    4 x! C) u1 k2 p# K1 j. s* W
20. F(1,3)=@exp(-sin(F(1,1)))*F(1,1)  g4 E# x/ B4 }: t. l
    
21. ( x5 b\" S+ v9 e6 O
    
22. '定义龙格库塔法的权重参数% ]$ a% L4 o/ D& h
    
23. scalar k1' [& D6 n# m8 w\" `$ `
    
24. scalar k2  W' K5 Y8 L4 v  T2 [; s. K3 ~
    
25. scalar k37 Q1 v$ S4 Y( Z0 A
    
26. scalar k4
    * F9 V! O+ p3 A% n; {
27. 
    3 J9 x3 h\" {# K) |
28. '定义权重的过程量
    0 R+ H2 n# E( v1 \2 x
29. scalar w1
    4 N0 K6 K2 @6 K) L
30. scalar w2
    ' i9 {\" _/ h6 B- J' ], }
31. scalar w30 V, Y# I4 M1 a/ M
    
32. scalar w4
    ! f\" _% n: {5 w/ ?
33. 
    2 ^$ _9 Y1 F\" S$ M/ K
34. '程序主体* j  d. ^& ^3 I  Y
    
35. for !k=1 to M step 1
    4 [# z) b- u. D\" }& E6 e
36.   F(!k+1,1)=F(!k,1)+h; D0 d; j* K- B' M: d\" {5 }
    
37.   '调用常微分方程计算权重
    1 q. T0 P7 \! _
38.   call obj(w1,F(!k,1),F(!k,2)) \" |\" c0 n! W$ P9 d( N- Q
    
39.     k1=w1*h/ ]% ~4 @, E7 g& V0 F\" z0 Q) \
    
40.   call obj(w2,F(!k,1)+h/2,F(!k,2)+k1/2)4 \9 S- X( x; N, k' \/ \' Y
    
41.     k2=w2*h
      q0 V5 w( Y+ B, D2 i
42.   call obj(w3,F(!k,1)+h/2,F(!k,2)+k2/2)& q! l) r) t, f, G0 {9 d
    
43.     k3=w3*h4 j( P, p8 {8 C1 [\" U& d
    
44.   call obj(w4,F(!k,1)+h,F(!k,2)+k3)
    ; R) G, _2 b+ a* _/ H& I. A
45.     k4=w4*h$ a6 _8 s( O! d; |& Z\" S6 v
    
46.   '计算函数估计值
    * n0 a. D# U\" d2 A& J3 {
47.   F(!k+1,2)=F(!k,2)+(k1+2*k2+2*k3+k4)/6
    + f  {' e; a8 [9 k+ v1 a& t
48.   '计算函数解析值& c2 @# K3 ]  o, Z( ^
    
49.   F(!k+1,3)=@exp(-sin(F(!k+1,1)))*F(!k+1,1)2 ]# U& v\" |5 Y
    
50. next\" K% r0 x. l9 Q\" E
    
51. + R5 a3 K% {4 L  Y9 x
    
52. '显示最终结果
    8 k# k+ o4 J# [0 s* a0 o% V! m3 S
53. freeze F.xyline
    ! O6 O/ U- O7 m* R! u1 ?: G. z
54. freeze F! o8 g, j7 @6 I( i! u$ S# i( o! {8 i
    
55. 
    ! s, K  d/ ^6 W- F( i, d
56. '定义常微分方程
    5 y* A3 n( T2 ]) u\" p+ e0 T
57. subroutine obj(scalar dydx,scalar x,scalar y)
    $ T) h\" f4 e8 l+ }, {0 a: w  P. V  L
58.   dydx=-y*@cos(x)+@exp(-@sin(x))
    0 N: Y: z( O5 C
59. endsub% D1 Z\" X% s2 F3 M% F. r0 H' k# g
    

复制代码

运行后求得结果如下：

1 |' C0 i0 g; p- I/ e6 e2 d
- c- }  ?2 q2 j# A其中C2列是数值解，C3列是解析解，比较之下，这二者之间无明显差异。


  r" k+ g5 ^' e; C

7 K1 ~7 F: a) U% A) e- L