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

liwenhui

( g* ^& k$ u/ o/ s6 x8 P# }
EViews除了能解决计量经济学的估计问题以外，还提供一个编程环境用以解决复杂的问题。经过调试，我在EViews中实现了用龙格库塔方法求解常微分方程的数值解，供大家交流。
4 q* j' U& V: c, Q( |演示中，我使用了如下常微分方程作为测试：
4 P6 _- ]3 ]) @. ~* X, T

) ]. R; i/ ~# Y0 r0 N# ~' a9 x
这个方程的解析通解是：
& \* i2 z4 h$ s9 A

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

1. '用龙格库塔法求解常微分方程dydx=-y*cos(x)+exp(-sin(x)),y(0)=0在区间[0,10*pi]的数值解
   , k3 _; N, E( z. T: x. o
2. '已知这个微分方程的解析解 y=exp(-sin(x))*x
   ' W$ \& }! F/ u
3. 
   3 R' `% X/ p' Q# q3 a* T
4. '生成一个workfile作为基本的数据容器
   5 F) H4 R3 ?& H: A\" n! s  M2 f
5. wfcreate (wf=temp) u 1000+ e/ p7 ?# d0 a- z! r+ ^7 R/ {
   
6. 
   1 O$ {$ r. |5 y7 T, Q! ]
7. '定义常量  E1 v9 R1 k\" {1 Y! v! u+ A
   
8. scalar pi=3.14159
   ; o8 ?' U9 f9 I* u$ W
9. scalar a=0        '定义自变量下限
   % q/ F* [+ ~7 T( J1 J, I# c
10. scalar b=10*pi     '定义自变量上限
    ! n1 h9 ?( A, _8 v) W
11. scalar  M=500       '定义步数8 s1 V7 z1 Y\" o* j: K
    
12. scalar h=(b-a)/M   '计算每步之间的间隔; n; T6 b' ?8 s9 f5 M7 \4 ?2 ^
    
13. 
    + j\" w+ h$ v! A. E' Q& g
14. '定义一个矩阵来储存计算数据，其中第一列储存自变量数据，第二列储存因变量数据,第三列储存解析解的值用以作为比较9 x; s0 [4 A8 R& M/ e. o
    
15. matrix(M+1,3) F   e6 T7 q  U( Q( g4 ]2 M6 z
    
16. \" t5 h. O& a+ D; B
    
17. '矩阵的第一行储存初值问题的初始条件: b' [/ {6 E$ d- m5 A; S2 n
    
18. F(1,1)=0
    : L# H7 d8 G\" e9 K) K) G, X+ D
19. F(1,2)=0* L0 H3 Y6 N% Z0 p& X8 T# h) L
    
20. F(1,3)=@exp(-sin(F(1,1)))*F(1,1)# z  G% b7 F3 o. i* s$ H
    
21. 
    & A: p3 O+ ~3 k# X7 U: h  b, J
22. '定义龙格库塔法的权重参数\" Q$ y/ I# ~* f% B3 |
    
23. scalar k1
    . D0 K9 r6 ~) ?& i  M7 [2 t
24. scalar k2& I- }# f% {: V\" y# [% z
    
25. scalar k3
    $ \. K  C9 m9 e) g8 ?; K5 V
26. scalar k4
    1 Z8 V8 W9 s4 V' A  J
27. 
    6 o2 e' S. g4 l8 A/ _% l- o: S
28. '定义权重的过程量( b/ d4 g! V* s/ t+ y0 f. G, q
    
29. scalar w1- f9 U- J: U8 t/ F  m  R
    
30. scalar w2
    ( j  c, p- S# ~  |3 r9 e. i
31. scalar w3+ a( z$ H1 g7 O. {\" z: J
    
32. scalar w4
    9 i4 }! y, K7 L: R0 m6 D3 j9 f
33. 5 |* `3 J0 o9 z5 v2 ?8 `9 B
    
34. '程序主体& ~  C# G4 M4 H! z. h: N
    
35. for !k=1 to M step 1+ K1 K% l' U+ i% k1 |, {; C
    
36.   F(!k+1,1)=F(!k,1)+h
    # ^0 z3 i' t2 i8 v; k& E. t
37.   '调用常微分方程计算权重' s8 m( E5 [4 H& L. a2 \
    
38.   call obj(w1,F(!k,1),F(!k,2)) ! p7 [4 q2 S- X6 _
    
39.     k1=w1*h0 v, j2 `# e! n! Z
    
40.   call obj(w2,F(!k,1)+h/2,F(!k,2)+k1/2)
    - r+ j& X! o6 T2 ]% e
41.     k2=w2*h
    7 ~; Z4 J0 x. b
42.   call obj(w3,F(!k,1)+h/2,F(!k,2)+k2/2)- ~0 s6 l# }\" S, }! Y5 r
    
43.     k3=w3*h3 J/ U5 }$ D8 ]( @: D
    
44.   call obj(w4,F(!k,1)+h,F(!k,2)+k3)
    2 v- {5 ]5 Y. e$ Q1 t
45.     k4=w4*h
    % b7 F% P$ [/ w* l
46.   '计算函数估计值
    + L2 p8 r5 Z  g7 Q; o% @! Y% w
47.   F(!k+1,2)=F(!k,2)+(k1+2*k2+2*k3+k4)/6
    / }/ _( E- i( _! f
48.   '计算函数解析值
    , v0 h- G' [$ O0 k0 R& u' q9 c5 L
49.   F(!k+1,3)=@exp(-sin(F(!k+1,1)))*F(!k+1,1)
    % H3 e0 j* b9 ?- @& Y/ u
50. next9 e; w6 R1 S6 H1 b\" x( h' `
    
51. . v- ?; Z9 l  H5 `: V
    
52. '显示最终结果
    : P3 s$ q* H* X& A8 H6 H
53. freeze F.xyline3 R9 P# j\" u& t. a+ _% U; {2 m5 a
    
54. freeze F\" r; }. B9 {  a( t2 X
    
55. 5 h7 G3 K9 e( X  C5 y$ C3 C
    
56. '定义常微分方程
    \" e. L! P& i; R% U9 p; v
57. subroutine obj(scalar dydx,scalar x,scalar y)! {4 {* G8 I3 U: s/ G2 Q6 `. M' |
    
58.   dydx=-y*@cos(x)+@exp(-@sin(x))2 T4 j3 J: X; [6 ]8 _
    
59. endsub
    1 q4 X\" a2 U0 ?. E\" B' s

复制代码

运行后求得结果如下：

6 C, {0 R' m7 @6 Z

6 }7 B/ X' A) [+ ~# |# R# T其中C2列是数值解，C3列是解析解，比较之下，这二者之间无明显差异。
& }. i7 I( @" W. v# c. c

8 q7 j# g+ \2 F, u+ w% `
% g! \* D+ q* N6 F2 g/ s' p
7 q, I  }3 t3 s1 I