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

liwenhui

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

7 ^  G0 N. i( F) m' `
- b  P7 g4 g& v+ p这个方程的解析通解是：

使用“龙格库塔方法”，编制的EViews程序如下：

1. '用龙格库塔法求解常微分方程dydx=-y*cos(x)+exp(-sin(x)),y(0)=0在区间[0,10*pi]的数值解2 d# s4 z( Q; j
   
2. '已知这个微分方程的解析解 y=exp(-sin(x))*x' u- k' r2 J( \: o# Z2 R0 E9 {% i3 [
   
3. 9 y. V/ L4 o/ \1 I, \  \
   
4. '生成一个workfile作为基本的数据容器
   \" R, \% H. n; ~8 x- I9 U
5. wfcreate (wf=temp) u 10009 D/ ^/ c+ x& }# M& E' I. z# G
   
6. 1 ?7 u! C& W, B6 w  \' P
   
7. '定义常量+ l/ D1 k8 [6 z5 _
   
8. scalar pi=3.14159
   + \( \* w; V; E. d\" q  G3 G- [7 M
9. scalar a=0        '定义自变量下限9 S# E: t; B& |; h- j: N1 ?1 u
   
10. scalar b=10*pi     '定义自变量上限
    : f3 b) i! D/ N: c. ]7 d
11. scalar  M=500       '定义步数8 j9 R! d, O+ Y: R7 y\" r0 ^
    
12. scalar h=(b-a)/M   '计算每步之间的间隔  b% Z) @+ `( Y4 h, W, C
    
13. 
    3 P8 k\" S4 d% w) Y$ P. ]6 q
14. '定义一个矩阵来储存计算数据，其中第一列储存自变量数据，第二列储存因变量数据,第三列储存解析解的值用以作为比较! t) A5 R+ b5 Q$ Z  V0 A
    
15. matrix(M+1,3) F & L! ~8 V2 W% H; [1 g! {
    
16. 
    / X- A3 f( ?, x
17. '矩阵的第一行储存初值问题的初始条件( E6 ~2 g+ E+ M\" p/ {# a
    
18. F(1,1)=0+ w- L1 v- c3 |) a' `4 E+ I
    
19. F(1,2)=0
    ( `$ a4 G! k& P7 S) Q& s# A
20. F(1,3)=@exp(-sin(F(1,1)))*F(1,1)
      N% M& `8 Z+ d. e: w
21. - Z  w  V, d6 e, j$ f; O
    
22. '定义龙格库塔法的权重参数2 F3 S# x# T5 q# {# n0 ?
    
23. scalar k1
    ) F5 @5 H7 `\" g0 A4 V3 s
24. scalar k2
    . ?2 L$ {! t  L/ ^\" z, B3 Z
25. scalar k3
    ) [6 J% g3 e' l. I) A$ S8 Q- d
26. scalar k4
    0 m! Y; G2 O: a/ R* i6 g1 q/ ~
27. 4 |0 J  O\" n  _! g
    
28. '定义权重的过程量# s- ^8 a\" V3 e: g
    
29. scalar w18 h( ?) B\" `7 M7 _- J2 s8 `
    
30. scalar w20 N: {1 g$ i) [& v+ w+ W. f+ L  t
    
31. scalar w3) e6 Y, _% M9 C8 g1 n
    
32. scalar w4
    1 W4 G+ B- _  ?2 p
33. \" _- m4 r9 h, [\" W6 ?
    
34. '程序主体
    6 `9 o4 \  Z$ ]1 a; F
35. for !k=1 to M step 1; a( C( }& M' m1 u+ k\" a
    
36.   F(!k+1,1)=F(!k,1)+h
    . U\" o3 {: R4 L0 n, Y$ P  u# l
37.   '调用常微分方程计算权重. ^5 I2 _# |/ R) w
    
38.   call obj(w1,F(!k,1),F(!k,2)) 6 j# {- a, ]$ g2 n# A
    
39.     k1=w1*h' `. f% i& s9 T& i& [' F
    
40.   call obj(w2,F(!k,1)+h/2,F(!k,2)+k1/2)/ ]1 p$ j1 p6 S7 ^9 s
    
41.     k2=w2*h5 ]3 L8 {% U! V7 _& d
    
42.   call obj(w3,F(!k,1)+h/2,F(!k,2)+k2/2)& ~' _7 K; h3 g' o, B
    
43.     k3=w3*h! ]  V' M( X7 i$ ~6 q4 N# n3 n
    
44.   call obj(w4,F(!k,1)+h,F(!k,2)+k3)7 f. h% W; c+ q
    
45.     k4=w4*h
      s2 f, i0 D3 s! w9 o  O( ^5 d
46.   '计算函数估计值$ J1 j! {$ _, R0 l6 \
    
47.   F(!k+1,2)=F(!k,2)+(k1+2*k2+2*k3+k4)/6( T& ]! N  ]8 Q/ Z. E& z4 V
    
48.   '计算函数解析值; u9 q! l5 Z' N3 w) }& g9 |) c' ]) L
    
49.   F(!k+1,3)=@exp(-sin(F(!k+1,1)))*F(!k+1,1)
    $ e$ W  y5 c8 f0 U
50. next% t( L# [6 Q6 f: e( q0 r
    
51. 4 S1 l) i- o# [5 Q: S
    
52. '显示最终结果
    \" k$ f* G: u( a( j- Q! N* L
53. freeze F.xyline
    2 `3 s1 Z% A9 P1 r* J
54. freeze F
    7 M3 g\" I; H6 J
55. 
    8 m5 k2 y2 `0 {# w
56. '定义常微分方程& e\" ?, M' o2 ^\" m& _
    
57. subroutine obj(scalar dydx,scalar x,scalar y)
    8 e1 f; e% K$ Y8 ~9 U9 F\" y; }. F
58.   dydx=-y*@cos(x)+@exp(-@sin(x))/ M\" l' F/ {\" s3 e7 h2 b
    
59. endsub
    ! G0 c( |7 i/ }0 L

复制代码

运行后求得结果如下：
8 y. w$ \- X8 s7 k

+ x! ^0 T8 j; r8 a! g
其中C2列是数值解，C3列是解析解，比较之下，这二者之间无明显差异。
' H0 i% J1 L6 z: x8 C6 }
) G* K* S+ D8 O( F0 H# L3 ~# U
7 H! g/ L, G0 Q4 r+ H9 F  i
3 M8 H7 @. }" D
/ c! K" @2 d; J$ b3 @0 @