第十章 常微分方程组求解.docx
- 文档编号:2993173
- 上传时间:2023-05-05
- 格式:DOCX
- 页数:34
- 大小:70.05KB
第十章 常微分方程组求解.docx
《第十章 常微分方程组求解.docx》由会员分享,可在线阅读,更多相关《第十章 常微分方程组求解.docx(34页珍藏版)》请在冰点文库上搜索。
第十章常微分方程组求解
9.1常微分方程(组)的MATLAB符号求解
9.1.1用MATLAB求常微分方程(组)的通解
调用格式一:
S=dsolve('eqn','var')
调用格式二:
S=dsolve('eqn1','eqn2',...,'eqnm','var')
9.1.2用MATLAB求常微分方程(组)的特解
调用格式三:
S=dsolve('eqn','condition1',…,'conditionn','var')
调用格式四:
S=dsolve('eqn1','eqn2',...,'eqnm','condition1','condition2',…,'var')
9.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;
9.2欧拉(Euler)方法的MATLAB程序
9.2.1向前欧拉公式及其误差估计
例9.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)')
9.2.2向前欧拉方法的兩种MATLAB程序
(一)向前欧拉公式的MATLAB主程序1
用向前欧拉公式求解常微分方程初值问题的数值解及其截断误差公式的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;
例9.2.2用向前欧拉公式求解初值问题
,
分别取
,并将计算结果与精确解作比较,写出在每个子区间
上的局部截断误差公式,画出数值解与精确解在区间
上的图形.
解输入程序
>>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),利用向前欧拉公式求出的与Xi(i=1,2)对应的数值解Yi(i=1,2),Yi的绝对误差jwYi(i=1,2)和相对误差xwYi(i=1,2)(略),每个子区间
上的局部截断误差公式Ren2=1/200*dy2,近似解Y与精确解的图形.
(二)向前欧拉公式的MATLAB主程序2
用向前欧拉公式求解常微分方程初值问题的数值解及其截断误差公式的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;
例9.2.3用向前欧拉公式求解初值问题
,取
,并将计算结果与精确解作比较,写出在每个子区间
上的局部截断误差公式,画出数值解与精确解在区间
上的图形.
解输入程序
>>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,利用向前欧拉公式求出的与X,X1对应的数值解Y,Y1及其绝对误差向量jwY,jwY1和相对误差向量xwY,xwY1(略),每个子区间
上的局部截断误差公式Ren2=1/200*dy2,数值解Y与精确解的图形.
9.2.3向后欧拉方法的MATLAB程序
用向后欧拉方法求解常微分方程初值问题的数值解的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]
例9.2.4用向后欧拉公式求解区间
上的初值问题
的数值解,取步长
,并与精确解作比较,在同一个坐标系中作出图形.然后再取
,观察数值解与精确解误差的变化,说明
与误差的关系.
解输入程序
>>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(略).
9.3改进的欧拉方法的MATLAB程序
9.3.1梯形公式的MATLAB程序
(一)梯形公式的MATLAB程序
用梯形公式求解常微分方程初值问题的数值解的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]
例9.3.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(略)及其数值解和精确解的图形.
例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等(略).
9.3.2改进的欧拉公式的MATLAB程序
用自适应改进的欧拉公式求解常微分方程初值问题的数值解的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=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(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]
例9.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
[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]
运行后屏幕显示取精度为
,分别用梯形公式和改进的欧拉公式求此解初值问题在区间
上的自变量X处数值解Yt,Y和精确解S1及其图形,步长H,Y的相对误差xiwY和绝对误差juwY(略).
9.5龙格-库塔方法
9.5.1二阶龙格-库塔方法的MATLAB程序
用二阶龙格—库塔方法求解常微分方程初值问题的数值解的MATLAB主程序
function[k,X,Y,fxy,wch,wucha,P]=RK2(funfcn,fun,x0,b,C,y0,h)
x=x0;y=y0;p=128;
n=fix((b-x0)/h);
fxy=zeros(p,1);
wucha=zeros(p,1);
wch=zeros(p,1);
X=zeros(p,1);
Y=zeros(p,length(y));
k=1;X(k)=x;Y(k,:
)=y';
%绘图.
clc,x,h,y
%计算
%fxy=fxy(:
);
fork=2:
n+1
x=x+h;a2=C(3);b21=C(4);
c1=C
(1);
c2=C
(2);
x1=x+a2*h;
k1=feval(funfcn,x,y);
y1=y+b21*h*k1;
k2=feval(funfcn,x1,y1);
fxy(k)=feval(fun,x);
y=y+h*(c1*k1+c2*k2);
X(k)=x;
Y(k,:
)=y;k=k+1;
plot(X,Y,'mh',X,fxy,'bo')
grid,xlabel('自变量X'),ylabel('因变量Y')
legend('用二阶龙格-库塔方法计算dy/dx=f(x,y),y(x0)=y0在[x0,b]上的数值解','y/dx=f(x,y),y(x0)=y0的精确解y=f(x)')
end
%计算误差.
fork=2:
n+1
wucha(k)=norm(Y(k-1)-Y(k));
wch(k)=norm(fxy(k)-Y(k));
end
X=X(1:
k);
Y=Y(1:
k,:
);
fxy=fxy(1:
k,:
);
n=1:
k;
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 第十章 常微分方程组求解 第十 微分 方程组 求解