十常微分方程组求解Word下载.docx
《十常微分方程组求解Word下载.docx》由会员分享，可在线阅读，更多相关《十常微分方程组求解Word下载.docx（53页珍藏版）》请在冰豆网上搜索。
的数值解，分别取
，并计算误差，画出精确解和数值解的图形.
解编写并保存名为Eulerli1.m的MATLAB计算和画图的主程序如下
functionP=Eulerli1（x0,y0,b,h>
n=（b-x0>
/h。
X=zeros（n,1>
Y=zeros（n,1>
k=1。
X（k>
=x0。
Y（k>
=y0。
fork=1:
n
X（k+1>
=X（k>
+h。
Y（k+1>
=Y（k>
+h*（X（k>
-Y（k>
k=k+1。
end
y=X-1+2*exp（-X>
plot（X,Y,'
mp'
X,y,'
b-'
grid
xlabel（'
自变量X'
ylabel（'
因变量Y'
title（'
用向前欧拉公式求dy/dx=x-y，y（0>
=1在[0,1]上的数值解和精确解y=x-1+2exp（-x>
legend（'
h=0.075时，dy/dx=x-y，y（0>
=1在[0,1]上的数值解'
精确解y=x-1+2exp（-x>
jwY=y-Y。
xwY=jwY./y。
k1=1:
n。
k=[0,k1]。
P=[k'
X,Y,y,jwY,xwY]。
在MATLAB工作窗口输入下面的程序
x0=0。
y0=1。
b=1。
h=0.0750。
P=Eulerli1（x0,y0,b,h>
h1=0.0075。
P1=Eulerli1（x0,y0,b,h1>
h1=0.0075时，dy/dx=x-y，y（0>
10.3.2向前欧拉方法的三种MATLAB程序
（一>
向前欧拉公式的MATLAB主程序1
用向前欧拉公式<
10.8）求解常微分方程初值问题<
10.5）的数值解及其截断误差公式的MATLAB主程序1
function[h,k,X,Y,P,REn]=Qeuler1（funfcn,x0,y0,b,n,tol>
x=x0。
h=（b-x>
/n。
y=y0。
X（k>
=x。
=y'
fork=2:
n+1
fxy=feval（funfcn,x,y>
delta=norm（h*fxy,'
inf'
wucha=tol*max（norm（y,'
1.0>
ifdelta>
=wucha
x=x+h。
y=y+h*fxy。
Y（k>
plot（X,Y,'
rp'
label（'
用向前欧拉<
Euler）公式计算dy/dx=f（x,y>
，y（x0>
=y0在[x0,b]上的数值解'
P=[X,Y]。
symsdy2,
REn=0.5*dy2*h^2。
例10.3.2用向前欧拉公式<
10.8）求解初值问题
，
分别取
，并将计算结果与精确解作比较，写出在每个子区间
上的局部截断误差公式，画出数值解与精确解在区间
上的图形.
解输入程序
subplot（2,1,1>
x0=0。
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,'
用向前欧拉公式计算dy/dx=8x-3y-7，y（0>
n=100时，dy/dx=8x-3y-7，y（0>
dy/dx=8x-3y-7，y（0>
=1在[0,1]上的精确解'
holdoff
jdwuc1=S1-Y1。
jwY1=S1-Y1。
xwY1=jwY1./S1。
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>
S1=8/3*x2-29/9+38/9*exp（-3*x2>
plot（x2,S1,'
n=10时，dy/dx=8x-3y-7，y（0>
jwY2=S1-Y2。
xwY2=jwY2./S1。
n1。
P2=[k'
x2,Y2,S1,jwY2,xwY2]
运行后屏幕显示分别取
时，所给的初值问题在
上的自变量
的值构成的数组Xi（i=1,2>
，利用向前欧拉公式<
10.8）求出的与Xi（i=1,2>
对应的数值解Yi（i=1,2>
，Yi的绝对误差jwYi（i=1,2>
和相对误差xwYi（i=1,2>
<
略），每个子区间
上的局部截断误差公式Ren2=1/200*dy2，近似解Y与精确解的图形.
（二>
向前欧拉公式的MATLAB主程序2
10.5）的数值解及其截断误差公式的MATLAB主程序2
function[k,X,Y,P,REn]=Qeuler2（funfcn,x0,y0,b,h>
n=fix（（b-x>
/h>
X=zeros（n+1,1>
Y=zeros（n+1,1>
=x+（k-1>
*h,
fxy=feval（funfcn,x,y>
=y+h*fxy
y=Y（k>
k=k+1,plot（X,Y,'
grid,xlabel（'
P=[k'
X,Y]。
例10.3.4用向前欧拉公式<
，取
b=2。
n=10。
h=2/10。
[k,X,Y,P,REn]=Qeuler2（@funfcn,x0,y0,b,h>
S1=1/6*（6+12*X+30*exp（2*X>
.^（1/2>
plot（X,S1,'
用向前欧拉公式计算dy/dx=y-x/（3y>
，y（0>
=1在[0,2]上的数值解'
n=10时，dy/dx=y-x/（3y>
dy/dx=y-x/（3y>
=1在[0,2]上的精确解'
jdwucY=S1-Y。
jwY=S1-Y。
xwY=jwY./Y。
X,Y,S1,jwY,xwY]
n1=100。
h1=2/100。
[k,X1,Y1,P1,Ren1]=Qeuler2（@funfcn,x0,y0,b,h1>
S2=1/6*（6+12*X1+30*exp（2*X1>
plot（X1,S2,'
n=100时，dy/dx=y-x/（3y>
jwY1=S2-Y1。
xwY1=jwY1./Y1。
P2=[k'
X1,Y1,S2,jwY1,xwY1]
运行后屏幕显示取
时，此初值问题在
的向量X，X1，利用向前欧拉公式<
10.8）求出的与X，X1对应的数值解Y，Y1及其绝对误差向量jwY，jwY1和相对误差向量xwY，xwY1<
上的局部截断误差公式Ren2=1/200*dy2，数值解Y与精确解的图形.
（三>
自适应向前欧拉公式的MATLAB主程序
用自适应向前欧拉公式求解常微分方程初值问题<
10.5）的数值解的MATLAB主程序2
function[H,X,Y,k,h,P]=QEuler（funfcn,x0,b,y0,tol>
%初始化.
pow=1/3。
ifnargin<
5|isempty（tol>
tol=1.e-6。
end。
6|isempty（trace>
trace=0。
h=0.0078125*（b-x>
y=y0（:
p=128。
H=zeros（p,1>
X=zeros（p,1>
Y=zeros（p,length（y>
Y（k,:
%绘图.
clc,x,h,y
%主循环.
while（x<
b>
&
（x+h>
x>
ifx+h>
b
h=b-x。
end
%计算斜率.
fxy=fxy（:
%计算误差，设定可接受误差.
%当误差可接受时重写解.
ifdelta<
ifk>
length（X>
X=[X。
zeros（p,1>
]。
Y=[Y。
zeros（p,length（y>
H=[H。
H（k>
=h。
Y（k,:
grid
%更新步长.
ifdelta~=0.0
h=min（h*8,0.9*h*（wucha/delta>
^pow>
if（x<
disp（'
Singularitylikely.'
x
H=H（1:
k>
X=X（1:
Y=Y（1:
k,:
n=1:
k。
P=[n'
H,X,Y]
10.3.4向后欧拉方法的MATLAB程序
用向后欧拉方法求解常微分方程初值问题<
10.5）的数值解的MATLAB主程序
function[X,Y,n,P]=Heuler1（funfcn,x0,b,y0,h,tol>
n=fix（（b-x0>
X=zeros（n+1,1>
Y=zeros（n+1,1>
k=1。
Y1（k,:
clc,x0,h,y0
%产生初值.
fori=2:
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,:
p>
Y（i,:
=p。
x0=x0+h。
y0=Y1（i,:
ro'
gridon
用向后欧拉公式计算dy/dx=f（x,y>
n+1>
Y=Y（1:
n+1,:
n=1:
n+1。
X,Y]
例10.3.6用向后欧拉公式求解区间
上的初值问题
的数值解，取步长
，并与精确解作比较，在同一个坐标系中作出图形.然后再取
，观察数值解与精确解误差的变化，说明
与误差的关系.
解输入程序
S1=dsolve（'
Dy=8*x-3*y-7'
y（0>
=1'
x'
tol=1.e-1。
h1=0.01。
[X1,Y1,n,P1]=Heuler1（@funfcn,x0,b,y0,h1,tol>
S2=8/3*X1-29/9+38/9*exp（-3*X1>
plot（X1,S2,'
h=0.01用向后欧拉公式计算dy/dx=8x-3y-7，y（0>
dy/dx=8x-3y-7，y（0>
juwY1=S2-Y1。
xiwY1=juwY1./Y1。
L=[P1,S2,juwY1,xiwY1]
h=0.05。
[X,Y,n,P]=Heuler1（@funfcn,x0,b,y0,h,tol>
S1=8/3*X-29/9+38/9*exp（-3*X>
plot（X,S1,'
h=0.05用向后欧拉公式计算dy/dx=8x-3y-7，y（0>
juwY=S1-Y。
xiwY=juwY./Y。
L=[P,S1,juwY,xiwY]
运行后屏幕显示用向后欧拉公式计算此初值问题在[0,2]上的自变量X处数值解Y和精确解S1及其图形，步长H,Y的相对误差xiwY和绝对误差juwY（略>
.
10.4改进的欧拉方法的MATLAB程序
10.4.2梯形公式的两种MATLAB程序
（一）梯形公式的MATLAB程序
用梯形公式求解常微分方程初值问题<
10.5）的数值解的MATLAB主程序1
function[X,Y,n,P]=odtixing1（funfcn,x0,b,y0,h,tol>
n+1
fx0y0=feval（funfcn,x0,y0>
=y0+h*fx0y0。
fxiyi=feval（funfcn,X（i>
=y0+h*（fxiyi+fx0y0>
/2。
fxip=feval（funfcn,X（i>
=y0+h*（fx0y0+fxip>
/2,P1=Y1（i,:
Y（i,:
=p1。
用梯形公式计算dy/dx=f（x,y>
例10.4.1用梯形公式求解区间
，取步长
，精度为
，并与精确解作比较，在同一个坐标系中画出图形.
b=2。
tol=0.1。
h=0.05。
[X,Yt,n,Pt]=odtixing1（@funfcn,x0,b,y0,h,tol>
S1=8/3*X-29/9+38/9*exp（-3*X>
holdoff
h=0.05,用梯形公式计算dy/dx=8x-3y-7，y（0>
juwYt=S1-Yt。
xiwYt=juwYt./Yt。
Lt=[Pt,S1,juwYt,xiwYt]
运行后屏幕显示取精度为
，分别用梯形公式和向前欧拉公式求解此初值问题在区间
上的自变量X处数值解Yi（i=t,q>
和精确解S1，步长H,Yi的相对误差xiwYi和绝对误差juwYi（略>
及其数值解和精确解的图形.
自适应梯形公式的MATLAB程序
用自适应梯形公式求解常微分方程初值问题<
function[H,X,Y,k,h,P]=odtixing2（funfcn,x0,b,y0,tol>
n=fix（（b-x0>
ifx+h>
%计算斜率.
y1=y+h*fxy。
fxy1=feval（funfcn,x,y1>
y2=y+h*fxy1。
y=（y1+y2>
go'
用自适应梯形公式计算dy/dx=f（x,y>
例10.4.2分别用自适应梯形公式和向前欧拉公式分别求解区间
，取精度为
，并与精确解作比较，在同一个坐标系中作出图形.
tol=1.e-1。
[Ht,X,Yt,k,h,Pt]=odtixing2（@funfcn,x0,b,y0,tol>
holdon,
[Hq,X,Yq,k,h,Pq]=QEuler（@funfcn,x0,b,y0,tol>
用自适应梯形公式计算dy/dx=8x-3y-7，y（0>
dy/dx=8x-3y-7，y（0>
用自适应梯形公式计算dy/dx=8x-3y-7,y（0>
和精确解S1及其图形，步长H等（略>
.
10.4.4改进的欧拉公式的MATLAB程序
用自适应改进的欧拉公式求解常微分方程初值问题<
function[H,X,Y,k,h,P]=gaiEuler（funfcn,x0,b,y0,tol>
trace=0。
y=y0（:
H=zeros（p,1>
X=zeros（p,1>
Y=
升级会员 