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

liwenhui

EViews除了能解决计量经济学的估计问题以外，还提供一个编程环境用以解决复杂的问题。经过调试，我在EViews中实现了用龙格库塔方法求解常微分方程的数值解，供大家交流。
( J: X0 c: G, Z8 w- B5 @演示中，我使用了如下常微分方程作为测试：
# C: G6 |* M$ ~& b6 }6 k

3 H9 _7 I: O3 ]$ U0 P
这个方程的解析通解是：

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

1. '用龙格库塔法求解常微分方程dydx=-y*cos(x)+exp(-sin(x)),y(0)=0在区间[0,10*pi]的数值解( b) p+ [2 n, z% M2 x
   
2. '已知这个微分方程的解析解 y=exp(-sin(x))*x
   4 b: V( u1 f$ n* U
3.   U  L# z& H. B
   
4. '生成一个workfile作为基本的数据容器
     J$ s% R$ g: P# k2 B8 _
5. wfcreate (wf=temp) u 1000
   6 |) l, \3 d% q3 |, l8 c8 i1 o0 R8 F
6. 
   2 k0 c+ F* c7 a7 \
7. '定义常量
   3 S- E$ ~: s. \2 M% X  X! ?; a
8. scalar pi=3.141595 W% Q7 O  T. J/ t\" @) b3 h5 Y\" M
   
9. scalar a=0        '定义自变量下限
   ) f- g$ ^1 V2 }' n9 I. L
10. scalar b=10*pi     '定义自变量上限
    6 u. [! e) L/ `8 B+ d% v
11. scalar  M=500       '定义步数' S7 _* A6 U- P3 P
    
12. scalar h=(b-a)/M   '计算每步之间的间隔- g8 T9 R& d\" n% B' p1 o3 T
    
13. 
    6 c2 j: G0 o' B2 |9 R) e! d- G
14. '定义一个矩阵来储存计算数据，其中第一列储存自变量数据，第二列储存因变量数据,第三列储存解析解的值用以作为比较. t. P' ]: f3 ~# k- a' q/ _
    
15. matrix(M+1,3) F
    , G7 z( ^. K( X2 T: `
16. 
    ' ~+ _; S2 h& W: ]% m- m
17. '矩阵的第一行储存初值问题的初始条件
    / v; }/ p6 [1 ^' C1 }( O
18. F(1,1)=0
    : i7 |# B+ N9 Q  d6 y
19. F(1,2)=0+ C' e! h5 V( b; _) V- w) w
    
20. F(1,3)=@exp(-sin(F(1,1)))*F(1,1)
      S3 o' f9 X! {8 c3 b# A0 B# @
21. 
    5 D- j0 v$ G' r7 Q* V/ x2 ?
22. '定义龙格库塔法的权重参数
      ~4 A+ C; a/ }9 j
23. scalar k1/ b7 u8 X\" ~2 a9 {% ~1 B
    
24. scalar k24 v! q$ P2 L& H\" _9 T4 m
    
25. scalar k3
    7 I' [' x* e2 Y! h; E6 |
26. scalar k4
    ( N: K( b- K9 e/ r' Y8 a5 O6 X
27. 
    4 A; c0 @2 a& k4 C1 {6 q# i. D
28. '定义权重的过程量- k! b# g( a( u1 U! E+ a
    
29. scalar w1+ j0 I. V1 Q0 t\" \% H. }/ @7 v
    
30. scalar w2' A% C  u! ^0 k8 s& t
    
31. scalar w39 I5 [# Q( ~/ U( L5 I5 i' e/ q
    
32. scalar w4' l0 o1 \9 b; z\" u4 X! K
    
33. 
    5 l# f* D6 g* J2 T2 r# ?0 l
34. '程序主体. C+ e2 G/ t1 k/ }* c; r& k
    
35. for !k=1 to M step 1
    4 ?. ]) T4 A. U+ J- _! P% Q, A
36.   F(!k+1,1)=F(!k,1)+h8 D: j\" U( [4 d) f
    
37.   '调用常微分方程计算权重
    * b+ R3 ^% s# c$ c
38.   call obj(w1,F(!k,1),F(!k,2))
    4 I7 }8 Q$ V- J8 e/ E
39.     k1=w1*h
    & @. h( i3 I$ R7 U
40.   call obj(w2,F(!k,1)+h/2,F(!k,2)+k1/2)7 q2 w0 Z) }  s# D7 _. |
    
41.     k2=w2*h5 _, g! B6 B2 x' I
    
42.   call obj(w3,F(!k,1)+h/2,F(!k,2)+k2/2); ~' l$ `9 t+ {) {- S/ ~- Z
    
43.     k3=w3*h
    1 g, M- c6 N1 F& B6 l7 `
44.   call obj(w4,F(!k,1)+h,F(!k,2)+k3)
    2 t6 H4 i3 h$ t1 @; q- u3 K( R
45.     k4=w4*h% \7 ]7 {0 e: ]' P) E8 Y
    
46.   '计算函数估计值$ x  R+ j) I0 }# V; H
    
47.   F(!k+1,2)=F(!k,2)+(k1+2*k2+2*k3+k4)/6+ Z( d9 W1 }9 }- K4 c9 {- J
    
48.   '计算函数解析值: {& N8 @' `6 ~% F4 |
    
49.   F(!k+1,3)=@exp(-sin(F(!k+1,1)))*F(!k+1,1)
    9 A: B3 R7 @! X7 [9 n$ F9 ]: c) q0 F
50. next
    - q0 o2 ^& S3 x: v
51. 1 u, {# _; ?! B
    
52. '显示最终结果
    1 w1 y7 L  d6 z& w; B
53. freeze F.xyline
    / X3 M2 k4 a; a1 C
54. freeze F9 o  L- P& d7 U( E0 J* i+ [
    
55. 0 ?0 P\" r' A. S. a
    
56. '定义常微分方程+ |, M0 y& J7 N$ w5 u. O
    
57. subroutine obj(scalar dydx,scalar x,scalar y)& Y) V; p* u; C6 |\" J
    
58.   dydx=-y*@cos(x)+@exp(-@sin(x))( U. \* P. ], d! B/ [3 s' G
    
59. endsub) Q5 \4 R6 o5 Z/ ?
    

复制代码

运行后求得结果如下：

1 N  k+ b$ A$ Q# ?( I) ^9 i
其中C2列是数值解，C3列是解析解，比较之下，这二者之间无明显差异。

* F9 s( H4 a8 ?& r1 a