十章 常微分方程组求解.docx
- 文档编号:14400258
- 上传时间:2023-06-23
- 格式:DOCX
- 页数:61
- 大小:178.31KB
十章 常微分方程组求解.docx
《十章 常微分方程组求解.docx》由会员分享,可在线阅读,更多相关《十章 常微分方程组求解.docx(61页珍藏版)》请在冰点文库上搜索。
十章常微分方程组求解
10.1常微分方程(组)的MATLAB符号求解
10.1.1用MATLAB求常微分方程(组)的通解
调用格式一:
S=dsolve('eqn','var')
调用格式二:
S=dsolve('eqn1','eqn2',...,'eqnm','var')
10.1.2用MATLAB求常微分方程(组)的特解
调用格式三:
S=dsolve('eqn','condition1',…,'conditionn','var')
调用格式四:
S=dsolve('eqn1','eqn2',...,'eqnm','condition1','condition2',…,'var')
10.1.3线性常微分方程组的解法
求解齐次线性常微分方程组
的MATLAB主程序
名为qcxxcwz.m的求解齐次线性常微分方程组
的MATLAB主程序如下:
function[X,E,V]=qcxxcwz(A)
symsx
E=eig(A);
[V,n]=eig(A);
X=exp(E*x)'*V;
10.3欧拉(Euler)方法的MATLAB程序
10.3.1向前欧拉公式及其误差估计
例10.2.1用欧拉方法求初值问题
的数值解,分别取
,并计算误差,画出精确解和数值解的图形.
解编写并保存名为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)
在MATLAB工作窗口输入下面的程序
>>h1=0.0075;P1=Eulerli1(x0,y0,b,h1)
legend('h1=0.0075时,dy/dx=x-y,y(0)=1在[0,1]上的数值解','精确解y=x-1+2exp(-x)')
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;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;
例10.3.2用向前欧拉公式(10.8)求解初值问题
,
分别取
,并将计算结果与精确解作比较,写出在每个子区间
上的局部截断误差公式,画出数值解与精确解在区间
上的图形.
解输入程序
>>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]
运行后屏幕显示分别取
时,所给的初值问题在
上的自变量
的值构成的数组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.8)求解常微分方程初值问题(10.5)的数值解及其截断误差公式的MATLAB主程序2
function[k,X,Y,P,REn]=Qeuler2(funfcn,x0,y0,b,h)
x=x0;n=fix((b-x)/h);
X=zeros(n+1,1);y=y0;
Y=zeros(n+1,1);k=1;
X(k)=x;Y(k)=y';
fork=2:
n+1
X(k)=x+(k-1)*h,
fxy=feval(funfcn,x,y),
Y(k)=y+h*fxy
y=Y(k),k=k+1,plot(X,Y,'rp'),
grid,xlabel('自变量X'),ylabel('因变量Y')
title('用向前欧拉(Euler)公式计算dy/dx=f(x,y),y(x0)=y0在[x0,b]上的数值解')
end
k1=1:
n;k=[0,k1];P=[k',X,Y];
symsdy2,
REn=0.5*dy2*h^2;
例10.3.4用向前欧拉公式(10.8)求解初值问题
,取
,并将计算结果与精确解作比较,写出在每个子区间
上的局部截断误差公式,画出数值解与精确解在区间
上的图形.
解输入程序
>>subplot(2,1,1)
x0=0;y0=1;b=2;n=10;h=2/10;
[k,X,Y,P,REn]=Qeuler2(@funfcn,x0,y0,b,h)
holdon
S1=1/6*(6+12*X+30*exp(2*X)).^(1/2);plot(X,S1,'b-')
title('用向前欧拉公式计算dy/dx=y-x/(3y),y(0)=1在[0,2]上的数值解')
legend('n=10时,dy/dx=y-x/(3y),y(0)=1在[0,2]上的数值解','dy/dx=y-x/(3y),y(0)=1在[0,2]上的精确解')
holdoff
jdwucY=S1-Y;jwY=S1-Y;
xwY=jwY./Y;k1=1:
n;k=[0,k1];
P1=[k',X,Y,S1,jwY,xwY]
subplot(2,1,2)
n1=100;h1=2/100;
[k,X1,Y1,P1,Ren1]=Qeuler2(@funfcn,x0,y0,b,h1)
holdon
S2=1/6*(6+12*X1+30*exp(2*X1)).^(1/2);plot(X1,S2,'b-')
legend('n=100时,dy/dx=y-x/(3y),y(0)=1在[0,2]上的数值解','dy/dx=y-x/(3y),y(0)=1在[0,2]上的精确解')
holdoff
jwY1=S2-Y1;xwY1=jwY1./Y1;k1=1:
n1;k=[0,k1];
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;
ifnargin<6|isempty(trace),
trace=0;
end;
x=x0;h=0.0078125*(b-x);y=y0(:
);p=128;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;y=y+h*fxy;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,'rp'),grid
xlabel('自变量X'),ylabel('因变量Y')
title('用向前欧拉(Euler)公式计算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(x
disp('Singularitylikely.'),x
end
H=H(1:
k);
X=X(1:
k);
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)/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]
例10.3.6用向后欧拉公式求解区间
上的初值问题
的数值解,取步长
,并与精确解作比较,在同一个坐标系中作出图形.然后再取
,观察数值解与精确解误差的变化,说明
与误差的关系.
解输入程序
>>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]
运行后屏幕显示用向后欧拉公式计算此初值问题在[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=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]
例10.4.1用梯形公式求解区间
上的初值问题
,取步长
,精度为
,并与精确解作比较,在同一个坐标系中画出图形.
解输入程序
>>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]
运行后屏幕显示取精度为
,分别用梯形公式和向前欧拉公式求解此初值问题在区间
上的自变量X处数值解Yi(i=t,q)和精确解S1,步长H,Yi的相对误差xiwYi和绝对误差juwYi(略)及其数值解和精确解的图形.
(二)自适应梯形公式的MATLAB程序
用自适应梯形公式求解常微分方程初值问题(10.5)的数值解的MATLAB主程序
function[H,X,Y,k,h,P]=odtixing2(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
%主循环.
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=y+h*fxy1;y=(y1+y2)/2;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,'go'),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(x
disp('Singularitylikely.'),x
end
H=H(1:
k);
X=X(1:
k);
Y=Y(1:
k,:
);
n=1:
k;P=[n',H,X,Y]
例10.4.2分别用自适应梯形公式和向前欧拉公式分别求解区间
上的初值问题
,取精度为
,并与精确解作比较,在同一个坐标系中作出图形.
解输入程序
>>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,
[Hq,X,Yq,k,h,Pq]=QEuler(@funfcn,x0,b,y0,tol),
holdoff
grid
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]上的数值解')
title('用自适应梯形公式计算dy/dx=8x-3y-7,y(0)=1在[0,2]上的数值解')
运行后屏幕显示取精度为
,分别用梯形公式和向前欧拉公式求解此初值问题在区间
上的自变量X处数值解Yi(i=t,q)和精确解S1及其图形,步长H等(略).
10.4.4改进的欧拉公式的MATLAB程序
用自适应改进的欧拉公式求解常微分方程初值问题(10.5)的数值解的MATLAB主程序
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=
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 十章 常微分方程组求解 微分 方程组 求解