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

liwenhui

EViews除了能解决计量经济学的估计问题以外，还提供一个编程环境用以解决复杂的问题。经过调试，我在EViews中实现了用龙格库塔方法求解常微分方程的数值解，供大家交流。
4 l: u% n0 ~. w% ^' g- S) @演示中，我使用了如下常微分方程作为测试：

这个方程的解析通解是：
& p) L0 ~# T& k3 D1 R

, p2 ^8 m! ?" U, A: I: i: P, s使用“龙格库塔方法”，编制的EViews程序如下：

1. '用龙格库塔法求解常微分方程dydx=-y*cos(x)+exp(-sin(x)),y(0)=0在区间[0,10*pi]的数值解
   # X: d, |9 d7 x1 U9 H
2. '已知这个微分方程的解析解 y=exp(-sin(x))*x
   : H4 L' v$ }6 ~, n: M$ H
3. 
   4 |9 X1 e4 W% a
4. '生成一个workfile作为基本的数据容器
   + G$ p7 a# L7 |  C+ T1 s6 A) v
5. wfcreate (wf=temp) u 1000
   ! L0 b2 k& p4 p& t
6. 
   7 {2 \; Z, F: L
7. '定义常量
   7 f9 v' [( N6 ?- b0 C  w5 l
8. scalar pi=3.141599 m4 s* L! A. R
   
9. scalar a=0        '定义自变量下限9 y' N3 m9 }1 d* f$ l3 m
   
10. scalar b=10*pi     '定义自变量上限3 w% |% c  f/ k# c
    
11. scalar  M=500       '定义步数
    . g& p* y- c9 |4 G( N7 g1 o
12. scalar h=(b-a)/M   '计算每步之间的间隔
    % e; E7 T( z% z  {/ I, t
13. 
    0 E) g8 U) W; ~! a1 s
14. '定义一个矩阵来储存计算数据，其中第一列储存自变量数据，第二列储存因变量数据,第三列储存解析解的值用以作为比较/ u. H- U/ g: r9 n) _$ f
    
15. matrix(M+1,3) F 0 l$ L8 @5 h# K& @6 {: H& z
    
16. 4 t5 M& ^* c. P# S3 C6 i
    
17. '矩阵的第一行储存初值问题的初始条件7 v\" d8 x+ r; R/ s2 N1 r( U4 U
    
18. F(1,1)=0' d\" P1 I. w+ A6 f9 b
    
19. F(1,2)=0
    % y5 E0 o4 F3 V4 U
20. F(1,3)=@exp(-sin(F(1,1)))*F(1,1)
    3 Q* ^5 w; x6 U) r
21. 
    : r/ @# @: a  z: |5 T+ F7 q
22. '定义龙格库塔法的权重参数
    + O/ w4 E% [: _
23. scalar k1) N% A6 p7 q, k4 j9 D) P
    
24. scalar k2
    8 R! k$ ?+ N% _\" U5 z. V5 J3 q
25. scalar k3
    \" _# ^$ L7 e* {) G
26. scalar k40 ?9 `( f' i% ?3 e4 f
    
27. 9 P0 z  M; W+ j) i7 v- c1 o
    
28. '定义权重的过程量; Z4 f$ `/ q$ o$ {4 O; |; e
    
29. scalar w1# B! E: U) A7 _9 v
    
30. scalar w2  H* d% Z9 Z8 `) \
    
31. scalar w3
    - e+ h* \7 O' w7 c* {
32. scalar w4
    * Q: `! t+ |- T. W# _
33. 3 @: Y( C/ C\" j: f
    
34. '程序主体# b5 t2 h% f& w4 U! P; b- P
    
35. for !k=1 to M step 1  N( ]8 t! L  P+ j6 M
    
36.   F(!k+1,1)=F(!k,1)+h7 n( x$ B1 m5 {) ~4 j8 x% b
    
37.   '调用常微分方程计算权重
    5 f5 ?8 ]( I; H& G& ^3 }$ G
38.   call obj(w1,F(!k,1),F(!k,2))
    . C4 D3 Z4 ]4 J( \1 w# x
39.     k1=w1*h
    $ e# q# Z8 O9 ]- Q3 U$ K5 u
40.   call obj(w2,F(!k,1)+h/2,F(!k,2)+k1/2)
    $ A4 @1 z& T4 z/ g; T$ g
41.     k2=w2*h0 l1 x1 e4 R6 [0 `- A
    
42.   call obj(w3,F(!k,1)+h/2,F(!k,2)+k2/2)5 k' B( I1 I( x& c1 L% s3 k- B
    
43.     k3=w3*h# u% [7 r6 B, [( e% {2 Q% x
    
44.   call obj(w4,F(!k,1)+h,F(!k,2)+k3); b+ g) [) i/ P* Q7 E
    
45.     k4=w4*h
    4 |9 x: c) d. f\" C* Y\" P! w3 p
46.   '计算函数估计值6 J  A1 i3 C/ J\" [/ |# a8 g, a
    
47.   F(!k+1,2)=F(!k,2)+(k1+2*k2+2*k3+k4)/6
    \" |4 Q/ f' {' r/ [
48.   '计算函数解析值1 @0 Q; @* ]# o. d/ |. ^$ M
    
49.   F(!k+1,3)=@exp(-sin(F(!k+1,1)))*F(!k+1,1)
    , H  Q& v% R- {- H
50. next- K. C, I  T6 P, i$ h* [
    
51. ' v% q, a\" ~% i3 G& k
    
52. '显示最终结果) b! {, `) M$ x, b* V9 F5 i
    
53. freeze F.xyline' r0 M) t/ f5 w0 R+ W+ m1 e+ x
    
54. freeze F, p% [; a3 Z( R* r
    
55. 
    \" ?. s4 Q, q$ X7 l
56. '定义常微分方程
    : }( ^* X, [6 h2 p2 E3 H
57. subroutine obj(scalar dydx,scalar x,scalar y)0 R. ^8 C. B  m/ \
    
58.   dydx=-y*@cos(x)+@exp(-@sin(x))
    # s0 a' ~8 ?- A
59. endsub
    , g  C# I5 H4 H, S7 k# X% Z: \

复制代码

运行后求得结果如下：
2 f  |& k2 j) Q2 f
# }% o9 p) |  h  Y6 N

其中C2列是数值解，C3列是解析解，比较之下，这二者之间无明显差异。

/ G8 k" w+ Q. D) z$ I$ x8 f
- f1 O( H! h6 A. l$ Z* s

, L5 H: M3 H+ F