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

liwenhui

1 \$ f- o' N/ ^, f( _/ F, |) p
EViews除了能解决计量经济学的估计问题以外，还提供一个编程环境用以解决复杂的问题。经过调试，我在EViews中实现了用龙格库塔方法求解常微分方程的数值解，供大家交流。
! ]5 g6 A, b/ ~/ {2 x$ w- y演示中，我使用了如下常微分方程作为测试：
7 Z5 m3 e, g7 L- j) H

( l& D. |3 {: L# z  E这个方程的解析通解是：
2 y) N; A+ d" Y9 E

3 {5 O' ~/ y: K0 k+ t/ L2 Q& _( S使用“龙格库塔方法”，编制的EViews程序如下：

1. '用龙格库塔法求解常微分方程dydx=-y*cos(x)+exp(-sin(x)),y(0)=0在区间[0,10*pi]的数值解
   1 E\" A& M% O) A  Z4 I6 _$ y: B
2. '已知这个微分方程的解析解 y=exp(-sin(x))*x, Q1 I+ i3 o, p! c: Q$ Z
   
3. 
   5 C: G7 N9 t3 S& @\" P
4. '生成一个workfile作为基本的数据容器# c! e  m$ Z! D6 A% P$ A\" J
   
5. wfcreate (wf=temp) u 10004 _% \, j( s& G2 H8 d) n
   
6. 8 O0 L2 B7 J- p  l
   
7. '定义常量# I, j, N: j2 ^! L
   
8. scalar pi=3.14159
   . O$ F9 j' S& y8 m! x/ J% q, O2 z
9. scalar a=0        '定义自变量下限
   0 T. N8 t- j1 V
10. scalar b=10*pi     '定义自变量上限. U9 Z+ p0 Z8 d8 K7 e: ?5 r- H
    
11. scalar  M=500       '定义步数: g\" o0 G5 L2 |' D0 I4 L
    
12. scalar h=(b-a)/M   '计算每步之间的间隔
      w2 R/ K* F0 n& I' d7 W! o* U# S
13. 
      ~$ m. Z! g$ }/ ^: U! F
14. '定义一个矩阵来储存计算数据，其中第一列储存自变量数据，第二列储存因变量数据,第三列储存解析解的值用以作为比较
    ' I- L) g- W+ r/ {( L+ D) s
15. matrix(M+1,3) F
    ! r* b9 l2 g' i0 P, m
16. 
    - P. x9 w3 f1 [& `2 M$ I+ a7 o
17. '矩阵的第一行储存初值问题的初始条件8 T0 L6 d  {  X) X6 P
    
18. F(1,1)=0% C8 Z( e! m6 v
    
19. F(1,2)=0
    : b% W. \! g( \
20. F(1,3)=@exp(-sin(F(1,1)))*F(1,1)1 m. b2 w! s) U( _
    
21. 7 |! B4 p% ?1 R0 X$ d: e
    
22. '定义龙格库塔法的权重参数3 G* t( X* f% C0 e1 c4 `  ^
    
23. scalar k1
    % a\" _0 o) x) ?  Q9 ]- t
24. scalar k2
    + H0 Q8 G$ H1 B/ |1 h
25. scalar k3
    ; Y  r# O: G  u
26. scalar k4! @* E# j- J( L9 G5 I
    
27. 
    8 m* z2 ?0 U/ |, {1 c  P& Z
28. '定义权重的过程量
    7 I5 R' Z2 L: K9 D8 J$ N
29. scalar w1
    4 |& H& ^) ]/ \\" `; n: k5 y: v
30. scalar w2+ L8 o\" D/ A1 P/ w) `/ b7 P\" t
    
31. scalar w3. V7 V1 q. ~9 U7 j; c+ J( b
    
32. scalar w4
    0 o5 l7 Y\" v/ s8 F, X- D) s3 `
33. 
    , h* S0 ?2 V' r6 M
34. '程序主体6 N% I9 I- Z, T, b
    
35. for !k=1 to M step 1  q  ]8 d* ^7 _4 ^\" C. P
    
36.   F(!k+1,1)=F(!k,1)+h; e+ U9 a\" j* s
    
37.   '调用常微分方程计算权重
    4 N' e$ a( r. u: \2 z/ i- G  d
38.   call obj(w1,F(!k,1),F(!k,2)) 8 p7 f. G\" g& \6 T\" M\" O, k
    
39.     k1=w1*h1 a4 C/ t; X7 O& A. N: O+ b
    
40.   call obj(w2,F(!k,1)+h/2,F(!k,2)+k1/2)
    0 _5 W* c$ ]6 |8 Y
41.     k2=w2*h6 {& `( z* k# f3 Z' [0 h
    
42.   call obj(w3,F(!k,1)+h/2,F(!k,2)+k2/2)
    7 J& c1 T# q, U
43.     k3=w3*h
    9 W& R8 {/ R; ?
44.   call obj(w4,F(!k,1)+h,F(!k,2)+k3)( Q6 W9 a' d- f0 `$ G
    
45.     k4=w4*h) `  f/ a+ W0 o7 s( D+ h7 w- W( E
    
46.   '计算函数估计值% l, F8 {5 V/ C  P$ ]. E& c
    
47.   F(!k+1,2)=F(!k,2)+(k1+2*k2+2*k3+k4)/68 P1 b3 U8 Q7 d7 i. Z( w2 t
    
48.   '计算函数解析值9 F5 ]7 h. s# x8 D4 Z: B
    
49.   F(!k+1,3)=@exp(-sin(F(!k+1,1)))*F(!k+1,1)4 V+ F2 Z3 x- c. E  V; o
    
50. next' j- B4 h+ _& @! Q1 J% X
    
51. 6 F& D\" Q1 s- y
    
52. '显示最终结果
    + @: _\" b\" b# `/ V- t
53. freeze F.xyline
    7 u\" q8 q# ~3 b
54. freeze F
    1 t6 ]2 @% L2 q: O
55. ; T5 W4 d1 d# c+ l% S# G
    
56. '定义常微分方程% H$ v: U3 i) T' y# a( ?! ]
    
57. subroutine obj(scalar dydx,scalar x,scalar y)\" i. N4 G6 `  J/ z
    
58.   dydx=-y*@cos(x)+@exp(-@sin(x))
    # c$ Q& R: j' V% v\" J$ o\" ]
59. endsub1 B# B: x2 E( s2 R
    

复制代码

运行后求得结果如下：
, I6 f7 f/ ^. W& B. S% g
* u  Q. S2 D6 f, W* I6 c

2 F8 f: ^' ^' {8 q$ C. U' V. ?
6 u4 y3 O1 }. S其中C2列是数值解，C3列是解析解，比较之下，这二者之间无明显差异。
- x+ J' u/ b7 u0 [( ]
( A  l% q6 X( R7 I/ N& u
- E. a+ r2 R- s7 n
% `% _5 W2 h( L* @
4 a, _$ R& U1 K4 w& J