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

liwenhui

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

2 Y' S) f6 w/ v+ W! k: q$ w
2 Y+ y9 a2 [! M) W* c5 i这个方程的解析通解是：
; `2 N/ X+ n" p* s2 u

: c, s: _# b/ J使用“龙格库塔方法”，编制的EViews程序如下：

1. '用龙格库塔法求解常微分方程dydx=-y*cos(x)+exp(-sin(x)),y(0)=0在区间[0,10*pi]的数值解2 v/ f6 [9 @+ e$ A1 k/ Q
   
2. '已知这个微分方程的解析解 y=exp(-sin(x))*x
   2 L0 |; F$ i! p/ p$ `/ p( z
3. $ ?6 o+ F1 Y. H
   
4. '生成一个workfile作为基本的数据容器& _* a: v( f) m+ {. T: i
   
5. wfcreate (wf=temp) u 1000\" {4 D- o  n& @4 {
   
6. 
   + c4 Q3 x- d; {3 f# x
7. '定义常量: a9 O+ G; g. o; }) a2 z
   
8. scalar pi=3.14159
   9 o: I+ z, X& V; V; V6 P  }
9. scalar a=0        '定义自变量下限
   ' l6 _+ \& G) x/ {
10. scalar b=10*pi     '定义自变量上限- K7 Q) X: o0 `+ f* r( r( x& A5 B
    
11. scalar  M=500       '定义步数
    , x3 f9 e9 a. o! |, R7 o' N8 g
12. scalar h=(b-a)/M   '计算每步之间的间隔; Q- R0 m# X# @# A
    
13. ' ~* A' ?9 e% {( l\" c
    
14. '定义一个矩阵来储存计算数据，其中第一列储存自变量数据，第二列储存因变量数据,第三列储存解析解的值用以作为比较
    \" S' E' {, u& C) S7 a
15. matrix(M+1,3) F . e3 V% m3 e  Z+ V7 Q9 p2 p
    
16. 
    4 ?# L* }. c  t9 W* C; a\" f
17. '矩阵的第一行储存初值问题的初始条件
    ; r) h: d- H1 g3 }3 c/ D4 T' X
18. F(1,1)=05 o1 [2 ~3 l5 n+ `
    
19. F(1,2)=0: [( @% K) A! I9 t6 T
    
20. F(1,3)=@exp(-sin(F(1,1)))*F(1,1)1 _! [  j, ~+ B9 a! p; G- N
    
21. 
    0 s/ _9 y2 S- _/ {  d
22. '定义龙格库塔法的权重参数
    \" R' I3 t# `3 d# N
23. scalar k1
    0 a. w4 ]& o, [3 E7 Q! f
24. scalar k26 |% A9 r# j- |- |( m
    
25. scalar k3* R9 F6 O4 J/ ]; h, p
    
26. scalar k45 `6 k- Q: R: e4 c
    
27.   B/ [' D$ P- H' Q/ K9 v+ E
    
28. '定义权重的过程量: }# y) o! L! \3 Y
    
29. scalar w1) l- y8 q: T0 w% e& q  r
    
30. scalar w20 E+ M# U& o  U9 e0 f3 `
    
31. scalar w3
    0 G; v) E5 l6 ]: O' Q\" |% r
32. scalar w4
    - z' Z0 k- n% O) z& E+ B\" m& }
33. \" b, G3 _% z( ?# n2 O
    
34. '程序主体
    # C\" t5 {& [+ f
35. for !k=1 to M step 1
      W# q, r; n4 d8 b: U# e7 r: A6 ]% i; j* g
36.   F(!k+1,1)=F(!k,1)+h) C8 E% F) D' p\" n# j+ H: u& \
    
37.   '调用常微分方程计算权重
    4 D3 z# c2 L/ Q3 y& G8 Q( X% Y
38.   call obj(w1,F(!k,1),F(!k,2))
    0 G1 o% q: }3 D$ K
39.     k1=w1*h% a  z& Y$ T, }: t
    
40.   call obj(w2,F(!k,1)+h/2,F(!k,2)+k1/2)) ]) v5 F# ~. {4 @# p
    
41.     k2=w2*h0 t$ d* ^- f3 O- H
    
42.   call obj(w3,F(!k,1)+h/2,F(!k,2)+k2/2)
    8 R7 f4 Y, m- l8 @/ T( J
43.     k3=w3*h
    \" @, C. `4 y7 H: V& \
44.   call obj(w4,F(!k,1)+h,F(!k,2)+k3)
    . g2 ?/ I. o* z% A0 n
45.     k4=w4*h) Q1 Y7 r3 f! Y$ Z  ~
    
46.   '计算函数估计值
    6 o; A9 y  S, m$ @: x# r
47.   F(!k+1,2)=F(!k,2)+(k1+2*k2+2*k3+k4)/69 r8 N4 K- [' h3 C
    
48.   '计算函数解析值: H\" [\" h  i  I  r  G+ w
    
49.   F(!k+1,3)=@exp(-sin(F(!k+1,1)))*F(!k+1,1)
    # ^3 v! M8 F% s; ~0 Q; F& D7 t1 R% X
50. next
    4 w; s3 g7 n# \5 ^/ o; ]\" w\" F0 l
51. 
    8 [: Z\" ?+ A2 i0 v; R  }# n5 r
52. '显示最终结果
    % `$ i- y; m# ?
53. freeze F.xyline, ~$ \) o- l, b* T! ]
    
54. freeze F7 w+ @: O0 ?0 k  e
    
55. 
    0 D7 H( f\" X. V( {  i9 Q
56. '定义常微分方程& w& C& O8 s2 B$ p' B! c7 I
    
57. subroutine obj(scalar dydx,scalar x,scalar y)7 `8 h- [. J  h4 ?7 U7 [
    
58.   dydx=-y*@cos(x)+@exp(-@sin(x))! M8 `* M+ k# M4 k\" d. V# |
    
59. endsub
    8 T2 Z8 \- s5 @4 O$ I

复制代码

运行后求得结果如下：

2 w9 \0 b2 T7 X6 c* x1 e6 M9 V$ A

% _9 Z$ e3 w5 L其中C2列是数值解，C3列是解析解，比较之下，这二者之间无明显差异。
1 z' F# t/ {  h; o! b

+ ]7 _! t/ o/ c