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

liwenhui

EViews除了能解决计量经济学的估计问题以外，还提供一个编程环境用以解决复杂的问题。经过调试，我在EViews中实现了用龙格库塔方法求解常微分方程的数值解，供大家交流。
; v' Y' {8 W3 b( M8 U0 F演示中，我使用了如下常微分方程作为测试：

( [( ^/ a4 O. S+ \' f* v( `这个方程的解析通解是：
+ @- H8 m1 P% I, [

9 V+ z3 {6 J/ }/ ^- E
使用“龙格库塔方法”，编制的EViews程序如下：

1. '用龙格库塔法求解常微分方程dydx=-y*cos(x)+exp(-sin(x)),y(0)=0在区间[0,10*pi]的数值解3 g9 @2 S4 R  m' w9 d' E
   
2. '已知这个微分方程的解析解 y=exp(-sin(x))*x% n/ r8 a8 A4 a7 }7 I; O
   
3. 
   & Z1 V\" t0 Z* Y8 C# M7 Y
4. '生成一个workfile作为基本的数据容器
   3 x, O4 k1 R# v$ s7 a
5. wfcreate (wf=temp) u 1000
   , q# z; }, q5 ?4 D7 U  v! _6 ^
6. + R$ P. @% E3 v( @* p9 ?
   
7. '定义常量. A4 D3 ~0 ?. v. G' o$ C9 M
   
8. scalar pi=3.14159. X\" U- b3 _2 S+ t5 T
   
9. scalar a=0        '定义自变量下限
   7 l& _# h- q, X
10. scalar b=10*pi     '定义自变量上限
    4 e7 ~+ x( Y; m: h: g( G, k# V, Z0 w
11. scalar  M=500       '定义步数
    ( F1 N* [4 m- E3 x# E' D8 `
12. scalar h=(b-a)/M   '计算每步之间的间隔
    2 G2 r0 H8 S5 N# `$ [4 b+ @' V. z
13. 3 Z9 P' \0 L. r
    
14. '定义一个矩阵来储存计算数据，其中第一列储存自变量数据，第二列储存因变量数据,第三列储存解析解的值用以作为比较
    ) c, F& F, @2 _$ V, G: _$ ^
15. matrix(M+1,3) F
    $ C9 i4 `9 D  [2 P, x1 U/ a4 S0 L
16. , n$ ]% S7 l% m( u6 U5 t9 x, n
    
17. '矩阵的第一行储存初值问题的初始条件
    6 f+ u6 t7 Q: Z
18. F(1,1)=06 D1 A; R( i9 b- |! \- _
    
19. F(1,2)=0
    ! k. l4 C6 l: |  ^) B! j
20. F(1,3)=@exp(-sin(F(1,1)))*F(1,1)  p8 k+ r  N' H+ W\" u( l
    
21. $ j- O' }) r; I; |  N+ }
    
22. '定义龙格库塔法的权重参数
    ! u) L$ O3 [, K! X3 {
23. scalar k1
    9 m7 R. z\" B& e9 [( k
24. scalar k29 z8 K+ F& A; ~1 c  n
    
25. scalar k3
    \" d\" A3 A* S. r& R5 S\" X+ p. W
26. scalar k41 Z( P3 G\" z7 i& x0 u: o$ u8 u% A
    
27. 0 S- o7 c2 n0 [7 [7 _
    
28. '定义权重的过程量$ d( Z# `* F/ v1 o: ]
    
29. scalar w1  c* B+ B0 I! i
    
30. scalar w2
    ' x/ s# _. s# E: r) ^& A
31. scalar w3
    + o& l\" [+ t0 m( A. e1 B1 I
32. scalar w4% ~# X1 N( f- y9 D0 J8 o
    
33. 1 H: H1 r# q; F. d
    
34. '程序主体* F$ ~9 ?  F- g# _+ v. M6 s& K
    
35. for !k=1 to M step 1) `9 I9 g$ I2 V5 E5 D
    
36.   F(!k+1,1)=F(!k,1)+h0 X: n\" Y) L6 A5 n. Q( a; i( x3 ^/ H5 W
    
37.   '调用常微分方程计算权重7 h: r: s& {# O4 s) d/ R* s* a2 Q
    
38.   call obj(w1,F(!k,1),F(!k,2)) / O5 T+ u% F% B( u4 y* {
    
39.     k1=w1*h
    0 Y5 l* {4 g( S; \$ q
40.   call obj(w2,F(!k,1)+h/2,F(!k,2)+k1/2)
    8 x3 |$ o1 `6 f) r3 N
41.     k2=w2*h
    8 N6 z  Q5 N* s% E8 P7 j$ o; N
42.   call obj(w3,F(!k,1)+h/2,F(!k,2)+k2/2)8 x0 z6 g$ e- ^5 A
    
43.     k3=w3*h
    \" q/ J- J4 @5 o  e& N( j
44.   call obj(w4,F(!k,1)+h,F(!k,2)+k3)# D- ]  }7 w/ i( s
    
45.     k4=w4*h
    5 z7 M: l7 J' r- L# ?3 [5 }& m
46.   '计算函数估计值6 W1 `7 S: g9 Q* ^
    
47.   F(!k+1,2)=F(!k,2)+(k1+2*k2+2*k3+k4)/69 j1 \3 T% d- u( ~5 H4 t) _. |
    
48.   '计算函数解析值4 g1 K, b/ R7 P( [
    
49.   F(!k+1,3)=@exp(-sin(F(!k+1,1)))*F(!k+1,1)
    3 L' c7 ~9 W9 Z% C; U
50. next
    ) x- k0 g; b' s2 c  j0 z
51. 
    - N- q$ Y4 @3 V- ^
52. '显示最终结果: X) n! b! p5 i) ~
    
53. freeze F.xyline  h1 [+ p: f9 k5 `& L, `- Q
    
54. freeze F
      H& I5 U\" f% A1 w5 p1 S, x
55. 
    # c+ k- X. t* S4 \1 r( ~
56. '定义常微分方程
    4 x$ E- n+ q5 P- O- ^- n7 D
57. subroutine obj(scalar dydx,scalar x,scalar y)
    - N1 O$ d0 g( P\" ]1 j
58.   dydx=-y*@cos(x)+@exp(-@sin(x))
    + o% `5 w6 J0 R7 a\" ]) G! `7 l
59. endsub
    3 ^4 Z. W4 J& Y+ h% e

复制代码

运行后求得结果如下：

2 [: N& }& J- {, _其中C2列是数值解，C3列是解析解，比较之下，这二者之间无明显差异。

+ F4 ^2 _, Y; S5 i0 `3 A4 A
" m  ?4 \/ ~3 R+ Y' E: l7 I

% x( b: a' r$ R: L