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

liwenhui

) s& s3 R- m' U7 O3 C
+ k  }- z1 z* K7 s: A0 mEViews除了能解决计量经济学的估计问题以外，还提供一个编程环境用以解决复杂的问题。经过调试，我在EViews中实现了用龙格库塔方法求解常微分方程的数值解，供大家交流。
演示中，我使用了如下常微分方程作为测试：
. v- a3 U9 h" x/ Q3 E7 E3 F

这个方程的解析通解是：

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

1. '用龙格库塔法求解常微分方程dydx=-y*cos(x)+exp(-sin(x)),y(0)=0在区间[0,10*pi]的数值解0 ?8 d% _. `% ~6 U1 t3 u2 V; `
   
2. '已知这个微分方程的解析解 y=exp(-sin(x))*x\" v6 ^& N, u6 u$ _: {( h; R! \
   
3. / Q  L$ V  ^9 L! ?\" ?  E
   
4. '生成一个workfile作为基本的数据容器& @. a! Q8 y* `' V- G# U( `
   
5. wfcreate (wf=temp) u 1000
   . O' ]6 w+ Z# Y- [: g6 z. f
6. & ~* E' a, Z. x3 R% u4 i
   
7. '定义常量
   ( A! z9 U9 C+ C+ ^
8. scalar pi=3.14159
   ) G* ~0 m* A- A0 O  C, T
9. scalar a=0        '定义自变量下限0 f! T! b/ S# j) j* z( c) y- w5 U
   
10. scalar b=10*pi     '定义自变量上限
    : K; T% a+ Q8 w/ a\" x1 g
11. scalar  M=500       '定义步数
    2 {9 S+ Z( \' `
12. scalar h=(b-a)/M   '计算每步之间的间隔& H+ d1 C\" {! I; B% _
    
13. 
    6 s# y& Z# R0 u! J* t. ^
14. '定义一个矩阵来储存计算数据，其中第一列储存自变量数据，第二列储存因变量数据,第三列储存解析解的值用以作为比较+ `5 ]! g$ W2 e. [, j
    
15. matrix(M+1,3) F 4 u, k! b: ?+ ~! f/ d
    
16. 5 |( {  M' Z\" n9 s2 ]
    
17. '矩阵的第一行储存初值问题的初始条件
    ! {6 ^( }\" c! v. V
18. F(1,1)=0
    ) G) E' S9 l6 [8 u. Y9 E# M- W
19. F(1,2)=0
    ' p4 A% Z4 O8 ]' D4 X
20. F(1,3)=@exp(-sin(F(1,1)))*F(1,1)7 Z+ ?( n# r4 k
    
21. $ k8 G; j& E; X* t
    
22. '定义龙格库塔法的权重参数
    8 [* F2 w0 O, {
23. scalar k1
    7 a: w- x2 z: }1 [
24. scalar k2% P6 l3 a( V4 k: g2 z# Q
    
25. scalar k3\" m, @( s5 H5 Q0 |# G
    
26. scalar k42 A1 _2 y* K% B6 W
    
27. ) _. ?9 r$ I5 D: [& v! d8 e. A
    
28. '定义权重的过程量
    9 [( w- P\" P5 Y9 L; A) H
29. scalar w1
    , t\" P0 P1 g5 m4 j$ y0 f0 T
30. scalar w2/ X$ @( ?8 c# |1 _* K' K0 G& \6 F
    
31. scalar w3
    7 y# o5 W- B8 w# T1 F( W: y- d
32. scalar w4
    . P: Z: P. ~1 o: r( Q
33. * U: E+ f% X$ T2 c5 v  l# L
    
34. '程序主体
    4 W0 x+ x' P0 E- X0 y! H
35. for !k=1 to M step 1
    & E! X, w; I8 [  H* G- w
36.   F(!k+1,1)=F(!k,1)+h) H5 D4 P/ f) l  B( S3 w3 i# W
    
37.   '调用常微分方程计算权重
    ) E& r; Z4 G0 E/ I( M, }& U* x
38.   call obj(w1,F(!k,1),F(!k,2)) 9 r+ u% G  @8 \\" m. a; T
    
39.     k1=w1*h; E  Z' G$ n. @
    
40.   call obj(w2,F(!k,1)+h/2,F(!k,2)+k1/2)
    2 K8 x4 w0 U* F) V+ V0 R
41.     k2=w2*h
    & d) Q4 W. |( k( |; G
42.   call obj(w3,F(!k,1)+h/2,F(!k,2)+k2/2)
    , f9 R4 {% ^0 l* {( @
43.     k3=w3*h1 v1 o6 N, _* \5 n4 ?+ W
    
44.   call obj(w4,F(!k,1)+h,F(!k,2)+k3)
    * q4 [1 W, `% n! ]\" W
45.     k4=w4*h
    & F: T% Q2 g! O
46.   '计算函数估计值
    , G' X( k9 v; `$ g4 f; U\" C1 H1 J. `
47.   F(!k+1,2)=F(!k,2)+(k1+2*k2+2*k3+k4)/6. p# ~3 E' L# p% a* j3 [1 L! I
    
48.   '计算函数解析值. O1 O\" X6 o- v. z) I
    
49.   F(!k+1,3)=@exp(-sin(F(!k+1,1)))*F(!k+1,1)
    6 d1 r5 t! O) s: o3 C
50. next\" e% c3 W; v6 t; ?1 H  W8 D, `
    
51. 
    3 T\" |7 J$ V  E
52. '显示最终结果% o' L! g5 \! X* C  l+ Z
    
53. freeze F.xyline
    0 e- u4 c6 K2 v$ ?- [/ N
54. freeze F
    2 M& |: e7 s6 |' _
55. 
    ( |9 X& N  u1 g5 C7 d+ [8 ^
56. '定义常微分方程
      V! J6 ^: g$ N
57. subroutine obj(scalar dydx,scalar x,scalar y)
    8 D7 k/ Z5 C' j7 t& Y$ \9 J2 z# X
58.   dydx=-y*@cos(x)+@exp(-@sin(x))* |: G- o& ]6 U( l
    
59. endsub
    . v! H, h- t# w# j! F7 T

复制代码

运行后求得结果如下：
2 W" n" Y9 g5 ?) |' `1 R( R

/ U  d4 D' @" k! Z/ k3 a9 H
# x4 s8 p0 ~1 C3 O4 h其中C2列是数值解，C3列是解析解，比较之下，这二者之间无明显差异。

) S( x( J' S+ r1 F