(2)常微分方程的差分方法.doc

实验三：常微分方程的差分方法程序设计与验证
李郭应数0810809501024
实验结果
1.改进的欧拉方法
（1）建立euler函数文件(euler.m)
euler.m文件程序如下：
function[x,y]=euler(dyfun,xspan,y0,h)
%dyfun为函数f(x,y)
%xspan为求解区间[x0,xN]
%h为步长
%x返回节点
%y返回数值解
x=xspan(1):h:xspan(2);
y(1)=y0;
forn=1:length(x)-1
k1=feval(dyfun,x(n),y(n));
y(n+1)=y(n)+h*k1;
k2=feval(dyfun,x(n+1),y(n+1));
y(n+1)=y(n)+h*(k1+k2)/2;
end
x=x';
y=y';
（2）在工作窗口输入如下程序
tic
formatlong
dyfun=inline('y-2*x/y');
[x,y]=euler(dyfun,[0,1],1,0.2);
[x,y]
toc
（3）结果
步长时，结果为：
01.0000
0.10001.0959
0.20001.1841
0.30001.2662
0.40001.3434
0.50001.4164
0.60001.4860
0.70001.5525
0.80001.6165
0.90001.6782
1.00001.7379
步长时，结果为：
01.0000
0.20001.1867
0.40001.3483
0.60001.4937
0.80001.6279
1.00001.7542
2.四阶经典龙格-库塔方法
(1)建立nark4函数文件
nark4.m程序如下
function[x,y]=nark4(dyfun,xspan,y0,h)
%dyfun为函数f(x,y)
%xspan为求解区间[x0,xN]
%h为步长
%x返回节点
%y返回数值解
x=xspan(1):h:xspan(2);
y(1)=y0;
forn=1:length(x)-1
k1=feval(dyfun,x(n),y(n));
k2=feval(dyfun,x(n)+h/2,y(n)+h/2*k1);
k3=feval(dyfun,x(n)+h/2,y(n)+h/2*k2);
k4=feval(dyfun,x(n+1),y(n)+h*k3);
y(n+1)=y(n)+h*(k1+2*k2+2*k3+k4)/6;
end
x=x'；
y=y';
（2）在工作窗口输入如下程序
dyfun=inline('y-2*x/y');
[x,y]=nark4(dyfun,[0,1],1,0.1);
[x,y]
（3）结果
步长时，结果为：
01.0000
0.10001.0954
0.20001.1832
0.30001.2649
0.40001.3416
0.50001.4142
0.60001.4832
0.70001.5492
0.80001.6125
0.90001.6733
1.00001.7321
步长时，结果为：
01.0000
0.20001.1832
0.40001.3417
0.60001.4833
0.80001.6125
1.00001.7321

分析讨论
由结果可知，改进的欧拉公式和四阶经典龙格库塔方法均能求解常微分方程的初值问题，对于改进的欧拉格式，得到的结果比更精确，同时，四阶经典龙格库塔方法也是如此，且对于实验所给的题目，时，四阶经典龙格库塔方法得到的结果和真实值最为接近。
aa=1.73205080756888;
plot(x1,y1,'r',x2,y2,'b')
holdon
plot(aa,'mh')
legend('euler','nark4','精确值')
gridon

图一：时，两种差分方法结果与真实值之间的比较
步长为时，两种常微分方程的差分方法结果比较
同理可以得到步长为时，两种常微分方程的差分方法结果比较的图形。

图二：时，两种差分方法结果与真实值之间的比较