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

liwenhui

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

! w% m% H+ X+ U( k* I3 c+ ^这个方程的解析通解是：

, @4 z5 O7 ^2 _" F' t0 f$ W* D; E使用“龙格库塔方法”，编制的EViews程序如下：

1. '用龙格库塔法求解常微分方程dydx=-y*cos(x)+exp(-sin(x)),y(0)=0在区间[0,10*pi]的数值解# R0 E! s  Y7 H$ ?
   
2. '已知这个微分方程的解析解 y=exp(-sin(x))*x# W* e1 `* e. `0 T2 Z6 g( B1 g
   
3. 
   : |/ b: N: H8 X- e. r9 C% `
4. '生成一个workfile作为基本的数据容器% d3 f( ^; g1 l$ i
   
5. wfcreate (wf=temp) u 1000/ [/ A) Z+ l4 z1 c! z& v7 E
   
6. 7 _! V2 q+ m: z$ m# y5 X! Q
   
7. '定义常量
   - g3 J# x9 G9 P# O* ?) Q
8. scalar pi=3.14159
   : N8 M( \) [3 c: S: k! C: N
9. scalar a=0        '定义自变量下限
   ! i\" P# U5 T% K4 ?
10. scalar b=10*pi     '定义自变量上限
    6 k$ J) d, c- Y: V8 n
11. scalar  M=500       '定义步数) a$ h% y# v( g4 U, n' J
    
12. scalar h=(b-a)/M   '计算每步之间的间隔
    3 r6 w' c$ m& g, D
13. % D3 z8 V, V% F$ L
    
14. '定义一个矩阵来储存计算数据，其中第一列储存自变量数据，第二列储存因变量数据,第三列储存解析解的值用以作为比较8 u, A0 l% v) l7 ]% b# @6 V2 l
    
15. matrix(M+1,3) F ) t% t+ a5 m1 O$ x0 Q  @
    
16. 
    0 b% i! N3 l+ x! ]
17. '矩阵的第一行储存初值问题的初始条件: P7 w5 ~. s: E5 d) s
    
18. F(1,1)=06 ^6 s/ z/ [5 i7 j
    
19. F(1,2)=01 J2 l5 I\" Q) B  R+ c\" ]\" y: {
    
20. F(1,3)=@exp(-sin(F(1,1)))*F(1,1)$ Q/ ^& ^* M2 w
    
21. 
    / p' c' E  }) J5 S7 s/ Q: i
22. '定义龙格库塔法的权重参数5 H( u# `. w2 c( X' l
    
23. scalar k1
    4 X! Y7 R. q7 X8 L2 `* i
24. scalar k2
    4 ]2 ~3 ~! l5 u5 X/ g5 G
25. scalar k3
    + t8 l- T: d8 P$ s3 q% {
26. scalar k4
    ( r0 ~2 z* s: `
27. 
    8 A* y- B8 A+ L& u  I2 ~) F
28. '定义权重的过程量
    9 P; d# _5 I. T' j& R7 P
29. scalar w1
    + X: [& Q0 q4 a$ s
30. scalar w2
    + g3 m, R1 F) A$ o
31. scalar w3+ {) u, Q0 t0 n$ q2 F) q% x, G
    
32. scalar w4( X$ c- @8 m+ t3 T# g9 C% J( X8 e
    
33. 
    1 a9 o+ G( M9 O, ~# ?- s
34. '程序主体5 p5 x: e( K% ^8 G4 _
    
35. for !k=1 to M step 1% a1 P: Q! j: }& {( S4 o
    
36.   F(!k+1,1)=F(!k,1)+h
    : i* m4 }$ {( M7 [2 K) C
37.   '调用常微分方程计算权重
    : E. |0 c% Z0 f) m$ U6 P( t0 J
38.   call obj(w1,F(!k,1),F(!k,2))
    $ U$ X3 `  Q8 ?! d8 K3 v& Q. D
39.     k1=w1*h
    , Q- K\" r2 ~0 Y
40.   call obj(w2,F(!k,1)+h/2,F(!k,2)+k1/2)$ N6 O\" k2 x$ `6 v7 e$ j; d
    
41.     k2=w2*h
    + f8 D& l0 Y: a% Y, ^; _# G0 T
42.   call obj(w3,F(!k,1)+h/2,F(!k,2)+k2/2)
    9 U: {' D$ U% I: Q' T8 a
43.     k3=w3*h
    ( `8 t2 c% V1 {9 i
44.   call obj(w4,F(!k,1)+h,F(!k,2)+k3)
    : n3 e, [1 m$ G# {
45.     k4=w4*h- i) g; E3 {; M, q$ g2 Z. f
    
46.   '计算函数估计值
    5 M( J3 w4 l0 d9 B
47.   F(!k+1,2)=F(!k,2)+(k1+2*k2+2*k3+k4)/6\" G8 \* F& b2 v+ g$ K( V+ t( ~
    
48.   '计算函数解析值8 c, j) C$ e2 v2 I6 t
    
49.   F(!k+1,3)=@exp(-sin(F(!k+1,1)))*F(!k+1,1)
    + ~* U; R; c3 q' h, r1 @% U
50. next
    1 c' R1 O% u( `8 O$ e
51. 
    % W- }9 w, B6 a1 z
52. '显示最终结果
    9 Q\" j6 p6 ~, @7 L
53. freeze F.xyline
    - T% I; M  x* [/ n0 O\" p
54. freeze F4 J. y6 c& R# Q4 P+ ?
    
55. 2 c% x\" M\" o3 I( I* s
    
56. '定义常微分方程
    $ |3 F8 T/ p- f  S9 [
57. subroutine obj(scalar dydx,scalar x,scalar y)
    1 S' _! `/ a# r6 z
58.   dydx=-y*@cos(x)+@exp(-@sin(x))0 e% U& d# j8 f$ r: ^0 W8 _6 `
    
59. endsub
    + y4 w( B3 G% j2 t) N

复制代码

运行后求得结果如下：
& s2 s. h. L' {

' l; w/ o; O  O: @: Q
其中C2列是数值解，C3列是解析解，比较之下，这二者之间无明显差异。

7 k6 o9 \5 w( d* b# M  e# _/ I' m
. o* |+ H& [1 l! i! T) K7 W0 l
. V2 c5 n8 j$ `( M2 B