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

liwenhui

6 B  P' K! I/ n" D% E: {
EViews除了能解决计量经济学的估计问题以外，还提供一个编程环境用以解决复杂的问题。经过调试，我在EViews中实现了用龙格库塔方法求解常微分方程的数值解，供大家交流。
5 b3 L% e* e& P1 C  u5 c演示中，我使用了如下常微分方程作为测试：

9 F$ T- r6 A+ o' m( U! Z
这个方程的解析通解是：

% |7 x% G7 Q8 {/ H2 j
使用“龙格库塔方法”，编制的EViews程序如下：

1. '用龙格库塔法求解常微分方程dydx=-y*cos(x)+exp(-sin(x)),y(0)=0在区间[0,10*pi]的数值解7 K; n( B% @- K$ L0 }; G. A0 h
   
2. '已知这个微分方程的解析解 y=exp(-sin(x))*x
   ' c& }; ?7 T# m6 P; o
3. , Z9 ?$ F6 e% f! p3 V/ w
   
4. '生成一个workfile作为基本的数据容器! K! w1 T\" z) r- i2 b4 `\" ], P; u
   
5. wfcreate (wf=temp) u 1000
   # j; C' F1 x: ]4 a+ f
6. 
   \" G# `& Z8 r. s$ U/ C/ s8 W
7. '定义常量+ j\" B3 ?, p/ @) y; `
   
8. scalar pi=3.14159
   ! {/ S0 X9 k# Q, |* W2 V4 j7 X
9. scalar a=0        '定义自变量下限9 A* n( u- P0 ~- x- _
   
10. scalar b=10*pi     '定义自变量上限
    6 x: O4 ?1 G  }( @
11. scalar  M=500       '定义步数
    ! j- X; L5 D' a/ m9 A8 s
12. scalar h=(b-a)/M   '计算每步之间的间隔1 {# w( k0 i6 j
    
13. + ~! i  m, A& C6 o3 m9 a; ]$ B2 y: M
    
14. '定义一个矩阵来储存计算数据，其中第一列储存自变量数据，第二列储存因变量数据,第三列储存解析解的值用以作为比较/ x( W2 R6 e9 Z. ~8 a
    
15. matrix(M+1,3) F
    * y9 [) e9 ?8 X# m- @$ e- I
16. 
    8 h% g& C1 s1 g5 m
17. '矩阵的第一行储存初值问题的初始条件0 G# e0 J  }& r( l( F8 y0 F
    
18. F(1,1)=0* d% X. {/ M1 W9 y- K
    
19. F(1,2)=0  G5 m2 A. d3 h% a# M. Q
    
20. F(1,3)=@exp(-sin(F(1,1)))*F(1,1)9 {6 [4 A% t* q' G
    
21. 
    # e7 F\" g* K  C- H
22. '定义龙格库塔法的权重参数
    - ~1 g: j4 A* [
23. scalar k1
    2 n$ @, ~6 K5 \) K$ I/ U
24. scalar k2
    & _- f. u) r* s4 U
25. scalar k31 o: e# x9 O) F, W+ D
    
26. scalar k4  V( y5 s. q' ^1 t1 J
    
27. 
    # ^\" k) b# _3 m+ R. A
28. '定义权重的过程量  W; @% J  r. j* d' |, H7 w: a
    
29. scalar w1  ?& X. [1 N) K5 R5 ]
    
30. scalar w2& O: C$ t  B3 j3 n) R( B% ?& y
    
31. scalar w3% c: C. C8 z% c7 C
    
32. scalar w40 O/ r9 y& t% K) i\" x) X$ l
    
33. 
    , A3 T, v; ~* @. B
34. '程序主体
    - m, n6 R, ~\" b# f  y% t
35. for !k=1 to M step 1- A# ^6 v6 ^# V3 v  i) Y
    
36.   F(!k+1,1)=F(!k,1)+h  n! E  o! S& @6 S# t\" L
    
37.   '调用常微分方程计算权重; i/ a- `0 |7 j4 p\" T
    
38.   call obj(w1,F(!k,1),F(!k,2))
    6 Q+ h2 Z  p: _; ]\" ~1 i0 ?5 r\" b
39.     k1=w1*h
    / w$ T\" m; b* \4 _
40.   call obj(w2,F(!k,1)+h/2,F(!k,2)+k1/2)! q# {8 o\" x- e5 O* D
    
41.     k2=w2*h: ~- Q- N\" Y: A8 J+ Z2 u0 E
    
42.   call obj(w3,F(!k,1)+h/2,F(!k,2)+k2/2)2 v0 o) Q/ M5 j
    
43.     k3=w3*h* e$ W8 W* a0 w7 s
    
44.   call obj(w4,F(!k,1)+h,F(!k,2)+k3). V$ a. H- R3 r  X8 r0 w1 B; q
    
45.     k4=w4*h
    ) p+ [7 \' Z2 [& {/ e) M/ w
46.   '计算函数估计值3 q, r1 l8 L& K/ W6 `% z8 l
    
47.   F(!k+1,2)=F(!k,2)+(k1+2*k2+2*k3+k4)/6
    ' R; O2 Q: v# z\" X; V: [
48.   '计算函数解析值
    / B) n4 k: v9 f. o8 k\" g7 C
49.   F(!k+1,3)=@exp(-sin(F(!k+1,1)))*F(!k+1,1)
    ( Q* N% F1 _: H
50. next+ _8 z4 x5 N. z, Q& p
    
51. 
    , |6 A5 J; o8 o; i# G
52. '显示最终结果
    8 M9 B, a% N: H  g9 ?
53. freeze F.xyline\" s; u; |8 ^\" p# o\" w6 _) f\" q# d
    
54. freeze F' _* R: d- L3 [: s% w# X
    
55. $ E- B3 v+ l+ n$ v( z; D
    
56. '定义常微分方程' c  X# D$ n1 B) s) ]\" F4 S
    
57. subroutine obj(scalar dydx,scalar x,scalar y)7 \0 ~5 [! g' }
    
58.   dydx=-y*@cos(x)+@exp(-@sin(x))
    ( Y\" a$ J+ z2 M6 @( u* X
59. endsub: b. [* v# w& l& Z* v; r
    

复制代码

运行后求得结果如下：
- m! h( ?' s9 h1 a/ K% _# t
7 g# d  c1 T! d2 w1 q

5 {0 B6 E" Y! b( M2 S
5 |+ x- ^' Y4 \+ G其中C2列是数值解，C3列是解析解，比较之下，这二者之间无明显差异。
' n" n7 O% y: E, d  M# X; f



* A- }1 G5 Y( t1 h$ o% y  V