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

liwenhui

EViews除了能解决计量经济学的估计问题以外，还提供一个编程环境用以解决复杂的问题。经过调试，我在EViews中实现了用龙格库塔方法求解常微分方程的数值解，供大家交流。
; [3 y6 G) {- f0 x; d( u8 E演示中，我使用了如下常微分方程作为测试：

) {3 A: b5 q/ C4 F; V8 b) A这个方程的解析通解是：

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

1. '用龙格库塔法求解常微分方程dydx=-y*cos(x)+exp(-sin(x)),y(0)=0在区间[0,10*pi]的数值解% r; G3 F: O/ K8 O
   
2. '已知这个微分方程的解析解 y=exp(-sin(x))*x2 `  X9 }6 [& s0 V
   
3. ( u# f$ ~% N: [. V5 Q) n
   
4. '生成一个workfile作为基本的数据容器0 J  f( [( ~; |! ~1 K
   
5. wfcreate (wf=temp) u 1000
   , U( C4 w! \3 ^' [( ?
6. ( b$ \* W7 X. r2 v7 L
   
7. '定义常量
   1 e8 E6 M+ h+ J; F% f
8. scalar pi=3.14159
   6 `/ v' s' ?( `! }% g+ q
9. scalar a=0        '定义自变量下限
   ( K+ J& V/ l+ @( h
10. scalar b=10*pi     '定义自变量上限
    9 y' c9 j/ f; M
11. scalar  M=500       '定义步数
    7 {$ Z) X) S5 B
12. scalar h=(b-a)/M   '计算每步之间的间隔
    , y! \8 z. k+ h
13. 
    & @0 }3 x9 `/ ^/ m$ w8 w+ X
14. '定义一个矩阵来储存计算数据，其中第一列储存自变量数据，第二列储存因变量数据,第三列储存解析解的值用以作为比较
    8 P  W1 S2 v2 i1 T; E
15. matrix(M+1,3) F
    ( D2 u3 V8 k1 Z7 u! T
16.   W3 {3 C# x6 }& o8 [5 Y
    
17. '矩阵的第一行储存初值问题的初始条件$ H% O: z2 |* e) }5 x
    
18. F(1,1)=0
    # c7 j6 `2 X* R3 n% B, ?
19. F(1,2)=0
    7 s6 J4 u/ T6 i
20. F(1,3)=@exp(-sin(F(1,1)))*F(1,1)9 @  J' z. M+ x$ L& Z6 _3 `0 Y7 \% |
    
21. 
    # Q- \& ]+ }8 L! O
22. '定义龙格库塔法的权重参数
    : P! J) a6 d2 @! v\" i7 g\" N0 B0 _6 A# U
23. scalar k1
    ( G: w$ C# F. x\" I/ h
24. scalar k2
    + S- y* d/ k7 @9 Z
25. scalar k34 y8 K8 R; [; @. ~8 T
    
26. scalar k4
    0 S% E( U7 X& J5 C4 R
27. 
    # x  \, Y1 l, E1 K. H
28. '定义权重的过程量
    7 |% x3 q3 ^( S+ D0 |\" @2 h3 V
29. scalar w1
    ( n# }  h\" s8 o, D: u9 s7 n\" @) y2 Z
30. scalar w26 _\" O( f& \- |5 I
    
31. scalar w3
    ) b/ G2 U: y- O8 f$ ?
32. scalar w4/ N/ W. q! F4 q
    
33. . ?  o1 d3 w' f. d2 ]
    
34. '程序主体% ]) }- S5 Y  X1 j4 T, ]% Z
    
35. for !k=1 to M step 1
    ; k; \. h3 L: ^$ m
36.   F(!k+1,1)=F(!k,1)+h
      u* A# |5 A, M. K/ x' A
37.   '调用常微分方程计算权重! p9 X$ ~1 D\" C  ^
    
38.   call obj(w1,F(!k,1),F(!k,2)) - H- O9 X' ~2 V7 x  U
    
39.     k1=w1*h5 i' ?- b% Z9 w$ ?2 p  N
    
40.   call obj(w2,F(!k,1)+h/2,F(!k,2)+k1/2)
    8 C: @3 A, k% L9 x# W+ T
41.     k2=w2*h
    6 G2 B) Y, b- H! f' S
42.   call obj(w3,F(!k,1)+h/2,F(!k,2)+k2/2)
    : e* |6 {2 t2 C  b
43.     k3=w3*h- N$ g1 _. @% G! y
    
44.   call obj(w4,F(!k,1)+h,F(!k,2)+k3)8 u5 R3 ^9 @, v( w0 i6 f; y' q% j
    
45.     k4=w4*h
    6 o% L& v3 b2 d7 ^. R4 G
46.   '计算函数估计值; m: z( m2 E% S, Z1 U$ \7 R
    
47.   F(!k+1,2)=F(!k,2)+(k1+2*k2+2*k3+k4)/6
    9 m* y$ G9 a( _9 C9 b( B
48.   '计算函数解析值' s7 o9 `8 t1 C! K
    
49.   F(!k+1,3)=@exp(-sin(F(!k+1,1)))*F(!k+1,1)
    & S% K3 `% U; L. b$ [) j
50. next4 \- I2 U7 X\" a* {8 @7 \
    
51. 3 R; j8 H$ b! l) H2 B
    
52. '显示最终结果0 y! J/ }! C# n0 @. O
    
53. freeze F.xyline
    $ s6 V, K) {7 @
54. freeze F) Y: V+ r& p9 X+ p
    
55. 
    % W( Y9 W& R# ^, F& n
56. '定义常微分方程
    $ {5 w' C( l3 c\" u  }4 X
57. subroutine obj(scalar dydx,scalar x,scalar y)
    ; D5 T* P7 f! ^5 g' R: \/ |
58.   dydx=-y*@cos(x)+@exp(-@sin(x))$ ?\" x' c5 n0 [/ `
    
59. endsub
    8 L' }0 F( E1 Y0 x* w# k\" o

复制代码

运行后求得结果如下：

$ F' k1 x1 q0 a+ h* ^
其中C2列是数值解，C3列是解析解，比较之下，这二者之间无明显差异。
: X, ~& Z( H: ?" F/ ^* H7 g
9 j& T# y! B; y- U$ h0 k4 b1 k  i