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

liwenhui

7 m4 K0 z0 b4 O+ N6 \; U5 _EViews除了能解决计量经济学的估计问题以外，还提供一个编程环境用以解决复杂的问题。经过调试，我在EViews中实现了用龙格库塔方法求解常微分方程的数值解，供大家交流。
演示中，我使用了如下常微分方程作为测试：
' l. C7 b8 H& A2 U) T) I

' o) m9 i8 ^1 t: H& A$ n1 h
8 [1 U' P: `9 s/ j3 V' H! F" P这个方程的解析通解是：

6 A% I  }1 r! V7 L/ L9 T
使用“龙格库塔方法”，编制的EViews程序如下：

1. '用龙格库塔法求解常微分方程dydx=-y*cos(x)+exp(-sin(x)),y(0)=0在区间[0,10*pi]的数值解3 @1 T7 [4 _9 B\" F0 M+ A8 p# j& C
   
2. '已知这个微分方程的解析解 y=exp(-sin(x))*x
   4 E/ Q4 O- r8 [/ ?0 Z7 w
3.   r! K- r3 D5 ?
   
4. '生成一个workfile作为基本的数据容器$ b; g) g; s2 Y: f2 a+ B
   
5. wfcreate (wf=temp) u 1000/ `- g* c0 N& [8 d7 t: U
   
6. 
   ( w3 c3 ^; H! o( t
7. '定义常量+ u, F8 W& N* w& e: d/ \, B\" l
   
8. scalar pi=3.14159\" e* U3 |7 X' h. c  i
   
9. scalar a=0        '定义自变量下限( o\" Q0 N! a% V, F
   
10. scalar b=10*pi     '定义自变量上限
    . b$ E) i* D6 Q$ j1 N' y
11. scalar  M=500       '定义步数
    0 e# ^# ?. S- l: B* f6 z% w
12. scalar h=(b-a)/M   '计算每步之间的间隔
    3 q/ _; v8 G) C1 q8 l$ u
13. $ s' [- n, R\" S& }, h8 m) l& z# h
    
14. '定义一个矩阵来储存计算数据，其中第一列储存自变量数据，第二列储存因变量数据,第三列储存解析解的值用以作为比较/ c% l  Z$ p$ M( d
    
15. matrix(M+1,3) F # C6 p2 L& ~\" {
    
16. 
      S$ I% L2 T  A1 [' }
17. '矩阵的第一行储存初值问题的初始条件
    7 u+ Q9 f  N* z9 m: ^# d$ J  @% y
18. F(1,1)=0) Z6 h) A7 m/ q: W
    
19. F(1,2)=0/ ~, |2 f$ \# a6 L% c. @6 T! ~$ o
    
20. F(1,3)=@exp(-sin(F(1,1)))*F(1,1)& I( s. @, }8 e) w
    
21. & d! ?# i# o; R. L0 @
    
22. '定义龙格库塔法的权重参数
    , v\" T8 t1 Q) R
23. scalar k1
    : Q) x) a5 u\" r5 r
24. scalar k2
    ! M, q) Q; k9 v: H
25. scalar k3
    9 ]\" M; [3 @, B\" q2 V
26. scalar k4
    - k% X7 X( I1 M' P# M0 u& x
27. 
    * @2 C* C8 C7 @
28. '定义权重的过程量& Q% W4 E9 `, d
    
29. scalar w1
    ' E( i* O8 X6 n' [\" _- A
30. scalar w2
    8 C1 n+ z# M: Y
31. scalar w3; A) n3 \% W3 O
    
32. scalar w4& r7 M6 ]5 Y# e3 d5 a# I5 O) p
    
33. 
    + B) Z+ l1 W$ W& l0 L# t
34. '程序主体- u: K7 o% x% l/ G\" o- y# l
    
35. for !k=1 to M step 1- j4 h2 Z* J; U- D3 E
    
36.   F(!k+1,1)=F(!k,1)+h9 @* b6 X) e3 [; v
    
37.   '调用常微分方程计算权重$ b( @7 T& M9 J1 s# Q9 o
    
38.   call obj(w1,F(!k,1),F(!k,2)) 0 Q! O4 S+ S) x\" b4 d9 L! g- R
    
39.     k1=w1*h
    1 L1 C8 }. `1 b# `
40.   call obj(w2,F(!k,1)+h/2,F(!k,2)+k1/2)  w9 {, q( F/ {( n
    
41.     k2=w2*h# M3 @& ]1 s* J! t# f  O# w
    
42.   call obj(w3,F(!k,1)+h/2,F(!k,2)+k2/2)* G$ S+ T- V: P
    
43.     k3=w3*h
    $ t$ L5 B6 w\" _
44.   call obj(w4,F(!k,1)+h,F(!k,2)+k3)
    ) E! ^) N8 r. x% d3 F\" S; v% x( W8 @
45.     k4=w4*h
    % A' w' y+ j4 [. G
46.   '计算函数估计值
    : k\" _5 B% K; f$ \+ N
47.   F(!k+1,2)=F(!k,2)+(k1+2*k2+2*k3+k4)/6
    ! X+ Y4 ^  Y: F; H0 @
48.   '计算函数解析值
    5 W) ~& B9 R: a
49.   F(!k+1,3)=@exp(-sin(F(!k+1,1)))*F(!k+1,1)
    0 v/ D; @, {1 _4 V5 @+ @
50. next
    # A! K/ q! L  c/ e, C, g
51. 9 Q' m: V5 s) v- n* W
    
52. '显示最终结果7 |\" f  M. F7 ]0 I
    
53. freeze F.xyline
    ' h+ e! n. |1 b9 }* E( w\" r
54. freeze F& [0 ^7 F4 w  z8 \8 n# [0 \
    
55. # m/ u. H( p6 @\" Z
    
56. '定义常微分方程. J* P4 ~* V$ @/ F4 B) \8 r$ Q
    
57. subroutine obj(scalar dydx,scalar x,scalar y)* U( u$ W: T& e: @3 b) x1 i
    
58.   dydx=-y*@cos(x)+@exp(-@sin(x))
    , O+ ^9 I1 w6 J0 ^' Z, L  j6 P
59. endsub
    # b( v8 l; N! y. W( ?  q; _

复制代码

运行后求得结果如下：
/ V9 H+ N7 r0 L" y

其中C2列是数值解，C3列是解析解，比较之下，这二者之间无明显差异。

7 ?+ j( O9 j: _( M/ d3 b8 R
( u, L2 r- V/ H, W

. @( F( i9 x- u