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

liwenhui

EViews除了能解决计量经济学的估计问题以外，还提供一个编程环境用以解决复杂的问题。经过调试，我在EViews中实现了用龙格库塔方法求解常微分方程的数值解，供大家交流。
3 {* `) Y9 l, P+ ^演示中，我使用了如下常微分方程作为测试：
" u; v/ V# l1 r. F& K

# G8 @* u% t. B4 T& }* n9 E& N
6 z2 q: K" L4 s% C9 w这个方程的解析通解是：
0 S% L) J1 j6 h! k/ M# Z- x1 [. g7 \

5 y, B+ ?& \! t使用“龙格库塔方法”，编制的EViews程序如下：

1. '用龙格库塔法求解常微分方程dydx=-y*cos(x)+exp(-sin(x)),y(0)=0在区间[0,10*pi]的数值解
   & s0 ?: J% b  b( {0 I. Q
2. '已知这个微分方程的解析解 y=exp(-sin(x))*x
   # R6 B6 P' I$ B\" Z, `\" ^, j3 ^
3. 
   6 V3 k) f0 J2 o  `5 Y: k2 }
4. '生成一个workfile作为基本的数据容器
   \" f% B' W- g* S) Y# {2 |  k% _
5. wfcreate (wf=temp) u 1000) Y1 @% l3 U, G! q3 F
   
6. 1 S3 n* _# J4 t' C) C
   
7. '定义常量
   \" h& k* Q& Q9 l4 t. i4 L, A
8. scalar pi=3.14159
   2 L\" ~5 Q& {* L( R6 G
9. scalar a=0        '定义自变量下限  s2 b! J# J\" G
   
10. scalar b=10*pi     '定义自变量上限- ?: k( b* M' N; v% ?3 V1 `: S
    
11. scalar  M=500       '定义步数& j$ D* n, U# b, P& H
    
12. scalar h=(b-a)/M   '计算每步之间的间隔% i, I+ k! l$ X' H8 ^4 v
    
13. 
    + \; V+ Q( M% A$ u
14. '定义一个矩阵来储存计算数据，其中第一列储存自变量数据，第二列储存因变量数据,第三列储存解析解的值用以作为比较\" s4 y\" [- K; V: {; k5 `6 T
    
15. matrix(M+1,3) F
    % ~; y0 G\" s7 ]3 e  W9 {; A
16. ( V3 ~6 M5 g9 O# C
    
17. '矩阵的第一行储存初值问题的初始条件
    ( y' U0 I/ N0 L\" [
18. F(1,1)=0. {5 ], j' f& u9 G
    
19. F(1,2)=0
    6 ]# x4 V5 K9 {, S8 ]
20. F(1,3)=@exp(-sin(F(1,1)))*F(1,1)
    2 k- _8 }+ [' M. X\" M
21. 
    , S  a$ q4 v+ d/ E' L5 N1 s\" i+ o  k
22. '定义龙格库塔法的权重参数
    3 Q' u3 ^0 N: J# E4 a6 c. W
23. scalar k1  D$ I- A  J' d( M2 k- m
    
24. scalar k2
    9 e# C9 b% j) {/ u8 H0 G
25. scalar k3) v. M! I( H' f. ]. ~
    
26. scalar k4- s  f& D4 ^, K+ z
    
27. 
    7 ]6 h1 C; ^- d/ x* s. }1 p
28. '定义权重的过程量4 G  r\" @0 l8 K
    
29. scalar w1
    ; x  \/ L8 ]+ v* I* @7 q% v6 k
30. scalar w2% t\" T, t2 C0 |9 J
    
31. scalar w3
    5 N5 o, e' J3 D0 s/ o9 P7 e
32. scalar w49 B+ w5 v9 U: J5 s6 L
    
33. 
    ; d1 X\" z: {/ K2 P! M$ V9 }- ^- N
34. '程序主体: N+ u# n! f1 c- n
    
35. for !k=1 to M step 1) K. y: H  @5 M
    
36.   F(!k+1,1)=F(!k,1)+h3 _1 _) A  m1 O
    
37.   '调用常微分方程计算权重
    , z( K\" M' e& R( m7 t. [' `
38.   call obj(w1,F(!k,1),F(!k,2))
    : i2 r( g/ L) z- \' f; r
39.     k1=w1*h
    9 K( g8 X8 b% u7 H0 P. r
40.   call obj(w2,F(!k,1)+h/2,F(!k,2)+k1/2)6 i) h. F8 X) ]6 F$ u
    
41.     k2=w2*h
      K3 Z7 b% s\" v. l\" L  x+ C/ O. u
42.   call obj(w3,F(!k,1)+h/2,F(!k,2)+k2/2)7 E* I& X% U* i, c# H* O* N\" L
    
43.     k3=w3*h3 C6 {5 V5 z0 T# q: E* I
    
44.   call obj(w4,F(!k,1)+h,F(!k,2)+k3)/ i* v' F4 F' o3 o, r+ d
    
45.     k4=w4*h
    3 M  ]# D) Z; t& o* r/ o
46.   '计算函数估计值
    3 ]9 C( Q/ Q; h6 ]* g( A
47.   F(!k+1,2)=F(!k,2)+(k1+2*k2+2*k3+k4)/6
    ' V& \+ b4 ^# q: u5 ~$ H( O
48.   '计算函数解析值) [8 p- `/ o' A* t5 {' H1 c
    
49.   F(!k+1,3)=@exp(-sin(F(!k+1,1)))*F(!k+1,1)- X; w) m4 _8 N9 K
    
50. next3 ]. @) k6 f! H5 O: N
    
51. ' }5 c$ l; t3 a7 N2 c7 j
    
52. '显示最终结果
    2 _  y# i8 @: ^& |' R$ [
53. freeze F.xyline
    8 q) O( w& I; r, d/ X# K6 Q
54. freeze F\" ?0 m. r9 P. ^, W0 @3 E2 ]# e1 A
    
55. 
    0 A; }\" N2 U, T% Y: Z
56. '定义常微分方程) q- v# G) b5 d0 Y- d0 O0 A( t
    
57. subroutine obj(scalar dydx,scalar x,scalar y)% a0 c3 X( x+ I
    
58.   dydx=-y*@cos(x)+@exp(-@sin(x))
    : [. u9 C: g7 ^/ B' [# A
59. endsub6 ]5 n0 X' O9 O5 ~7 p1 r
    

复制代码

运行后求得结果如下：
& M9 \7 m4 J/ }, N" o% J
+ a6 T  v# x4 d2 w1 p

+ }/ |8 M  ^! D" t; S/ R- f; B- z其中C2列是数值解，C3列是解析解，比较之下，这二者之间无明显差异。


  w: W- l% c: d( m9 F4 t$ B! u

  K; o# B  b- w  ^: }! L' j+ j