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

liwenhui

% h$ K% x$ d3 e. V3 A
1 s* B: ^6 ?2 z+ [, o# n% s- ^1 e& GEViews除了能解决计量经济学的估计问题以外，还提供一个编程环境用以解决复杂的问题。经过调试，我在EViews中实现了用龙格库塔方法求解常微分方程的数值解，供大家交流。
6 o6 a" F9 r3 \9 h# q演示中，我使用了如下常微分方程作为测试：
1 s$ r1 G! `! X; e, y* T" K

8 x" _3 d2 g1 U% _/ m; b* L0 {这个方程的解析通解是：

. }" _' {0 v8 |
2 T8 e* o; m; V1 j* s使用“龙格库塔方法”，编制的EViews程序如下：

1. '用龙格库塔法求解常微分方程dydx=-y*cos(x)+exp(-sin(x)),y(0)=0在区间[0,10*pi]的数值解- x$ L\" [1 \( W. Z- {
   
2. '已知这个微分方程的解析解 y=exp(-sin(x))*x
   . L! g# ?, ^1 A2 Z; s( [
3. 2 d$ X3 W& s2 {7 V  x
   
4. '生成一个workfile作为基本的数据容器
   ) r( F# z/ ?( x9 N: p* G
5. wfcreate (wf=temp) u 10001 [5 W7 Z* ^! {) s/ S
   
6. 
   6 x3 j2 l8 [1 k6 R
7. '定义常量, ]5 r7 _/ w/ E6 B* S, g
   
8. scalar pi=3.14159
   ( w: F# `5 c5 ^9 j
9. scalar a=0        '定义自变量下限/ V7 C/ D: A4 E# q
   
10. scalar b=10*pi     '定义自变量上限
    4 q; {9 M9 ]# _
11. scalar  M=500       '定义步数
    % H* [; d; H; o9 F* ]
12. scalar h=(b-a)/M   '计算每步之间的间隔: y/ r4 f+ ?% |: z+ s* a' ]0 B
    
13. 
    $ q4 H3 T; }: a% ?' `5 y) P9 U8 i
14. '定义一个矩阵来储存计算数据，其中第一列储存自变量数据，第二列储存因变量数据,第三列储存解析解的值用以作为比较
    4 M/ s: F7 t0 r5 E9 U7 M1 p
15. matrix(M+1,3) F ( _2 `% |* k  Y5 d
    
16. 
    / j! ?9 p0 ~. ]\" k
17. '矩阵的第一行储存初值问题的初始条件
    # P5 k& J1 j5 v1 N# |\" ~( B* R3 w
18. F(1,1)=0# F6 }1 ?+ Z/ l) Y
    
19. F(1,2)=0% G% F- @) D. C\" _0 F7 |$ [4 X, h
    
20. F(1,3)=@exp(-sin(F(1,1)))*F(1,1)% H4 N\" t7 F/ n: o- o4 D& @% w
    
21. 
    - K+ w! L\" o8 o0 b1 y
22. '定义龙格库塔法的权重参数
    4 x7 x2 v4 u8 M, K* W. A
23. scalar k1
    # Q& o6 G2 g: ~
24. scalar k2' S' h& W- R5 r' O
    
25. scalar k3
    $ [. k3 h9 Z6 K
26. scalar k4' d1 L: c2 ~, U( e  y\" w
    
27. \" s\" c& E6 d& ]6 a' {; x# E' t
    
28. '定义权重的过程量\" p& V1 m3 d% K, q\" M4 u
    
29. scalar w1% b6 F- U& J: d3 d\" K
    
30. scalar w2
    ' K( P% E! G3 G- `
31. scalar w3
    + H- ?* ~6 ]\" c* s6 f; d( l
32. scalar w4$ U4 n' W, n, f- }
    
33. 
    4 E2 ]\" z% d' A: i* ~
34. '程序主体0 C8 Q1 [  [5 h2 Y
    
35. for !k=1 to M step 1' j\" U! s; `( A. _
    
36.   F(!k+1,1)=F(!k,1)+h
    ( m& M+ h\" T  V( c+ ^/ W/ ~# k
37.   '调用常微分方程计算权重8 P& i1 Q1 K8 ~6 r# e2 m1 @
    
38.   call obj(w1,F(!k,1),F(!k,2)) # g2 C# p/ u; G2 n- _' J9 y6 `6 H. G4 @
    
39.     k1=w1*h
    + y5 g7 d: b! y& O: |1 k' e
40.   call obj(w2,F(!k,1)+h/2,F(!k,2)+k1/2)7 K: E; \5 q; d3 v! [: |
    
41.     k2=w2*h4 o5 z& d  J7 \# O6 @! i8 }% L
    
42.   call obj(w3,F(!k,1)+h/2,F(!k,2)+k2/2)
    8 E5 k6 A& L* I& B6 d
43.     k3=w3*h0 W9 |$ n% B+ R! y- [8 k& Y0 a$ s
    
44.   call obj(w4,F(!k,1)+h,F(!k,2)+k3); y) ?7 O0 {2 x) o0 u+ Y9 g
    
45.     k4=w4*h
    & Q1 U9 y% O4 h$ M( V8 r
46.   '计算函数估计值+ h/ y2 L9 V# l) `* m' d
    
47.   F(!k+1,2)=F(!k,2)+(k1+2*k2+2*k3+k4)/6
    5 h/ A; L% w\" R1 y! z$ t
48.   '计算函数解析值
    ; X% E3 L0 `& G$ P# }4 P7 ^
49.   F(!k+1,3)=@exp(-sin(F(!k+1,1)))*F(!k+1,1): z! _$ x6 W3 U) ^
    
50. next: k& r7 n9 h/ `/ G
    
51. 5 i3 f; I' V6 ?: J\" f. B
    
52. '显示最终结果
    ) W$ Y# D: Q+ p2 {+ X
53. freeze F.xyline
    0 l2 V2 V/ r) x7 T- N
54. freeze F: M* E/ ~1 ^5 v& \4 y\" a8 }
    
55. 
    : G7 r$ p! K4 u# r2 f% X
56. '定义常微分方程
      F  A\" B3 A) i) F7 w
57. subroutine obj(scalar dydx,scalar x,scalar y)- A0 \8 n& ]& F6 G\" p) n3 ?
    
58.   dydx=-y*@cos(x)+@exp(-@sin(x))3 s& {$ ?# n9 p3 }: Y- z6 g! n
    
59. endsub
    ) Y6 Q9 c% G# D; E1 d

复制代码

运行后求得结果如下：

- d' p. x6 W4 P& s' Y% M; n: N9 r其中C2列是数值解，C3列是解析解，比较之下，这二者之间无明显差异。
. }7 X" B) r* n, u  O8 C

5 S; s$ v# J1 q+ V6 ~) }+ y+ l
5 n5 ^6 q% D( E  D6 K2 E' S