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

liwenhui

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

( ?7 ?& X7 [. g这个方程的解析通解是：
7 f! Y# c$ Q3 g" {+ d

" H" D' y! G1 r, u+ [. |
2 c# w, k1 Z2 e/ X0 r: V使用“龙格库塔方法”，编制的EViews程序如下：

1. '用龙格库塔法求解常微分方程dydx=-y*cos(x)+exp(-sin(x)),y(0)=0在区间[0,10*pi]的数值解
   $ @* ^; c  X; V( J; p3 m
2. '已知这个微分方程的解析解 y=exp(-sin(x))*x
   ; R+ }* \( ^! M- P  n, E
3. 
   * F6 @5 @( ?+ t, Z! l
4. '生成一个workfile作为基本的数据容器
   ! Y7 D) n5 K- ^' ^( O\" Q
5. wfcreate (wf=temp) u 10006 U; F  F; a; K: H
   
6. 
   + Z; h$ J/ ?$ {- B* x
7. '定义常量( f1 {; G) ^\" n- g& J4 ]7 {! s
   
8. scalar pi=3.14159& R0 U' F5 J2 q1 f3 Z. A
   
9. scalar a=0        '定义自变量下限/ s& ^8 m  R( D  s
   
10. scalar b=10*pi     '定义自变量上限2 g1 f; I  [\" G2 Z( R
    
11. scalar  M=500       '定义步数' ^, u; M3 |6 Z\" l/ |( g/ E7 T
    
12. scalar h=(b-a)/M   '计算每步之间的间隔
    1 K9 B\" E  Y! i: H2 W/ f0 ^* V+ q
13. 
    - c/ F! G: ]- I; ?9 S
14. '定义一个矩阵来储存计算数据，其中第一列储存自变量数据，第二列储存因变量数据,第三列储存解析解的值用以作为比较
    \" Q- s  m3 E# \
15. matrix(M+1,3) F ; F4 J: ^6 u3 m4 @+ f
    
16. ' X* t% w+ }: D; I' i% H: R
    
17. '矩阵的第一行储存初值问题的初始条件
    ; z\" D; z2 M1 O& t
18. F(1,1)=03 @; |( H+ h+ w# W- _3 @6 {
    
19. F(1,2)=0
    9 f8 c! u: G* J
20. F(1,3)=@exp(-sin(F(1,1)))*F(1,1)
    4 u# a# L+ u4 ^+ o
21. 
    $ w# E: Q; K3 D9 B& c- @  u! C
22. '定义龙格库塔法的权重参数
    9 B# H$ m3 t% A
23. scalar k1
    ) n3 Y) p) D: J% b5 k9 _  L  E# B
24. scalar k2# R. D& m+ K; A/ c  S$ ~9 O7 Y
    
25. scalar k3: |\" S/ j9 ?& \' x. n
    
26. scalar k4
    ( H' j& r1 {) ^, J' S9 w
27. - {+ n9 ^, Z7 z8 I- D7 W3 J
    
28. '定义权重的过程量
    ; }- C1 j& p6 p3 U4 Q9 s
29. scalar w1
    8 a3 e  ?3 _) B, k5 |- J4 \7 K
30. scalar w2
    - o2 Q5 Z\" b) p3 ~5 O# `, R$ `
31. scalar w30 G& C& N; g7 U( p, K3 k
    
32. scalar w4/ g, `8 X- f# E: [6 P) c\" g
    
33. 8 V! ]( ^\" v6 ^6 G5 ]. c
    
34. '程序主体
    4 s) x( ]+ G# z  S! t
35. for !k=1 to M step 1
    ; c% {( `2 @7 W
36.   F(!k+1,1)=F(!k,1)+h
    ( ~3 ^- R- E( i& S
37.   '调用常微分方程计算权重4 C  F/ D: E0 M5 C4 z+ J1 b! W
    
38.   call obj(w1,F(!k,1),F(!k,2)) 8 Q! E/ a  |7 I) d( \$ @- E* p\" P
    
39.     k1=w1*h) ~6 C# N# l0 C! N
    
40.   call obj(w2,F(!k,1)+h/2,F(!k,2)+k1/2)( j- T) @. ]0 K: j
    
41.     k2=w2*h, \; r$ x( R# P# K8 O& O3 B
    
42.   call obj(w3,F(!k,1)+h/2,F(!k,2)+k2/2)
    + r( Z. Y- p' E5 S2 n5 s
43.     k3=w3*h
    * s4 ]5 ?8 L, W) c: Q1 L
44.   call obj(w4,F(!k,1)+h,F(!k,2)+k3)
    ; @5 C4 A2 A4 m) F% |0 {) P: Y
45.     k4=w4*h+ P3 p& j# f0 I5 o- [3 I
    
46.   '计算函数估计值  l4 o  v6 ~# L7 {
    
47.   F(!k+1,2)=F(!k,2)+(k1+2*k2+2*k3+k4)/64 A9 X9 D  W, ^5 Z% L& W% T
    
48.   '计算函数解析值
    8 w6 y/ X& `, @5 a! `
49.   F(!k+1,3)=@exp(-sin(F(!k+1,1)))*F(!k+1,1)
    + w, F/ X8 ^1 b3 y, V/ Q' O1 @3 Z
50. next4 K( t+ L0 W. a, C% x
    
51. ( k4 d2 d; m9 a3 Q/ \
    
52. '显示最终结果. v( J/ q/ y8 F: r, j
    
53. freeze F.xyline
    : R\" L8 W& p4 `5 ?
54. freeze F& s$ G/ {) {* t
    
55. 7 [6 \* H, G2 I7 r6 O\" o
    
56. '定义常微分方程) e* r0 U8 W\" x
    
57. subroutine obj(scalar dydx,scalar x,scalar y)* v8 l7 D. r0 }( _2 |: T
    
58.   dydx=-y*@cos(x)+@exp(-@sin(x))
    ) L2 J/ @, W* _
59. endsub& ~' e- G- S' V2 W8 k
    

复制代码

运行后求得结果如下：

3 k; ~! K% H+ r% T
其中C2列是数值解，C3列是解析解，比较之下，这二者之间无明显差异。

8 w4 ^% K/ C/ l# E; h: O

, u- P$ O5 D) Q$ O  X1 S1 ~