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

liwenhui

6 i& x, w0 ^: D( o7 E, D+ x- q. TEViews除了能解决计量经济学的估计问题以外，还提供一个编程环境用以解决复杂的问题。经过调试，我在EViews中实现了用龙格库塔方法求解常微分方程的数值解，供大家交流。
' @. F: P9 @; T( D) y演示中，我使用了如下常微分方程作为测试：
9 z& z1 w" m; H8 @# Y$ v

这个方程的解析通解是：
- z2 a; [* r& J) I" ?% i. M$ u/ O

! A" F9 ^. d- X
使用“龙格库塔方法”，编制的EViews程序如下：

1. '用龙格库塔法求解常微分方程dydx=-y*cos(x)+exp(-sin(x)),y(0)=0在区间[0,10*pi]的数值解' @% c# D! J2 G) ~
   
2. '已知这个微分方程的解析解 y=exp(-sin(x))*x( Z) p  x, V( w4 K6 A% r6 R
   
3. 
   ' i0 v9 |8 @2 m$ L+ }
4. '生成一个workfile作为基本的数据容器4 m  {, |3 X! E  z& ?
   
5. wfcreate (wf=temp) u 1000/ D* t6 ?' I/ I% [; T5 L( r. O
   
6. ( j% V% k5 b9 @; Z5 O
   
7. '定义常量
   4 L4 C0 S- _+ ~1 ~( A- N0 a
8. scalar pi=3.14159
   $ \\" e' r& `- O
9. scalar a=0        '定义自变量下限
   5 n. d& n3 E6 R( Q% T9 G$ B1 y! t) E
10. scalar b=10*pi     '定义自变量上限
    - p  d+ x6 b! V' f) D4 g$ ~9 |
11. scalar  M=500       '定义步数1 l; P- `, q( `2 C9 ]1 d7 G% |
    
12. scalar h=(b-a)/M   '计算每步之间的间隔
    $ z) R& d/ ?# }\" W
13. * s5 ^/ o9 ~\" L  c3 L' L- Z\" Y
    
14. '定义一个矩阵来储存计算数据，其中第一列储存自变量数据，第二列储存因变量数据,第三列储存解析解的值用以作为比较
    % F7 Y; b- N1 Y1 h4 C; [
15. matrix(M+1,3) F % E* g. v/ A5 i4 |  j2 M
    
16. . W( Y, i' j: \6 l+ T: G5 k/ A- j
    
17. '矩阵的第一行储存初值问题的初始条件. {  ^5 F* \1 |* c
    
18. F(1,1)=0  x( a0 s  J5 J; L: Z5 j8 t; ^
    
19. F(1,2)=0
    8 A. b+ J9 m  b, R
20. F(1,3)=@exp(-sin(F(1,1)))*F(1,1)
    6 Z) n/ M3 l; M- ?7 f
21. ; o! X* b$ \7 T+ T0 o
    
22. '定义龙格库塔法的权重参数. E9 e: g* \0 ]. X3 \9 P\" S
    
23. scalar k1
    . q+ u& i9 r1 G9 k* _. k1 X$ S
24. scalar k2
    - s) ~6 s1 s6 k+ c* `
25. scalar k3
    , N7 ~$ m7 ]! ^. ~
26. scalar k4
    5 G* K0 X* C/ c
27. 
    + c! T7 G- l/ R- O. j
28. '定义权重的过程量
    2 z6 k3 }% D  }$ m& D% F
29. scalar w1
    3 C6 G; a\" R. L  p0 r$ W  ?/ z/ q
30. scalar w2
    ' B$ Z- j, O) V1 ]1 y' W
31. scalar w3) h) c) W% z$ c$ s+ ~
    
32. scalar w4
    2 q/ s; T# T3 B: {  b) Y' ?6 O
33. 
    * p0 E/ V9 d( M& f* n; R
34. '程序主体, M\" d% Z5 b# J6 ^( m* r) }# X
    
35. for !k=1 to M step 1. M# \3 J4 ?5 [; T
    
36.   F(!k+1,1)=F(!k,1)+h, I7 E  q\" B( W  @7 d
    
37.   '调用常微分方程计算权重5 z\" z' }0 ?& s( U2 N
    
38.   call obj(w1,F(!k,1),F(!k,2)) , _. I4 }! J2 D4 h) }+ b
    
39.     k1=w1*h
    . R; i4 R9 a; Z7 _\" @
40.   call obj(w2,F(!k,1)+h/2,F(!k,2)+k1/2); P: d* g- n8 j* C\" h
    
41.     k2=w2*h- X3 \+ j8 a& }& H\" x
    
42.   call obj(w3,F(!k,1)+h/2,F(!k,2)+k2/2): L\" l7 b1 B' \, c/ [! p
    
43.     k3=w3*h3 b* S. G, e1 D# D! h: ?% A
    
44.   call obj(w4,F(!k,1)+h,F(!k,2)+k3)- F3 r( x& a/ f
    
45.     k4=w4*h\" y: a. g$ `/ g
    
46.   '计算函数估计值
    ' f% h; W\" }, U9 z4 s3 N
47.   F(!k+1,2)=F(!k,2)+(k1+2*k2+2*k3+k4)/6$ P1 ~4 V$ x- }7 a, L
    
48.   '计算函数解析值
    / b: F- o\" i. w6 H
49.   F(!k+1,3)=@exp(-sin(F(!k+1,1)))*F(!k+1,1)
    + E: q) r+ i0 o
50. next3 N: w! a, r1 t0 a0 F0 q- a& z
    
51. 
    4 I4 O) C. v7 q2 Q$ u+ \
52. '显示最终结果
    4 |0 q+ E  [  V. q/ D& e: D
53. freeze F.xyline) q* G! d2 L1 y
    
54. freeze F
    0 I( e5 J$ C* H) C5 V3 u
55. ! V6 I% t4 _; s0 W
    
56. '定义常微分方程
    8 e5 G& s, d- b9 F* S* z
57. subroutine obj(scalar dydx,scalar x,scalar y)
    ) o; V7 e) Q- b5 r& C# z4 ]/ U
58.   dydx=-y*@cos(x)+@exp(-@sin(x))- ^/ y3 x2 Y' I
    
59. endsub) d1 l: q) R* Z, `8 e* s
    

复制代码

运行后求得结果如下：
& Y$ v( F4 ~4 M: @1 }- P) W
! N- P, p5 Y8 V. f/ n  s6 [8 O- @

其中C2列是数值解，C3列是解析解，比较之下，这二者之间无明显差异。

. e  H* H, d4 R5 Q- t
) Z! U$ w  a6 w
) m" H. x; j# t8 k