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

liwenhui

. n' M. Z8 |; T0 @4 Y
" K5 p5 g5 R: W$ SEViews除了能解决计量经济学的估计问题以外，还提供一个编程环境用以解决复杂的问题。经过调试，我在EViews中实现了用龙格库塔方法求解常微分方程的数值解，供大家交流。
. w, d$ q4 _+ D0 C! G演示中，我使用了如下常微分方程作为测试：
, u6 c" i0 Q. x% ]7 H) Q

, ^% d% x$ I2 c1 B- Q
这个方程的解析通解是：

2 Z& k; @0 ]( }/ S/ U. i
使用“龙格库塔方法”，编制的EViews程序如下：

1. '用龙格库塔法求解常微分方程dydx=-y*cos(x)+exp(-sin(x)),y(0)=0在区间[0,10*pi]的数值解5 O% |! r0 I1 N! I9 u
   
2. '已知这个微分方程的解析解 y=exp(-sin(x))*x
   * L: U1 h$ c* D7 J
3. 
   5 e6 T: w# q0 u, D, H
4. '生成一个workfile作为基本的数据容器5 X1 E; R2 K8 n2 o6 j9 S
   
5. wfcreate (wf=temp) u 10007 T# r3 \9 V4 H& V; E- U
   
6. % U  H* m! f0 }( [1 N
   
7. '定义常量0 B8 Z2 l, z  P
   
8. scalar pi=3.14159# C: q/ T+ `; E7 |& j9 D% G  u
   
9. scalar a=0        '定义自变量下限$ b& [- q* ^# ?
   
10. scalar b=10*pi     '定义自变量上限6 X9 |/ [; L, p2 o
    
11. scalar  M=500       '定义步数
    # @( d! r$ p1 j& i% d. s7 X
12. scalar h=(b-a)/M   '计算每步之间的间隔
    + Z& T: k0 p7 X& a
13. 
    . F& I, g\" ^8 J4 r
14. '定义一个矩阵来储存计算数据，其中第一列储存自变量数据，第二列储存因变量数据,第三列储存解析解的值用以作为比较
    7 M' S2 E/ p+ [5 z# p
15. matrix(M+1,3) F
    . [/ L$ d6 b- f/ |* I
16. 
    ) f6 S* R4 _2 p' u9 M
17. '矩阵的第一行储存初值问题的初始条件
    + [\" S3 t# s, S8 z6 P4 o3 ]
18. F(1,1)=0
    3 b% G# @6 C4 S
19. F(1,2)=0
    4 Q\" R! e$ ^( j
20. F(1,3)=@exp(-sin(F(1,1)))*F(1,1)4 P' t' q+ [$ a( \' Q, A
    
21. 
    & Q( X2 I7 A% M, F
22. '定义龙格库塔法的权重参数3 I) A6 J, E& }! @: k8 x7 ]  q
    
23. scalar k1% T9 J1 L9 w; F% ]  M+ p$ i) R2 [) u
    
24. scalar k2# a  \2 w# J$ H; P
    
25. scalar k3
    / [( Q' I4 B7 V
26. scalar k4
    4 X7 i1 l1 P! u6 k
27. 1 G* A# F\" E5 U( K: `
    
28. '定义权重的过程量1 t. R4 s9 h1 R+ y3 @0 @; I( {\" ~
    
29. scalar w15 J' ^5 j1 M\" z+ E; q
    
30. scalar w2
    , _' C  U# h. n) i1 P; e4 C$ O' s+ \
31. scalar w31 c4 C/ G3 E& Z9 p5 M1 v
    
32. scalar w4
    7 _9 o% X! f1 r) b. [1 Y\" v
33. # ]+ C: f4 W: c  k* p& s
    
34. '程序主体
    6 A- w* L4 n* R! A: K
35. for !k=1 to M step 1) P2 J' l  ~- b/ f7 q& M7 i& A
    
36.   F(!k+1,1)=F(!k,1)+h' N# k1 g3 b/ W: o3 Z* Q
    
37.   '调用常微分方程计算权重
      M6 x$ {# y) m$ t6 ~\" X! _
38.   call obj(w1,F(!k,1),F(!k,2))
    , U. T0 t% }- W7 j
39.     k1=w1*h6 l( U( c8 F7 P4 C- ]! ^
    
40.   call obj(w2,F(!k,1)+h/2,F(!k,2)+k1/2)2 C6 `5 E& j; Q. A, s' A) s
    
41.     k2=w2*h- i5 U9 w( o! s; C
    
42.   call obj(w3,F(!k,1)+h/2,F(!k,2)+k2/2)$ c$ j7 ^# m  y# f, s' \  _9 l' I) M
    
43.     k3=w3*h- E/ ]) t4 i. G
    
44.   call obj(w4,F(!k,1)+h,F(!k,2)+k3)- X) O9 _2 f& h3 e' b: C0 n6 {
    
45.     k4=w4*h
    5 A  g0 k/ V% S+ B, Z/ }8 k# x6 x
46.   '计算函数估计值' z  `- ^% d  t\" {' ?\" P0 k
    
47.   F(!k+1,2)=F(!k,2)+(k1+2*k2+2*k3+k4)/6
    2 U3 P1 P  t! {9 l5 z1 ^9 R, g
48.   '计算函数解析值
    + k/ {% z; q6 U
49.   F(!k+1,3)=@exp(-sin(F(!k+1,1)))*F(!k+1,1)
    $ P, _% ^\" A; ]! L8 z+ l! F
50. next& y% ^. B4 I. K: S- g
    
51. * s\" [& Z. o! q( I' Y3 V- ?) Y: j
    
52. '显示最终结果/ U& u4 @* p9 W\" p
    
53. freeze F.xyline2 t1 ^4 U5 n3 U  K1 U
    
54. freeze F
      Y) D! K5 a5 \
55. 
    \" Z\" G, z6 K+ r& {3 v# T\" ?
56. '定义常微分方程
    5 k: q# ~! o3 t* \. d
57. subroutine obj(scalar dydx,scalar x,scalar y)
    ( e- _& z3 ~' @  o
58.   dydx=-y*@cos(x)+@exp(-@sin(x))
    ' X4 w* i* {* ?0 ]7 P/ V3 t$ Y
59. endsub
    ; r  s- h% }5 c5 ]% M  q' S

复制代码

运行后求得结果如下：
( C1 r; e( \  H* w+ `
: `' n" F/ B  N, F

: a8 K6 x% P6 o2 _) F+ S) S其中C2列是数值解，C3列是解析解，比较之下，这二者之间无明显差异。

4 @- a0 W0 V0 d: z0 G
% K# u! ]0 f( o" ~/ m7 {* t5 d9 T