在EViews中实现数值求解常微分方程（ODE）
本帖最后由 liwenhui 于 2016-12-6 15:41 编辑EViews除了能解决计量经济学的估计问题以外，还提供一个编程环境用以解决复杂的问题。经过调试，我在EViews中实现了用龙格库塔方法求解常微分方程的数值解，供大家交流。
演示中，我使用了如下常微分方程作为测试：
这个方程的解析通解是：
使用“龙格库塔方法”，编制的EViews程序如下：'用龙格库塔法求解常微分方程dydx=-y*cos(x)+exp(-sin(x)),y(0)=0在区间的数值解
'已知这个微分方程的解析解 y=exp(-sin(x))*x
'生成一个workfile作为基本的数据容器
wfcreate (wf=temp) u 1000
'定义常量
scalar pi=3.14159
scalar a=0 '定义自变量下限
scalar b=10*pi '定义自变量上限
scalar M=500 '定义步数
scalar h=(b-a)/M '计算每步之间的间隔
'定义一个矩阵来储存计算数据，其中第一列储存自变量数据，第二列储存因变量数据,第三列储存解析解的值用以作为比较
matrix(M+1,3) F
'矩阵的第一行储存初值问题的初始条件
F(1,1)=0
F(1,2)=0
F(1,3)=@exp(-sin(F(1,1)))*F(1,1)
'定义龙格库塔法的权重参数
scalar k1
scalar k2
scalar k3
scalar k4
'定义权重的过程量
scalar w1
scalar w2
scalar w3
scalar w4
'程序主体
for !k=1 to M step 1
F(!k+1,1)=F(!k,1)+h
'调用常微分方程计算权重
call obj(w1,F(!k,1),F(!k,2))
k1=w1*h
call obj(w2,F(!k,1)+h/2,F(!k,2)+k1/2)
k2=w2*h
call obj(w3,F(!k,1)+h/2,F(!k,2)+k2/2)
k3=w3*h
call obj(w4,F(!k,1)+h,F(!k,2)+k3)
k4=w4*h
'计算函数估计值
F(!k+1,2)=F(!k,2)+(k1+2*k2+2*k3+k4)/6
'计算函数解析值
F(!k+1,3)=@exp(-sin(F(!k+1,1)))*F(!k+1,1)
next
'显示最终结果
freeze F.xyline
freeze F
'定义常微分方程
subroutine obj(scalar dydx,scalar x,scalar y)
dydx=-y*@cos(x)+@exp(-@sin(x))
endsub
运行后求得结果如下：
其中C2列是数值解，C3列是解析解，比较之下，这二者之间无明显差异。
页:
[1]