微分方程数值.docx
《微分方程数值.docx》由会员分享，可在线阅读，更多相关《微分方程数值.docx（14页珍藏版）》请在冰点文库上搜索。
微分方程数值
微分方程数值解
指导教师：
李晓峰
学院：
应用科学学院
班级：
信息与计算科学131801
姓名：
王慧兰
学号：
201318030120
目录
一、实验名称错误!
未定义书签。
二、实验目的错误!
未定义书签。
三、实验内容错误!
未定义书签。
四、实验要求错误!
未定义书签。
五、实验步骤、源程序及运行截图错误!
未定义书签。
1.Euler法（向前）错误!
未定义书签。
1.1源程序错误!
未定义书签。
1.2实验步骤错误!
未定义书签。
1.3运行截图错误!
未定义书签。
2.Euler法（向后）错误!
未定义书签。
2.1源程序错误!
未定义书签。
2.2实验步骤错误!
未定义书签。
2.3运行截图错误!
未定义书签。
3.梯形迭代法公式错误!
未定义书签。
3.1源程序错误!
未定义书签。
3.2实验步骤错误!
未定义书签。
3.3运行截图错误!
未定义书签。
4.改进的Euler公式错误!
未定义书签。
4.1源程序10
4.2实验步骤错误!
未定义书签。
4.3运行截图11
一、实验名称
Euler法、梯形法求解初值问题。
二、实验目的
掌握求解常微分方程初值问题的单步法；比较不同算法所得的数值结果，体会步长缩小对问题解得影响。
三、实验内容
求解初值问题
的数值解。
四、实验要求
编写程序，步长取h=0.1，分别用Euler法、梯形法迭代格式、4阶RK方法求解该问题。
整理比较各节点处的数值解、误差、相对误差、体会算法对数值解梯度的改进。
（1）改变步长h=0.2和h=0.5，通过比较各节点处数值解，在不同的步长下，算法优劣结论一致。
（2）按同一算法不同步长整理数据，比较在同一节点按不同步长计算所得的数值解的误差变化趋势，体会缩小步长对数值解的影响。
（3）改进Euler法程序使其适用任意初值问题。
五、源程序、实验核心步骤运行截图
1、Euler法（向前）
1.1源程序
QEuler.m
function[h,k,X,Y,P,REn]=Qeuler1（funfcn,x0,y0,b,n,tol）
x=x0;h=（b-x）/n;X=zeros（n,1）;y=y0;
Y=zeros（n,1）;k=1;X（k）=x;Y（k）=y';
fork=2:
n+1
fxy=feval（funfcn,x,y）;
delta=norm（h*fxy,'inf'）;
wucha=tol*max（norm（y,'inf'）,1.0）;
ifdelta>=wucha
x=x+h;y=y+h*fxy;X（k）=x;Y（k）=y';
end
plot（X,Y,'rp'）
grid
label（'自变量X'）,ylabel（'因变量Y'）
title（'用向前欧拉（Euler）公式计算dy/dx=f（x,y），y（x0）=y0在[x0,b]上的数值解'）
end
P=[X,Y];
symsdy2,
REn=0.5*dy2*h^2;
1.2实验步骤
（1）建立并保存以funfcn.m文件命名的M文件函数
functionf=funfcn（x,y）
f=8*x-3*y-7;
（2）建立并保存以QEuler.m文件命名的M文件函数。
（3）输入程序：
subplot（2,1,1）
x0=0;y0=1;b=1-1.e-4;
n=100;tol=1.e-4;
[h1,k1,x1,Y1,P1,Ren1]=QEuler1（@funfcn,x0,y0,b,n,tol）
holdon
S1=8/3*x1-29/9+38/9*exp（-3*x1）,plot（x1,S1,'b-'）
title（'用向前欧拉公式计算dy/dx=8x-3y-7，y（0）=1在[0,1]上的数值解'）
legend（'n=100时，dy/dx=8x-3y-7，y（0）=1在[0,1]上的数值解','dy/dx=8x-3y-7，y（0）=1在[0,1]上的精确解'）
holdoff
jdwuc1=S1-Y1;jwY1=S1-Y1;
xwY1=jwY1./S1;k1=1:
n;k=[0,k1];
P1=[k',x1,Y1,S1,jwY1,xwY1]
subplot（2,1,2）
n1=10;[h2,k2,x2,Y2,P2,Ren2]=QEuler1（@funfcn,x0,y0,b,n1,tol）
holdon
S1=8/3*x2-29/9+38/9*exp（-3*x2）,plot（x2,S1,'b-'）
legend（'n=10时，dy/dx=8x-3y-7，y（0）=1在[0,1]上的数值解','dy/dx=8x-3y-7，y（0）=1在[0,1]上的精确解'）
holdoff
jwY2=S1-Y2;xwY2=jwY2./S1;k1=1:
n1;k=[0,k1];P2=[k',x2,Y2,S1,jwY2,xwY2]
1.3运行截图
图1-3
=10,100时，
在[0,1]上的数值解和精确解
2.向后Euler公式：
2.2源程序：
Heuler1.m
function[X,Y,n,P]=Heuler1（funfcn,x0,b,y0,h,tol）
n=fix（（b-x0）/h）;X=zeros（n+1,1）;Y=zeros（n+1,1）;
k=1;X（k）=x0;Y（k,:
）=y0;Y1（k,:
）=y0;
%绘图.
clc,x0,h,y0
%产生初值.
fori=2:
n+1
X（i）=x0+h;Y（i,:
）=y0+h*feval（funfcn,x0,y0）;
Y1（i,:
）=y0+h*feval（funfcn,X（i）,Y（i,:
））;
%主循环.
Wu=abs（Y1（i,:
）-Y（i,:
））;
whileWu>tol
p=Y1（i,:
）;
Y1（i,:
）=y0+h*feval（funfcn,X（i）,p）;
Y（i,:
）=p;
end
x0=x0+h;y0=Y1（i,:
）;
Y（i,:
）=y0;plot（X,Y,'ro'）
gridon
xlabel（'自变量X'）,ylabel（'因变量Y'）
title（'用向后欧拉公式计算dy/dx=f（x,y），y（x0）=y0在[x0,b]上的数值解'）
end
X=X（1:
n+1）;Y=Y（1:
n+1,:
）;n=1:
n+1;P=[n',X,Y]
2.1实验步骤：
（1）建立并保存以funfcn.m文件命名的M文件函数
functionf=funfcn（x,y）
f=8*x-3*y-7;
（2）建立并保存以Heuler1.m文件命名的M文件函数。
（3）输入程序
>>S1=dsolve（'Dy=8*x-3*y-7','y（0）=1','x'）
>>x0=0;y0=1;
b=2;tol=1.e-1;
subplot（2,1,1）
h1=0.01;
[X1,Y1,n,P1]=Heuler1（@funfcn,x0,b,y0,h1,tol）
holdon
S2=8/3*X1-29/9+38/9*exp（-3*X1）,plot（X1,S2,'b-'）
legend（'h=0.01用向后欧拉公式计算dy/dx=8x-3y-7，y（0）=1在[0,2]上的数值解','dy/dx=8x-3y-7，y（0）=1在[0,2]上的精确解'）
holdoff
juwY1=S2-Y1;
xiwY1=juwY1./Y1;
L=[P1,S2,juwY1,xiwY1]
subplot（2,1,2）
h=0.05;
[X,Y,n,P]=Heuler1（@funfcn,x0,b,y0,h,tol）
holdon
S1=8/3*X-29/9+38/9*exp（-3*X）,
plot（X,S1,'b-'）
legend（'h=0.05用向后欧拉公式计算dy/dx=8x-3y-7，y（0）=1在[0,2]上的数值解','dy/dx=8x-3y-7，y（0）=1在[0,2]上的精确解'）
holdoff
juwY=S1-Y;
xiwY=juwY./Y;
L=[P,S1,juwY,xiwY]
2.3实验结果图：
图10–10取h=0.05时，用向后欧拉公式求初值问题的数值解和精确解的图形
3、梯形迭代法公式
3.1源程序
odtixing1.m
function[X,Y,n,P]=odtixing1（funfcn,x0,b,y0,h,tol）
n=fix（（b-x0）/h）;
X=zeros（n+1,1）;
Y=zeros（n+1,1）;
k=1;X（k）=x0;
Y（k,:
）=y0;Y1（k,:
）=y0;
%绘图.
clc,x0,h,y0
%产生初值.
fori=2:
n+1
X（i）=x0+h;
fx0y0=feval（funfcn,x0,y0）;
Y（i,:
）=y0+h*fx0y0;
fxiyi=feval（funfcn,X（i）,Y（i,:
））;
Y1（i,:
）=y0+h*（fxiyi+fx0y0）/2;
%主循环.
Wu=abs（Y1（i,:
）-Y（i,:
））;
whileWu>tol
p=Y1（i,:
）,
fxip=feval（funfcn,X（i）,p）;
Y1（i,:
）=y0+h*（fx0y0+fxip）/2,P1=Y1（i,:
）,Y（i,:
）=p1;
end
x0=x0+h;y0=Y1（i,:
）;
Y（i,:
）=y0;plot（X,Y,'ro'）
gridon
xlabel（'自变量X'）,ylabel（'因变量Y'）
title（'用梯形公式计算dy/dx=f（x,y），y（x0）=y0在[x0,b]上的数值解'）
end
X=X（1:
n+1）;Y=Y（1:
n+1,:
）;
n=1:
n+1;P=[n',X,Y]
3.2实验步骤
（1）建立并保存以funfcn.m文件命名的M文件函数
functionf=funfcn（x,y）
f=8*x-3*y-7;
（2）建立并保存以QEuler.m文件命名的M文件函数。
（3）输入程序：
x0=0;y0=1;b=2;tol=0.1;h=0.05;
[X,Yt,n,Pt]=odtixing1（@funfcn,x0,b,y0,h,tol）
holdon
S1=8/3*X-29/9+38/9*exp（-3*X）;
plot（X,S1,'b-'）,holdoff
legend（'h=0.05,用梯形公式计算dy/dx=8x-3y-7，y（0）=1在[0,2]上的数值解','dy/dx=8x-3y-7，y（0）=1在[0,2]上的精确解'）
juwYt=S1-Yt;xiwYt=juwYt./Yt;Lt=[Pt,S1,juwYt,xiwYt]
3.3运行截图
4、改进的Euler公式
4.1源程序：
gaiEuler.m
function[H,X,Y,k,h,P]=gaiEuler（funfcn,x0,b,y0,tol）
%初始化.
pow=1/3;
ifnargin<5|isempty（tol）,
tol=1.e-6;
end;
ifnargin<6|isempty（trace）,
trace=0;
end;
x=x0;h=0.0078125*（b-x）;
y=y0（:
）;p=128;n=fix（（b-x0）/h）;
H=zeros（p,1）;X=zeros（p,1）;
Y=zeros（p,length（y））;k=1;
X（k）=x;Y（k,:
）=y';
%绘图.
clc,x,h,y
end
%主循环.
while（xx）
ifx+h>b
h=b-x;
end
%计算斜率.
fxy=feval（funfcn,x,y）;fxy=fxy（:
）;
%计算误差，设定可接受误差.
delta=norm（h*fxy,'inf'）;wucha=tol*max（norm（y,'inf'）,1.0）;
%当误差可接受时重写解.
ifdelta<=wucha
x=x+h;y1=y+h*fxy;fxy1=feval（funfcn,x,y1）;
fxy=fxy（:
）;y2=（fxy+fxy1）/2;y=y+h*y2;k=k+1;
ifk>length（X）
X=[X;zeros（p,1）];Y=[Y;zeros（p,length（y））];
H=[H;zeros（p,1）];
end
H（k）=h;X（k）=x;Y（k,:
）=y';plot（X,Y,'mh'）,grid
xlabel（'自变量X'）,ylabel（'因变量Y'）
title（'用改进的欧拉公式计算dy/dx=f（x,y），y（x0）=y0在[x0,b]上的数值解'）
end
%更新步长.
ifdelta~=0.0
h=min（h*8,0.9*h*（wucha/delta）^pow）;
end
end
if（xdisp（'Singularitylikely.'）,x
end
H=H（1:
k）;X=X（1:
k）;Y=Y（1:
k,:
）;n=1:
k;P=[n',H,X,Y]
4.2实验步骤：
（1）建立并保存以funfcn.m文件命名的M文件函数
functionf=funfcn（x,y）
f=8*x-3*y-7;
（2）建立并保存以gaiEuler.m文件命名的M文件函数。
（3）输入程序
x0=0;y0=1;b=2;tol=1.e-1;
[Ht,X,Yt,k,h,Pt]=odtixing2（@funfcn,x0,b,y0,tol）
holdon
S1=8/3*X-29/9+38/9*exp（-3*X）,
plot（X,S1,'b-'）
holdoff
holdon
[H,X,Y,k,h,P]=gaiEuler（@funfcn,x0,b,y0,tol）
holdoff
legend（'用梯形公式计算dy/dx=8x-3y-7，y（0）=1在[0,2]上的数值解','dy/dx=8x-3y-7，y（0）=1在[0,2]上的精确解','用改进的欧拉公式计算dy/dx=8x-3y-7，y（0）=1在[0,2]上的数值解'）
juwY=S1-Y;xiwY=juwY./Y;
L=[P,S1,juwY,xiwY]
3.3运行截图
升级会员 