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

liwenhui

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

6 }* h) }: M$ v
5 d8 h  h# V4 }, b4 J$ q9 U这个方程的解析通解是：

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

1. '用龙格库塔法求解常微分方程dydx=-y*cos(x)+exp(-sin(x)),y(0)=0在区间[0,10*pi]的数值解! ^/ l: P4 n2 `, V
   
2. '已知这个微分方程的解析解 y=exp(-sin(x))*x
   2 Q9 I# b8 |0 n$ u' _
3. \" L4 R. l8 W* T8 L, t
   
4. '生成一个workfile作为基本的数据容器
   & L/ V5 a% A; E1 K1 _! z( k
5. wfcreate (wf=temp) u 1000! |' Q0 \! i* f) u, z0 f4 G6 E) l
   
6. 
   # J& ]( `' n: A2 V) D% H
7. '定义常量, w/ K  a2 V\" P+ B; i2 M
   
8. scalar pi=3.14159
   : i. v% i& J* e' p
9. scalar a=0        '定义自变量下限# S& ^# G$ P\" W! ~
   
10. scalar b=10*pi     '定义自变量上限; m: c/ o6 z( x& L/ |
    
11. scalar  M=500       '定义步数
      z\" |  f1 R) k+ T& k5 a8 u8 M+ u. x
12. scalar h=(b-a)/M   '计算每步之间的间隔; v\" F5 }\" k+ K  a( _
    
13. ! Y  d6 ]9 u9 ~' `4 V\" X( D' V1 E
    
14. '定义一个矩阵来储存计算数据，其中第一列储存自变量数据，第二列储存因变量数据,第三列储存解析解的值用以作为比较) J: k/ x# l& i\" Q0 K. ?
    
15. matrix(M+1,3) F   s8 x: V+ g7 ?- J0 i
    
16. 
    8 l1 T! O1 V! _9 a5 m& w
17. '矩阵的第一行储存初值问题的初始条件
    1 J) u4 W* Y) X: t- i2 [
18. F(1,1)=0$ V# q/ Z+ q1 h. @7 @
    
19. F(1,2)=04 s\" p8 `: R7 n
    
20. F(1,3)=@exp(-sin(F(1,1)))*F(1,1)* q: g4 H9 M& u6 B
    
21. , u% f% \8 K+ n! I5 J2 L
    
22. '定义龙格库塔法的权重参数/ Q; S# x9 w/ F* o* [
    
23. scalar k19 V8 a: J8 Q9 Q- \+ s( \0 \
    
24. scalar k2
    * |* D( {' G$ D9 _
25. scalar k3
    $ i! a8 F\" _; P: L& k
26. scalar k4, d\" o; ^! [5 v/ @$ M. Y
    
27. 
    ; s- l/ e' n0 t8 M# Q5 Y
28. '定义权重的过程量  |9 o0 q' W1 W) H$ c' N
    
29. scalar w1
    - y\" t0 n* r7 J& g, F' F6 w
30. scalar w2
    / ]$ W. j: T( c3 W\" p2 n: S
31. scalar w31 T' t, v. Y% ]\" B9 y2 V* G
    
32. scalar w4
    ' g4 T4 ]5 x* \# Z! _1 r
33. 
    3 B) p  {3 ^: J& W: ^
34. '程序主体2 g: S: f: _' R\" g3 r
    
35. for !k=1 to M step 15 y/ v7 L2 w% j
    
36.   F(!k+1,1)=F(!k,1)+h  G\" o% u& v& b) F  N
    
37.   '调用常微分方程计算权重
    2 Y: C2 r1 U4 `- d$ s2 [* |
38.   call obj(w1,F(!k,1),F(!k,2)) : J. B, z& X6 ^9 K+ p+ Q7 ?
    
39.     k1=w1*h
    & X0 u1 j: r3 d4 L! [- w* c5 M' z- J
40.   call obj(w2,F(!k,1)+h/2,F(!k,2)+k1/2)9 r- ]; z: j8 C! R$ {2 T' Z  `
    
41.     k2=w2*h
    ) ^/ ^. t/ m6 z  Y. i' l\" Z. \
42.   call obj(w3,F(!k,1)+h/2,F(!k,2)+k2/2)
    ; ]. k  v& u- }% R. s
43.     k3=w3*h- l! V+ i5 ?6 q3 c6 \( L8 {
    
44.   call obj(w4,F(!k,1)+h,F(!k,2)+k3)$ ?, O. k; U1 F' T9 \% Q
    
45.     k4=w4*h
      {* Z. F) u1 G! A- J
46.   '计算函数估计值
    0 c3 J# ^( w0 Y/ N+ r  H- U: ]
47.   F(!k+1,2)=F(!k,2)+(k1+2*k2+2*k3+k4)/6
    ( o3 O1 C/ N8 k& S/ e
48.   '计算函数解析值4 S) ~8 }; |3 v+ ]4 ]9 _$ y
    
49.   F(!k+1,3)=@exp(-sin(F(!k+1,1)))*F(!k+1,1)
    3 z; P4 D3 l! k7 i) k- ]1 ~
50. next, N4 g) N* A  `* Q8 l' ~
    
51. ( h) z4 I# C  D; t( Q- Z! B+ S* N! v
    
52. '显示最终结果1 G# h9 V2 z, F; x\" J; A8 [2 j1 }
    
53. freeze F.xyline
    . `2 A% v\" o/ y3 V
54. freeze F
    * v( N& ?9 I& p# v& R; Y; W
55. ) ]$ j7 s4 i2 @, |/ k% \
    
56. '定义常微分方程0 M8 w\" u. Z' b: u+ T. S3 |
    
57. subroutine obj(scalar dydx,scalar x,scalar y)
    : A8 F6 j0 ?; ^6 G; E2 s
58.   dydx=-y*@cos(x)+@exp(-@sin(x))
    2 @4 m1 y$ n. T
59. endsub6 x( [& V! t* s+ h7 V( U
    

复制代码

运行后求得结果如下：
6 D; |7 n' i3 k1 s  `
" }# a  e& x3 s0 k" j

2 |* Z# `* ?7 D/ T# i& o, E
8 ]. g, m7 v: y其中C2列是数值解，C3列是解析解，比较之下，这二者之间无明显差异。

, u7 F( s# Q% f/ B
; \( @) d6 Y, x