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

liwenhui

. T6 _! m5 `# p: f  SEViews除了能解决计量经济学的估计问题以外，还提供一个编程环境用以解决复杂的问题。经过调试，我在EViews中实现了用龙格库塔方法求解常微分方程的数值解，供大家交流。
5 r- B  u+ T3 I* Q演示中，我使用了如下常微分方程作为测试：

这个方程的解析通解是：

5 G7 l0 @2 X5 Q2 G0 l0 d; P
使用“龙格库塔方法”，编制的EViews程序如下：

1. '用龙格库塔法求解常微分方程dydx=-y*cos(x)+exp(-sin(x)),y(0)=0在区间[0,10*pi]的数值解8 M5 q' b; p5 o$ z# R0 A
   
2. '已知这个微分方程的解析解 y=exp(-sin(x))*x
   + U2 X6 R3 G# @) O' ~% T* p: i
3. 
   4 h2 i# G5 e/ V3 q
4. '生成一个workfile作为基本的数据容器2 w5 I9 z# j8 P* L\" |' X\" Q
   
5. wfcreate (wf=temp) u 1000# L9 z5 y/ q5 o2 V8 Q8 d% }\" s3 B
   
6. ) o, L- w5 W8 ^: _$ h% C2 `
   
7. '定义常量' u4 |& b\" h9 o5 l
   
8. scalar pi=3.14159
   0 k, g- r# c1 d% i
9. scalar a=0        '定义自变量下限
   ) V% O( u9 K1 u* r2 j1 x) g; m3 ~
10. scalar b=10*pi     '定义自变量上限
    2 y3 y* `) j+ V  v4 o2 I9 p3 W$ _
11. scalar  M=500       '定义步数
    . Y5 K5 ~9 i* i1 F; |( ~- k
12. scalar h=(b-a)/M   '计算每步之间的间隔% J' M( e/ z8 h9 n6 R
    
13. 2 a9 Q: X$ P. m! l
    
14. '定义一个矩阵来储存计算数据，其中第一列储存自变量数据，第二列储存因变量数据,第三列储存解析解的值用以作为比较
    , |. m5 ^$ Q' L$ n% N2 X  |
15. matrix(M+1,3) F 1 ], h  h; i0 |& N% z9 d
    
16. # a5 {4 `) X9 K4 c5 T
    
17. '矩阵的第一行储存初值问题的初始条件
    0 u& [\" x+ U6 ?5 d( O5 v7 A
18. F(1,1)=0- S: N# b* L: z3 i
    
19. F(1,2)=0
    + M' ?0 ~0 s- ^9 ]- j+ y% V9 N( k
20. F(1,3)=@exp(-sin(F(1,1)))*F(1,1)
    ) M7 }9 ^7 G. B1 E7 Q1 ?5 U4 [
21. + p/ {: i0 Q+ d\" o
    
22. '定义龙格库塔法的权重参数
    8 b. L+ A5 g, |3 M( X* [
23. scalar k1
    ; \0 u4 P: `7 B' ^4 [, t7 H
24. scalar k2
    # s! A' Z, W+ i3 V
25. scalar k3
    8 L; U. `$ s9 @, x2 k7 b
26. scalar k4! D2 |: K- c* N
    
27.   P+ z3 F( E, n  {: N! D( v, {8 ?
    
28. '定义权重的过程量( X/ M/ g8 u\" H
    
29. scalar w1$ R  }% `- A# w# P\" o# B3 q\" X; j
    
30. scalar w27 p& G& e6 u- v5 r1 W# m5 q
    
31. scalar w3: b1 I$ P- j9 \
    
32. scalar w4
    0 \0 N: `  I, q! g6 G
33. 
    + w  L/ x6 i3 Y/ F\" V( |7 z
34. '程序主体9 x3 K( O  z) c; `! j7 ]( m
    
35. for !k=1 to M step 1
    5 W$ O8 l( ^1 Y2 k( Q! o
36.   F(!k+1,1)=F(!k,1)+h
    * R6 s! y+ }0 ^/ ?  _
37.   '调用常微分方程计算权重( Y) f0 E1 Q) r* x1 Y) ^# B0 _! V
    
38.   call obj(w1,F(!k,1),F(!k,2))
    9 j6 @8 z- C, q8 \5 D
39.     k1=w1*h/ e! W, c5 P9 }\" F' B0 ^; A
    
40.   call obj(w2,F(!k,1)+h/2,F(!k,2)+k1/2): w  Y2 Z. D  B8 Q+ e
    
41.     k2=w2*h\" Y4 b9 M# [1 q3 Y$ I
    
42.   call obj(w3,F(!k,1)+h/2,F(!k,2)+k2/2)
    * X+ Z5 L: {  M  w; o  q
43.     k3=w3*h
    ' \1 ^+ k# P2 F
44.   call obj(w4,F(!k,1)+h,F(!k,2)+k3)$ K) {, M4 s! A5 q
    
45.     k4=w4*h6 q! ?8 t. M, F7 {7 G) A8 Y
    
46.   '计算函数估计值! V7 _: Q0 j% I: h2 W* p- b
    
47.   F(!k+1,2)=F(!k,2)+(k1+2*k2+2*k3+k4)/69 z! S0 O; u/ g  f5 e9 v, j
    
48.   '计算函数解析值
    6 [7 Y( p0 u1 C\" C+ ~
49.   F(!k+1,3)=@exp(-sin(F(!k+1,1)))*F(!k+1,1)
    1 G+ B\" M: G; V' y, ]% ?3 q7 r% v
50. next
    3 C# g. j: R: n! L: P
51. ; S5 A( a3 A) e\" s/ [: h9 w
    
52. '显示最终结果9 W3 i+ A6 l8 k\" a- Q) W( s' y0 l+ g
    
53. freeze F.xyline
    $ ]1 g$ N( K8 {  ?; J: D& X! `
54. freeze F
    & d5 C: F/ D# i% |  l0 g$ _
55. $ B1 m0 N# z1 A! y# a9 H- L, ]8 L
    
56. '定义常微分方程
    ' e! p1 U, o; o$ l) B0 X
57. subroutine obj(scalar dydx,scalar x,scalar y)7 L+ X5 V' O, c8 Q8 n% m; Y
    
58.   dydx=-y*@cos(x)+@exp(-@sin(x))
    2 ?0 F, m) C! U) j. W
59. endsub! ~' Y3 Y( ]. h( ]3 p
    

复制代码

运行后求得结果如下：

( b0 U7 R- p6 W

其中C2列是数值解，C3列是解析解，比较之下，这二者之间无明显差异。

: m. Y/ l: S. q" o+ h