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

liwenhui

6 C5 E0 D4 t8 O$ ?/ T$ \4 m& M- o$ T3 MEViews除了能解决计量经济学的估计问题以外，还提供一个编程环境用以解决复杂的问题。经过调试，我在EViews中实现了用龙格库塔方法求解常微分方程的数值解，供大家交流。
* l  x3 F! k3 s2 c演示中，我使用了如下常微分方程作为测试：
; H1 j9 t. y% f# D5 @- P' b8 Q1 f

, E+ j+ _- K0 \这个方程的解析通解是：

/ J, p1 B; o5 R6 W5 I/ }; b$ r使用“龙格库塔方法”，编制的EViews程序如下：

1. '用龙格库塔法求解常微分方程dydx=-y*cos(x)+exp(-sin(x)),y(0)=0在区间[0,10*pi]的数值解
   9 x8 o\" h. E& \/ m
2. '已知这个微分方程的解析解 y=exp(-sin(x))*x$ t' P3 h$ V$ o, g8 H
   
3. 
   ; Y/ [: H2 N$ n\" B' c1 z7 w
4. '生成一个workfile作为基本的数据容器
   8 ~/ _5 Z' t4 y7 \  @7 v7 M( `
5. wfcreate (wf=temp) u 10001 {( K+ {2 }* N$ ?
   
6. 
   ! `) w9 n9 z( F4 g  F/ t' U
7. '定义常量4 H! e5 U, w' X6 z
   
8. scalar pi=3.14159- E/ d6 H5 M4 b) T# H3 J2 \
   
9. scalar a=0        '定义自变量下限
   ( Z' o$ s- l( _5 [7 r4 L
10. scalar b=10*pi     '定义自变量上限
    - Q( r. P6 }% I3 D7 d/ B' H
11. scalar  M=500       '定义步数3 O( x0 ]( ^\" c! z/ D6 _
    
12. scalar h=(b-a)/M   '计算每步之间的间隔2 b1 ~3 Y6 d. g1 a5 F- U( j) H$ d
    
13. 
    : s. Z/ V5 _0 K$ y# Z
14. '定义一个矩阵来储存计算数据，其中第一列储存自变量数据，第二列储存因变量数据,第三列储存解析解的值用以作为比较- y& U( y3 U5 [: k; o3 i0 s
    
15. matrix(M+1,3) F
    3 X$ E& y$ {; i- Z( \( ?
16. 1 c$ d0 ~  Q. d! }- O8 v: y
    
17. '矩阵的第一行储存初值问题的初始条件% h$ y; |' I9 I; p' L! v+ y
    
18. F(1,1)=0
    # ]. z6 l- n' s/ ^
19. F(1,2)=01 V% w3 e( A7 p
    
20. F(1,3)=@exp(-sin(F(1,1)))*F(1,1)
    , \5 g+ i8 \# N
21. 
    + [  {% _) t' ~8 }# |3 t4 r7 z
22. '定义龙格库塔法的权重参数3 R+ Z; J) X\" w2 s
    
23. scalar k1( i, G$ ?2 A  O- h
    
24. scalar k2\" Y/ g# s2 M) e$ w: ~0 R& o
    
25. scalar k3
    & V( X4 [* m0 V& B/ u
26. scalar k4
    + T5 g7 ~. Z6 G1 h: p- G
27. 
    6 Y0 Y9 J* r3 \0 Z
28. '定义权重的过程量$ I& U- ]( n$ {: {- H5 Q
    
29. scalar w1
    & {, n7 w7 d! e
30. scalar w2
    7 Q+ b0 `) T6 z+ u7 x' u8 o- X
31. scalar w3/ [7 {( a! x5 F: r4 F2 l  l& u: n
    
32. scalar w4
    2 t5 j: c0 N! Z( @, Y1 l3 s9 m! O
33. 
    5 ^- q; W! ^- E- e2 w7 x
34. '程序主体
    9 y5 ]7 R; b' e' O# G/ u7 Q5 m+ ?9 @
35. for !k=1 to M step 1- S. l: ^; \\" c# b, |- H
    
36.   F(!k+1,1)=F(!k,1)+h: F# i+ l% Z! k/ R\" v' o
    
37.   '调用常微分方程计算权重# {7 w1 L+ [. {6 `. Y1 f: ^
    
38.   call obj(w1,F(!k,1),F(!k,2)) ( D  Y3 a9 H3 r% s* c! Q8 |
    
39.     k1=w1*h
      w: p6 _1 E3 G# H  h( h# W
40.   call obj(w2,F(!k,1)+h/2,F(!k,2)+k1/2)
    ; J1 F7 `$ e4 L: B' O# Q
41.     k2=w2*h! n2 x. X3 j+ u) l; d: v& S5 S
    
42.   call obj(w3,F(!k,1)+h/2,F(!k,2)+k2/2)* S1 T7 n6 a. {- h. n
    
43.     k3=w3*h\" V' A0 e6 R. T. L+ J5 F* f
    
44.   call obj(w4,F(!k,1)+h,F(!k,2)+k3), L( f8 G1 ~8 X- |
    
45.     k4=w4*h
    + u3 {9 l3 N/ V
46.   '计算函数估计值
    7 A( X4 B( u+ P' Y) U
47.   F(!k+1,2)=F(!k,2)+(k1+2*k2+2*k3+k4)/6
    + ?6 ]\" g6 L2 g$ k% Z
48.   '计算函数解析值
    9 v) T, `+ E! Y7 T2 E
49.   F(!k+1,3)=@exp(-sin(F(!k+1,1)))*F(!k+1,1), L6 L/ k& N/ S3 Q7 C! W
    
50. next& ~9 r: k/ f$ I1 N
    
51. / _2 q' I' O* ~6 L: V* M- g
    
52. '显示最终结果
    $ {0 g& W. u* a1 Y8 e9 W
53. freeze F.xyline
    . o  F2 l2 I; ]1 O% h* I7 E
54. freeze F
    7 G5 U0 a6 O3 i\" h
55. 
    9 U0 Q9 o: X8 x
56. '定义常微分方程+ }( `. }  p, D& L' T# N/ M, _# u
    
57. subroutine obj(scalar dydx,scalar x,scalar y)
    ( j) \& J2 m- ]) P& }
58.   dydx=-y*@cos(x)+@exp(-@sin(x))' K9 |& E\" z0 ]/ k. d
    
59. endsub. w4 O/ Q: R/ z- O( Z! D, A
    

复制代码

运行后求得结果如下：
! u- n2 C( e* S4 g

( f9 n  u+ d' |& o3 [+ \9 ?3 D其中C2列是数值解，C3列是解析解，比较之下，这二者之间无明显差异。
+ _( T5 V) n, g, G+ W" f
5 C  Z4 o3 l1 U# V, j
0 l/ L; _! B! y: W3 g6 B; O" Z
! f" E1 d* S9 }! b