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

liwenhui

3 w4 ]# @8 P4 a( a0 `" T% eEViews除了能解决计量经济学的估计问题以外，还提供一个编程环境用以解决复杂的问题。经过调试，我在EViews中实现了用龙格库塔方法求解常微分方程的数值解，供大家交流。
演示中，我使用了如下常微分方程作为测试：
% a' n0 F- b4 s/ D2 \) @0 d5 e

$ p. r/ B, T  }4 s7 U) ^
这个方程的解析通解是：
# L: z6 ^+ E+ G2 i7 t3 U

7 K4 M; H6 E# c; u/ f
使用“龙格库塔方法”，编制的EViews程序如下：

1. '用龙格库塔法求解常微分方程dydx=-y*cos(x)+exp(-sin(x)),y(0)=0在区间[0,10*pi]的数值解; r: n: J( E2 n/ c, B
   
2. '已知这个微分方程的解析解 y=exp(-sin(x))*x4 J( s9 T3 j) A. q
   
3. 3 P, t, f0 L# D  u# Q- @# G* }# C& H
   
4. '生成一个workfile作为基本的数据容器
   3 t) m8 o8 N6 \! P9 w/ P
5. wfcreate (wf=temp) u 1000+ c  X7 o+ n1 _4 I5 H3 ~* E  z
   
6. 
   ( K( q' g5 ]# k5 X/ v, q
7. '定义常量
   5 l' u: {7 f3 @: f$ c
8. scalar pi=3.14159
   ! ^3 W0 B8 G4 b, w: J) {1 T4 f; A
9. scalar a=0        '定义自变量下限- O( f+ W2 [# M+ i$ C  A3 n
   
10. scalar b=10*pi     '定义自变量上限# L# Y8 C  ~4 t$ E\" V
    
11. scalar  M=500       '定义步数/ a: O8 u: D4 F0 B; o
    
12. scalar h=(b-a)/M   '计算每步之间的间隔
    : I\" h; b. M+ h: R
13. 
    ( z4 |5 n$ i$ e& _4 N
14. '定义一个矩阵来储存计算数据，其中第一列储存自变量数据，第二列储存因变量数据,第三列储存解析解的值用以作为比较2 V, L/ I& \& W
    
15. matrix(M+1,3) F
    4 D& X/ t( W, i: I4 ?
16. 
    6 I! k( I. @% _$ A! ?0 s
17. '矩阵的第一行储存初值问题的初始条件# S; ^) ]- D' c  `( [! n
    
18. F(1,1)=0
    ' R5 Q* ?& H: q. j7 @3 ?+ e
19. F(1,2)=0
    ! ~7 o/ f: l2 L# k, ~1 f
20. F(1,3)=@exp(-sin(F(1,1)))*F(1,1)
    ' }1 O  Z+ D# f
21. 
    6 c8 [) _1 |4 i8 [
22. '定义龙格库塔法的权重参数+ e) ?9 W4 _8 n& q# J: l
    
23. scalar k1
    5 `& l) X1 {/ X2 s& k; L( p' A, w
24. scalar k2
    8 p$ K* z& G' ]+ b4 m2 s5 w$ Y( a
25. scalar k3
    8 x/ e5 M5 h  i! H8 L; X( R' S
26. scalar k4
    , S0 m9 U, d5 G6 \, u
27. 
    . j' J+ N: x8 K4 o. P& G- I& i
28. '定义权重的过程量( N) Q8 W- I4 o\" {5 t# W* u6 d, D
    
29. scalar w1% D! |& N- d; R8 N
    
30. scalar w2
    ) p: m6 f6 T+ `\" b0 G+ K' Q+ i
31. scalar w3# }  U: ~% ]& }: d
    
32. scalar w4# x& W$ F8 s3 _1 W2 c& G- w: |
    
33. 9 q- Q+ p# a$ x; R7 Q0 N
    
34. '程序主体2 K3 j\" O, f! _5 j! D\" i& D! r
    
35. for !k=1 to M step 1
    6 a3 P4 x/ S- V; h& w! J: V
36.   F(!k+1,1)=F(!k,1)+h6 d3 }7 k7 z\" `) r
    
37.   '调用常微分方程计算权重# x6 u2 E% e$ C
    
38.   call obj(w1,F(!k,1),F(!k,2))
    7 v/ p' y6 |' S  L# r
39.     k1=w1*h& a4 s$ y  r0 Q. S1 V- M. `& b5 F
    
40.   call obj(w2,F(!k,1)+h/2,F(!k,2)+k1/2)/ {; X, D( ^$ Y
    
41.     k2=w2*h
    ; l% F! ]# ]! f' l: j! _
42.   call obj(w3,F(!k,1)+h/2,F(!k,2)+k2/2)3 M\" ]3 v- `8 l/ w, b, ^
    
43.     k3=w3*h
    * r\" e( j9 ~/ M9 h% ~: R/ D/ Z
44.   call obj(w4,F(!k,1)+h,F(!k,2)+k3)
    5 a/ Z7 Q; a0 n5 V: \2 I\" r
45.     k4=w4*h9 a6 f\" ^) R; \5 F; S) p
    
46.   '计算函数估计值/ ], p- u* {7 \8 X& G
    
47.   F(!k+1,2)=F(!k,2)+(k1+2*k2+2*k3+k4)/6# }* a) T6 z% G% a
    
48.   '计算函数解析值3 ^) u\" O& @1 X/ `\" V: M) M
    
49.   F(!k+1,3)=@exp(-sin(F(!k+1,1)))*F(!k+1,1)
    0 W( |+ l% c' e( h2 V
50. next
    ! z; Z2 L! Y5 H
51. 
    . g: j9 P0 v0 ^6 ^1 C' U. V& ?
52. '显示最终结果/ |5 d! Y\" A9 {. S+ I  y
    
53. freeze F.xyline
    8 @8 g) k0 t8 X\" w2 f
54. freeze F
    : l/ K& x8 o7 z, [
55. 
    % S9 P$ _% H  T: o! U4 Z
56. '定义常微分方程
    \" S1 t5 g\" ]4 N( P& l0 O
57. subroutine obj(scalar dydx,scalar x,scalar y)- L: |0 U6 {: P; n, ?9 ^
    
58.   dydx=-y*@cos(x)+@exp(-@sin(x))
    9 c! `: d! T( p7 u
59. endsub1 {) _  `' v  p$ J6 K
    

复制代码

运行后求得结果如下：
* n. z$ q5 f; k2 [: Q! C( l

6 R2 i* `( i( M1 e
其中C2列是数值解，C3列是解析解，比较之下，这二者之间无明显差异。
" _* T! K/ _# Q* A2 U( N/ ^( \! Q

& S8 Q/ d( h8 R4 c' ^# z
* Y) M. z% U* Z7 X% e% C
/ ^; R- i: q% Z: X* j