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

liwenhui

5 _% K( w, k8 s8 ~
1 Z9 E8 H: l/ r+ y7 P5 t% v7 b2 EEViews除了能解决计量经济学的估计问题以外，还提供一个编程环境用以解决复杂的问题。经过调试，我在EViews中实现了用龙格库塔方法求解常微分方程的数值解，供大家交流。
- O. r+ Q9 [2 ~( e演示中，我使用了如下常微分方程作为测试：

5 g8 {( r% R* q- G( B
这个方程的解析通解是：
) ~5 o1 W* F  I3 _$ s) H! a

# Q+ m, G7 N4 t" Z) Q使用“龙格库塔方法”，编制的EViews程序如下：

1. '用龙格库塔法求解常微分方程dydx=-y*cos(x)+exp(-sin(x)),y(0)=0在区间[0,10*pi]的数值解! l5 W7 M& r+ Q1 J  m, d% D' X1 j
   
2. '已知这个微分方程的解析解 y=exp(-sin(x))*x
   \" a/ `/ o) ~8 }% w
3. 
   9 S) v$ D4 W$ H6 J( u
4. '生成一个workfile作为基本的数据容器
   & l; b1 g5 S* Q# }# f
5. wfcreate (wf=temp) u 1000
   ) Q9 o: k3 g' V, J  r9 ^/ D
6. + G6 I2 }8 S  J8 n5 I- V
   
7. '定义常量$ D: p; X' o# a  [0 N6 v: i
   
8. scalar pi=3.14159
   ! Y+ R% |( `, @1 \. }' D+ r1 R
9. scalar a=0        '定义自变量下限
   % _# J\" o9 t: z$ @+ a/ I4 C
10. scalar b=10*pi     '定义自变量上限: ~. x% D% [3 F
    
11. scalar  M=500       '定义步数
    7 [0 H; Q6 e3 M3 ^* h$ C
12. scalar h=(b-a)/M   '计算每步之间的间隔
    0 n4 R4 e  E: R
13. 
    * Y* H! s* C- j; J# J
14. '定义一个矩阵来储存计算数据，其中第一列储存自变量数据，第二列储存因变量数据,第三列储存解析解的值用以作为比较
    * `' S  @9 R1 b3 |# A  G+ D
15. matrix(M+1,3) F
    # S6 q# M4 Q0 D
16. 
    5 v% s4 P( \, T. D$ J
17. '矩阵的第一行储存初值问题的初始条件
    # f7 _+ |7 e, k. O( W1 ]3 O
18. F(1,1)=0% {) S. n6 a4 @% j8 Y/ q
    
19. F(1,2)=0
    7 w\" \5 R3 F1 U9 V+ u, z' A; F
20. F(1,3)=@exp(-sin(F(1,1)))*F(1,1)
    2 w' N4 Q: k8 l7 g6 p
21. 8 L6 R0 W; i; g8 H* s
    
22. '定义龙格库塔法的权重参数
    7 ?# Z5 u9 `2 y\" ?8 V/ L
23. scalar k1+ t  g6 ^( f0 q$ _1 Q( ~% u# g
    
24. scalar k2* u8 b. P, D$ [% p$ n* }! D4 D
    
25. scalar k3' R2 h4 w) |: S- t9 h9 [
    
26. scalar k4, q2 a0 U7 R$ J2 M, X2 }7 l# f
    
27. 
    1 M\" P1 Z8 Y8 D' e+ k+ v
28. '定义权重的过程量3 w& @* B) O: Y& v8 J! t; S
    
29. scalar w16 d: X1 [/ `3 Z( Y) i\" t
    
30. scalar w2
    - ~$ Y( }: O* w1 \
31. scalar w3\" s7 F8 I# N, r2 e$ S4 v
    
32. scalar w4; V+ H& t: H, H& r
    
33. 
      \# q! v. I) N  A7 ^
34. '程序主体
    \" C; |0 t& b8 W8 }& R9 [& l* P
35. for !k=1 to M step 1# _\" q\" R* A+ ?- Z8 e\" ?
    
36.   F(!k+1,1)=F(!k,1)+h4 R% s% p5 [- X. e9 N3 _
    
37.   '调用常微分方程计算权重
    ( l( B: Q0 z3 x
38.   call obj(w1,F(!k,1),F(!k,2))
    3 r; V# C, n0 _1 P0 \
39.     k1=w1*h
    ' W2 t# C7 n* |/ r
40.   call obj(w2,F(!k,1)+h/2,F(!k,2)+k1/2)8 `2 D- d5 Z# q
    
41.     k2=w2*h% ]1 Y. x7 o& ^# L
    
42.   call obj(w3,F(!k,1)+h/2,F(!k,2)+k2/2)( P; d+ `8 p5 q/ X
    
43.     k3=w3*h
    + X1 U\" O+ N5 P/ u( l\" c* Y4 w
44.   call obj(w4,F(!k,1)+h,F(!k,2)+k3)# S6 R3 w2 N* u+ b; y* I/ e- J
    
45.     k4=w4*h
    * h5 x3 T\" [* }
46.   '计算函数估计值
    , I- P5 {( Q6 C4 H4 c* n
47.   F(!k+1,2)=F(!k,2)+(k1+2*k2+2*k3+k4)/68 ~8 ~) z; G3 ?8 }/ R
    
48.   '计算函数解析值
    ! B) y: H5 X1 r\" G1 q  `9 n# B
49.   F(!k+1,3)=@exp(-sin(F(!k+1,1)))*F(!k+1,1)! N$ U, W# V  R! G3 c
    
50. next\" {: X6 w# h6 ~6 y
    
51. % ]- I; I2 ?. [: @/ u
    
52. '显示最终结果8 ^# e5 V$ i; i9 e0 z
    
53. freeze F.xyline
    1 ]* E1 \- W' Q$ ?5 m6 c
54. freeze F7 w6 ?9 E1 l4 F1 j8 l
    
55. / j5 t\" n% ]' j1 Q
    
56. '定义常微分方程
    . P* N* w+ `. J/ i5 L/ l
57. subroutine obj(scalar dydx,scalar x,scalar y)
    ' B, k( j. X4 _& y: I- P5 L8 y
58.   dydx=-y*@cos(x)+@exp(-@sin(x))
    - t6 W- n. F3 Z2 Q, k2 \# U9 O
59. endsub- n9 q! N* o# g  ~* S6 A: [
    

复制代码

运行后求得结果如下：
4 t2 ]9 L* r5 s- r; I
/ w. r% i$ U; z4 Q+ P; O

4 p4 X8 a0 r! C其中C2列是数值解，C3列是解析解，比较之下，这二者之间无明显差异。
' K' m$ v) g0 r6 B8 A& a
3 H6 e1 i, S2 F3 X0 a" `

4 X' Y' e! E2 n( V
% v2 l( {9 p9 ^2 ], D3 ]