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

liwenhui

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

这个方程的解析通解是：
( s6 B- h  |0 p( \, Z8 f

* ^& ~& i8 L  q9 f( G4 O3 Q7 C
使用“龙格库塔方法”，编制的EViews程序如下：

1. '用龙格库塔法求解常微分方程dydx=-y*cos(x)+exp(-sin(x)),y(0)=0在区间[0,10*pi]的数值解
   & s& B; J9 k& j) e\" p8 [: M
2. '已知这个微分方程的解析解 y=exp(-sin(x))*x
   \" Y0 e1 T' I4 m& r% H
3. ! x3 g! v\" _9 j/ i, C1 F3 a# L
   
4. '生成一个workfile作为基本的数据容器9 t\" C& H* I+ g3 A0 Y' F
   
5. wfcreate (wf=temp) u 1000
   4 l! h% I3 v5 o\" m
6. 
   . V\" n0 {# ^$ b0 j3 m& g\" A
7. '定义常量/ g( u6 L: b; D9 i  p
   
8. scalar pi=3.14159! Z, h3 H4 y. s
   
9. scalar a=0        '定义自变量下限
   ) V* \' {4 N; t3 ~
10. scalar b=10*pi     '定义自变量上限: x8 _5 K  M- [\" x6 i
    
11. scalar  M=500       '定义步数( u# O) W# O7 Y3 V/ P
    
12. scalar h=(b-a)/M   '计算每步之间的间隔
    ' i! N; U% L8 x# C. `7 X
13. 5 l/ T9 t; T  o0 r, b5 d
    
14. '定义一个矩阵来储存计算数据，其中第一列储存自变量数据，第二列储存因变量数据,第三列储存解析解的值用以作为比较
    1 D* x6 M7 ?7 E0 Q  V\" O' F
15. matrix(M+1,3) F / c; M) [7 \% [9 i2 C
    
16.   G4 [. j6 N+ O* x4 z) R
    
17. '矩阵的第一行储存初值问题的初始条件
    5 h: ~7 V. a2 c0 S
18. F(1,1)=0
    % f) H% n( D5 d1 z  f
19. F(1,2)=0( \# v+ Q. S: i\" Q
    
20. F(1,3)=@exp(-sin(F(1,1)))*F(1,1)
    ; X4 i4 g- N0 L
21. - ~$ j8 P& l$ ]
    
22. '定义龙格库塔法的权重参数
    + L& K& F# W: ?: l5 C) p
23. scalar k1+ a; _/ Y4 v2 B4 i
    
24. scalar k2
    / v8 o9 q6 t\" y# u. x2 i
25. scalar k3
    $ ^6 s* U  ~$ V( y9 V( ?5 X) X# Q
26. scalar k4  y, c# M6 ]$ i& V& Y
    
27. % |% C& A( `- [  O
    
28. '定义权重的过程量* N9 \; U- _& j, t\" p, J
    
29. scalar w1& f9 P# i2 N/ b
    
30. scalar w25 T+ j, K; A$ U8 ^
    
31. scalar w3
    * \) c\" V, r5 U\" t
32. scalar w4% u, o. e, [8 p  A- o, V8 Z
    
33. 
    * n& i9 E/ T) `- _
34. '程序主体, o! m8 s4 F8 x3 y, {
    
35. for !k=1 to M step 1! M) @3 ~4 J% Z( |' h
    
36.   F(!k+1,1)=F(!k,1)+h1 O2 T\" q! I; X3 m
    
37.   '调用常微分方程计算权重
    . T; T& L4 M$ h) _  O
38.   call obj(w1,F(!k,1),F(!k,2)) 7 n( e% G( |5 S% v# y\" T5 I6 \
    
39.     k1=w1*h: R+ S* i, V6 I8 c2 e( I
    
40.   call obj(w2,F(!k,1)+h/2,F(!k,2)+k1/2); C' {  X( A4 y0 V+ E2 ?
    
41.     k2=w2*h
    4 C# ~# g9 Q2 W/ H! ^) S
42.   call obj(w3,F(!k,1)+h/2,F(!k,2)+k2/2)
    7 \9 k* G9 b! x6 p% V2 J6 k2 V
43.     k3=w3*h' i: B% l8 X- U) z
    
44.   call obj(w4,F(!k,1)+h,F(!k,2)+k3)6 C& U& @/ Y3 @1 e6 Y
    
45.     k4=w4*h
    # I, S7 D! J\" t2 G. J' _\" Z
46.   '计算函数估计值
    & j7 g2 j8 a' w4 i3 X! ?# v8 t
47.   F(!k+1,2)=F(!k,2)+(k1+2*k2+2*k3+k4)/6
    / T3 a, p: g# a\" F3 X0 H
48.   '计算函数解析值
    6 c2 |  p1 B0 G
49.   F(!k+1,3)=@exp(-sin(F(!k+1,1)))*F(!k+1,1)5 Z' l  _1 |! b7 F. Z
    
50. next; ^( t) F6 J; ~
    
51. 
    & a\" O; u( g  S# p9 l
52. '显示最终结果
    1 N7 E4 d* a8 ~
53. freeze F.xyline
    4 c/ s/ U% R3 w5 x3 j6 X6 a
54. freeze F
    ( a: ]5 ~- W! [' I' P0 r/ G4 T( U# Z& ~
55. 
    5 Z2 d, h: y. T; _- I. e1 |$ y8 m
56. '定义常微分方程+ R* f4 G6 |. v! E/ @
    
57. subroutine obj(scalar dydx,scalar x,scalar y)0 s+ ]3 G8 [9 w6 T/ Z6 u1 I
    
58.   dydx=-y*@cos(x)+@exp(-@sin(x))4 H. J! C6 U6 S- a4 j) e
    
59. endsub
    8 R3 w7 L. i, L7 U. L

复制代码

运行后求得结果如下：
! x) {  H, W; u
* a2 `7 Y9 D7 p6 }2 [

& C3 U4 X1 \) M
4 u4 U% a. h6 [  x, O8 O5 p* Q- r其中C2列是数值解，C3列是解析解，比较之下，这二者之间无明显差异。
: g4 K- \% Z5 u1 d
1 n& [- R0 b# I8 r4 a, Z+ `& e

) T! W& i( u+ `+ e