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

liwenhui

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

$ p( n/ k! i0 a( J$ }1 K这个方程的解析通解是：
8 ]' ?$ S! ]1 o* u

& N$ U, M0 Q% ]- _- z1 V" [
( Z; Z6 ]8 Y$ n使用“龙格库塔方法”，编制的EViews程序如下：

1. '用龙格库塔法求解常微分方程dydx=-y*cos(x)+exp(-sin(x)),y(0)=0在区间[0,10*pi]的数值解
   % P; _; |; ^1 o6 z% m
2. '已知这个微分方程的解析解 y=exp(-sin(x))*x* k  B1 Q( T7 N% j8 a9 B& n
   
3. 
   - r7 v\" A- F8 s/ |5 l
4. '生成一个workfile作为基本的数据容器
   7 ~; T0 Q3 {  e3 R& h
5. wfcreate (wf=temp) u 10008 @7 \* A- X. z. O. c/ b0 P
   
6. 
   - X) ?* h: m- [0 A& |\" v
7. '定义常量
   ( s$ B; N1 K! u, t2 A; M
8. scalar pi=3.14159
   / f6 }* Y: G1 V$ O, l% O
9. scalar a=0        '定义自变量下限
   3 a0 R$ b) y% j* m% S3 S& P
10. scalar b=10*pi     '定义自变量上限
    % X9 W4 \% m  z! H, G1 u
11. scalar  M=500       '定义步数
    6 b! |- n# Z' W! N1 E7 h
12. scalar h=(b-a)/M   '计算每步之间的间隔5 v8 J- H# A5 E( T\" w
    
13. 5 i9 u# D! R% \7 _* b6 U) H
    
14. '定义一个矩阵来储存计算数据，其中第一列储存自变量数据，第二列储存因变量数据,第三列储存解析解的值用以作为比较
    ' H0 }3 b0 V+ s# w
15. matrix(M+1,3) F
    7 z0 `8 e# ~6 V. P. l; I
16. 
    ; m0 a9 ?* Y) Z% o) A# y  l8 Q
17. '矩阵的第一行储存初值问题的初始条件( I# u* s5 v7 _  a3 i
    
18. F(1,1)=0
    : ^5 u3 l* x8 I
19. F(1,2)=0
    1 s7 I; j' ?. R% i4 m& S% g
20. F(1,3)=@exp(-sin(F(1,1)))*F(1,1)0 R* D) H1 @) D' H# L& u
    
21. 
    : O0 j9 |6 s' v0 x\" F  c- z0 Z: y
22. '定义龙格库塔法的权重参数
    ' y9 ^, d\" A# l$ i  u: S; H
23. scalar k1
    - l2 [8 B+ Q9 l4 p' z5 z9 i
24. scalar k23 Z% w; \7 W# S5 Z7 f\" g
    
25. scalar k3
    6 |8 \+ D' M/ }6 j% F  G
26. scalar k4, O+ f. F8 u/ L& [% j
    
27. 
    % a# F7 H+ R1 c9 B3 C
28. '定义权重的过程量+ W7 `; y% @4 ~9 ]( s  {/ w
    
29. scalar w1
    - u8 x% l- b3 f, H
30. scalar w2
    0 l! s. k  K; c+ _9 ]; L
31. scalar w37 o( \4 r; U\" K! s* e4 p' k
    
32. scalar w4
    / E; ?9 s1 _0 x1 f) s) e1 @\" ~
33. 
    , m( J( x# i) {8 U1 i
34. '程序主体- O8 q\" B8 x2 O' f+ k
    
35. for !k=1 to M step 14 O; A+ ^; e& S
    
36.   F(!k+1,1)=F(!k,1)+h6 x) c/ K\" Q0 |( n% O& p; q
    
37.   '调用常微分方程计算权重
    % r3 _# B! V0 F
38.   call obj(w1,F(!k,1),F(!k,2)) 0 p7 Y) d9 [2 L' b+ L
    
39.     k1=w1*h* v) E0 C& M! N: h' ~
    
40.   call obj(w2,F(!k,1)+h/2,F(!k,2)+k1/2)
    * M( ?7 w4 u* y! w; x* U\" `) [
41.     k2=w2*h
    3 P8 Z/ O+ ]6 s( S7 H4 M
42.   call obj(w3,F(!k,1)+h/2,F(!k,2)+k2/2)
    ' ~\" I% A* r) c& g' m0 U
43.     k3=w3*h
    ) @  L& G& A& ^7 X/ e: W
44.   call obj(w4,F(!k,1)+h,F(!k,2)+k3)
    2 m: K0 S5 z$ @6 H/ j1 E
45.     k4=w4*h1 k1 a; c1 N6 D& u5 D! |' a5 \6 r
    
46.   '计算函数估计值
    + C9 y  c5 q( ]  l1 s: b7 R
47.   F(!k+1,2)=F(!k,2)+(k1+2*k2+2*k3+k4)/6
    % }; L9 P5 ]; `5 d
48.   '计算函数解析值
    1 V8 P' ^, B& n  t
49.   F(!k+1,3)=@exp(-sin(F(!k+1,1)))*F(!k+1,1)& |/ r3 J; J( _$ [6 }6 o: `- w, ^
    
50. next6 c( u' {* k1 |. L3 }
    
51. , U/ S& m  r' T$ u8 s
    
52. '显示最终结果) x! F$ {/ S9 Y6 e& {
    
53. freeze F.xyline6 j6 \3 u  d+ M( c5 r9 u
    
54. freeze F
    : G4 {\" L5 D. Q, r  ^) m
55. 
    * t\" ~+ ~. Z- Q$ Y* U' Z! g( e4 v
56. '定义常微分方程7 U4 ?4 [2 N/ F/ \) k# ~6 W! ~! z' q
    
57. subroutine obj(scalar dydx,scalar x,scalar y)
    / b5 y7 l: X' D: E' Q3 r
58.   dydx=-y*@cos(x)+@exp(-@sin(x))! |, [3 D- g  W! O; U
    
59. endsub
    7 x+ l& b0 l/ T; C1 k8 ~+ x9 M

复制代码

运行后求得结果如下：
4 z% w& L$ |- \8 \4 @

1 j. v. z$ G' O. K8 m' U7 ^1 ]! ?其中C2列是数值解，C3列是解析解，比较之下，这二者之间无明显差异。

3 r. y& n. Y$ X8 X% Q) S" Z

' C, p( V6 ~3 V