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

liwenhui

$ j; |$ X. R' p: hEViews除了能解决计量经济学的估计问题以外，还提供一个编程环境用以解决复杂的问题。经过调试，我在EViews中实现了用龙格库塔方法求解常微分方程的数值解，供大家交流。
演示中，我使用了如下常微分方程作为测试：
( c! R3 l7 p0 z! N. m5 ^

7 |; {* V0 O; H) d8 g. |: I
这个方程的解析通解是：
5 d0 Z/ F6 l8 y+ J

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

1. '用龙格库塔法求解常微分方程dydx=-y*cos(x)+exp(-sin(x)),y(0)=0在区间[0,10*pi]的数值解
   ' i( ~7 E- o. P& X
2. '已知这个微分方程的解析解 y=exp(-sin(x))*x
   + l/ `, a/ w% g3 q8 B# G& f- ^
3. 
   & m& P. T9 Y# m7 R9 _5 @
4. '生成一个workfile作为基本的数据容器
   , V% c* l: f* d* Z+ a, `
5. wfcreate (wf=temp) u 1000' d/ ]( l; l2 Y% D( |! d. ^
   
6. 2 e$ i\" C8 R) z
   
7. '定义常量; W9 F\" v: |, [% z7 y: j
   
8. scalar pi=3.14159! h7 o- ?+ h6 n% l) x
   
9. scalar a=0        '定义自变量下限
     O. M5 y6 g0 K
10. scalar b=10*pi     '定义自变量上限8 d' m2 K0 ^* Z& t
    
11. scalar  M=500       '定义步数
    $ I2 J8 h% H) N, N) u& E
12. scalar h=(b-a)/M   '计算每步之间的间隔
    / b\" {8 A/ s) f1 j7 ^7 ~! M
13. 
    ' d& u4 S& n$ W' q
14. '定义一个矩阵来储存计算数据，其中第一列储存自变量数据，第二列储存因变量数据,第三列储存解析解的值用以作为比较
    . G# K* V0 F6 z# l3 U$ R& B2 h
15. matrix(M+1,3) F
    \" M# X\" D  c8 x' t! A3 H
16. 
    + W7 F( c/ y! A/ p% D
17. '矩阵的第一行储存初值问题的初始条件
    % `; p) x+ q! U
18. F(1,1)=0
    4 }5 w0 B) X. z! u6 C6 D0 ^
19. F(1,2)=03 Q  B4 J\" t( E) x
    
20. F(1,3)=@exp(-sin(F(1,1)))*F(1,1)1 G+ f% W# }$ J5 Z# Z. j
    
21. 
    9 }- f\" Q7 b! t, v! R1 N( l
22. '定义龙格库塔法的权重参数
    7 e, R7 x* W3 ^1 \  y* m6 m7 Z
23. scalar k1\" v% _6 z' s, k, R) C5 P
    
24. scalar k2
    5 l$ Y! \, |* S$ U1 z1 R0 L: l/ Y7 e
25. scalar k3
    + V$ M7 r& Q# \2 g2 X9 Q: t
26. scalar k45 u) B) \$ j' C# G\" E& Y9 t
    
27. 
    # F! R9 c! ~0 Q; R/ A
28. '定义权重的过程量
    / a, Z7 D* d3 o0 _. ?8 t* A3 q' N
29. scalar w1( Z\" D; p. y( K3 Z0 L1 K
    
30. scalar w20 F- W9 b\" q6 F
    
31. scalar w3
    ; u6 M) \9 N, H. u4 H) f# s+ w
32. scalar w4
    # P# t3 s6 \$ O  w4 B
33. 
    0 f( f/ t( O8 `$ Y# Z7 E
34. '程序主体+ D: A( e! X; {4 k
    
35. for !k=1 to M step 1
    - Y: }- u! ^2 \. g3 m' j
36.   F(!k+1,1)=F(!k,1)+h
    3 r' P! \+ y. _; m8 L
37.   '调用常微分方程计算权重
    ! V7 o; X( w  k
38.   call obj(w1,F(!k,1),F(!k,2)) - b1 f8 s- E% K- ^& ^; c& c9 @
    
39.     k1=w1*h% J; T: K6 x6 x* b+ r! g) h
    
40.   call obj(w2,F(!k,1)+h/2,F(!k,2)+k1/2)+ d7 x5 t2 o5 Y- t  Q$ X- ^
    
41.     k2=w2*h  Q  Q2 W/ {. I. M) R- T$ G7 \\" Q
    
42.   call obj(w3,F(!k,1)+h/2,F(!k,2)+k2/2)
    4 A2 z* T1 j( h6 f+ D4 ^' G
43.     k3=w3*h+ }2 p3 A+ z# I2 a7 e, x+ J, l1 l
    
44.   call obj(w4,F(!k,1)+h,F(!k,2)+k3)
    5 N( R9 m- d\" ?% m' d8 L0 s  i
45.     k4=w4*h& h; d# ?1 ?: ~( A
    
46.   '计算函数估计值
    7 b' v$ V1 u5 g3 K; o4 u2 S/ Z
47.   F(!k+1,2)=F(!k,2)+(k1+2*k2+2*k3+k4)/6
    6 E- v, D( i5 p+ g) @2 o
48.   '计算函数解析值
    \" r* v6 V5 o8 j5 h7 r/ n
49.   F(!k+1,3)=@exp(-sin(F(!k+1,1)))*F(!k+1,1)
    5 ]; S6 H8 p4 j2 l
50. next
    ' G( z# J5 I: v8 g# K
51. 
    ) @3 }! d7 n\" r0 _( [\" ]' \
52. '显示最终结果
    5 g: u# F, A, M' k9 Z1 f% h
53. freeze F.xyline( Q5 j/ O6 ?6 K9 @/ m' J( k
    
54. freeze F
    ' Z1 ^\" {. _& h9 C1 g# m0 s
55. 
    \" M2 J  I- W5 P\" H3 b# l+ s
56. '定义常微分方程
    \" e- d( O\" A' U9 B) K! o1 B7 {
57. subroutine obj(scalar dydx,scalar x,scalar y). T- D5 p1 Q' Q: ]# n( X
    
58.   dydx=-y*@cos(x)+@exp(-@sin(x))
    $ q) N/ l! E& b& A4 i$ y
59. endsub
    3 n, F! P3 ~3 G1 }$ B0 @

复制代码

运行后求得结果如下：

3 [6 Z5 t* A9 w; T: I" R  m

7 q# b3 o6 Z2 F4 r6 `- p; O3 K其中C2列是数值解，C3列是解析解，比较之下，这二者之间无明显差异。
( H# D( j  b/ C5 T/ f
& }' U6 A$ D' |1 W7 [