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

liwenhui

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

" _( Y  V2 e' l: x2 B+ [
% i6 I4 r* f# ^( X5 ?6 a这个方程的解析通解是：

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

1. '用龙格库塔法求解常微分方程dydx=-y*cos(x)+exp(-sin(x)),y(0)=0在区间[0,10*pi]的数值解1 m+ K4 d7 p. a9 z+ J
   
2. '已知这个微分方程的解析解 y=exp(-sin(x))*x
   . t! m4 i; g: w$ d: T! N6 B
3. & t; }5 `8 r5 y. d
   
4. '生成一个workfile作为基本的数据容器
   1 E8 R: w  X' r: c# I6 E
5. wfcreate (wf=temp) u 1000. ?# `8 I1 e- b$ Z, I
   
6. 
   0 N3 T' q# i/ ^$ s+ W
7. '定义常量% b! x8 [% {) Z1 S) }  q% _3 j
   
8. scalar pi=3.141591 ]0 v+ z/ u. Q/ ^
   
9. scalar a=0        '定义自变量下限
   0 h: M4 z' X% _
10. scalar b=10*pi     '定义自变量上限
    8 c9 X8 Z6 E! E- g6 |1 X4 d- S
11. scalar  M=500       '定义步数
    2 X1 u& C2 t; x7 e
12. scalar h=(b-a)/M   '计算每步之间的间隔
    0 a. N4 x3 Y  [( `6 p9 U4 o$ k
13. 
    0 v9 t) m\" W4 @6 p( X1 ^
14. '定义一个矩阵来储存计算数据，其中第一列储存自变量数据，第二列储存因变量数据,第三列储存解析解的值用以作为比较
    7 m# a9 B8 l: ^$ x6 K
15. matrix(M+1,3) F 9 l' L* \  t3 ]4 L; Z) ~  \
    
16. # ]4 H) K! y* Y2 Q0 G* @
    
17. '矩阵的第一行储存初值问题的初始条件
    1 d9 P; _/ U/ Q
18. F(1,1)=0
    8 X0 z8 ~5 ~7 A5 ~
19. F(1,2)=09 N: E0 d6 z4 k3 j
    
20. F(1,3)=@exp(-sin(F(1,1)))*F(1,1)
    & o/ p( w  v1 Q+ c
21. 
    ! A. I( g9 z9 V# e# X
22. '定义龙格库塔法的权重参数
    8 B. r/ `* p+ H- _\" `
23. scalar k1; t4 \& j, A: v9 @! d
    
24. scalar k2
    / z1 k\" d# T; b1 v  N% O* n8 I
25. scalar k3  W+ l$ T# h% E5 D6 |
    
26. scalar k45 p; X: k8 y3 U% L, `) n) `' V
    
27. $ Y/ e( j& F5 Q$ r0 L, @3 [\" d1 ^! z
    
28. '定义权重的过程量
    \" k& L. Z\" `\" O+ Y+ e
29. scalar w1+ V/ P: a0 g. M+ }  j
    
30. scalar w2
    5 Z6 X7 ]) k) v- P; K: |
31. scalar w3
    1 K8 ~2 p+ Z' f3 d+ I8 L
32. scalar w4
    $ |6 c0 K5 a, B9 J
33. 
      P3 ^+ ]* p/ U: S. `& h
34. '程序主体
    $ L3 `2 D/ X. h& i( w  S2 n( t
35. for !k=1 to M step 1
    . a  ^* F+ R, e6 \\" u
36.   F(!k+1,1)=F(!k,1)+h) R# J' p7 Q' e6 W3 o8 h
    
37.   '调用常微分方程计算权重
    % I+ n; R1 t. \6 k' ~+ `
38.   call obj(w1,F(!k,1),F(!k,2)) ( O3 h' C3 P3 x$ |$ e
    
39.     k1=w1*h
      T0 G( T, ^1 v! F! [. J/ ~
40.   call obj(w2,F(!k,1)+h/2,F(!k,2)+k1/2)$ y& J% B9 E$ i( n# G* s) |5 r9 y9 _
    
41.     k2=w2*h
    1 J( O$ H9 c4 s$ ?5 L- g, j. l
42.   call obj(w3,F(!k,1)+h/2,F(!k,2)+k2/2)
    \" g3 d\" |& ]. M; A7 K6 ], g
43.     k3=w3*h, t$ F/ ?# {  T3 S9 W% [  Z
    
44.   call obj(w4,F(!k,1)+h,F(!k,2)+k3)6 \, R1 N1 c; t; t, w4 C) @/ d
    
45.     k4=w4*h- n8 h3 ~. P) B) X8 t
    
46.   '计算函数估计值4 c- A) _: A\" @. A: N3 a
    
47.   F(!k+1,2)=F(!k,2)+(k1+2*k2+2*k3+k4)/6
    ( T: P( m/ ]) ~' b# N/ g
48.   '计算函数解析值% X# L1 H7 q4 D\" m% p
    
49.   F(!k+1,3)=@exp(-sin(F(!k+1,1)))*F(!k+1,1)' r* M5 y+ Z2 x
    
50. next
    \" A9 L- W- n' \: x2 B. Q
51. 
    3 m9 L: g& S2 J; [' D; S- F
52. '显示最终结果
    & |/ R/ u% j! ^% A# H
53. freeze F.xyline
    1 D0 u: M6 ~' [$ ?* `& d
54. freeze F9 D' m# D) q\" @; j
    
55. # }  e/ v6 i0 n7 N; \$ B+ n
    
56. '定义常微分方程2 N5 S* J7 Z3 o2 e/ {
    
57. subroutine obj(scalar dydx,scalar x,scalar y)
    & n. g  p( b: Y' s! }2 B
58.   dydx=-y*@cos(x)+@exp(-@sin(x))
    ! d* S0 F' T\" @
59. endsub
    * |3 A; r) \2 R# ~\" q. y3 R

复制代码

运行后求得结果如下：
) @) g& A% f. A

5 n* U6 ^) w$ n1 H2 O7 a其中C2列是数值解，C3列是解析解，比较之下，这二者之间无明显差异。


+ f$ K1 [6 Y0 {/ }4 s) n6 N2 ?

9 Y* z* g& N) j  |* Z7 s