第十章 常微分方程组求解.docx
- 文档编号:15323116
- 上传时间:2023-07-03
- 格式:DOCX
- 页数:22
- 大小:33.64KB
第十章 常微分方程组求解.docx
《第十章 常微分方程组求解.docx》由会员分享,可在线阅读,更多相关《第十章 常微分方程组求解.docx(22页珍藏版)》请在冰点文库上搜索。
第十章常微分方程组求解
第十章常微分方程(组)求解
第三篇第十章常微分方程求解 Matlab常微分方程求解一、求微分方程的解 相关函数及简介 1,dsolve(‘equ1’,’equ2’,…):
Matlab求微分方程的解析解。
equ1,equ2,…为方程。
写方程时用Dy表示y关于自变量的一阶导数,用D2y表示y关于自变量的二阶导数,依次类推。
2,simplify(s):
对表达式s使用maple的化简规则进行化简。
例如:
symsx simplify(sin(x)_+cos(x)_)ans=1 3,[r,how]=simple(s):
于Matlab提供了多种化简规则,simple命令就是对表达式s用各种规则进行化简,然后用r返回最简形式,how返回形成这种形式所用的规则。
例如:
symsx [r,how]=simple(cos(x)_-sin(x)_)r=cos(2*x)how=combine 4,[T,Y]=solver(odefun,tspan,y0),求微分方程的数值解。
其中的solver为命令 ode45,ode23,ode113,ode15s,ode23s,ode23t,ode23tb之一。
?
dy?
f(t,y),odefun是显式常微分方程:
?
?
dt?
?
y(t0)?
y0.在积分区间tspan?
[t0,tf]上,从t0到tf,用初始条件y0求解。
要获得问题在其他指定时间点t0,t1,t2,...上的解,则令tspan?
[t0,t1,t2,...,tf] 因为没有一种算法可以有效地解决所有的ODE问题,为此,Matlab提供了多种求解器solver,对于不同的ODE问题,采用不同的solver。
求解器solverODE类型 特点 说明 ode45 非刚性 单步算法:
4,5阶的Runge-Kutta 大部分场合的首选算法 方程;累计截断误差达(?
x)3 ode23 非刚性 单步算法:
2,3阶的Runge-Kutta 使用于精度较低的情形 3 方程;累计截断误差达(?
x) 170. 第三篇第十章常微分方程求解 ode113 非刚性 多步法:
Adams算法;高低精度均 计算时间比ode45短可到10~10。
ode23t 适度刚性 采用梯形算法 适度刚性情形 ode15s 刚性 多步法:
Gear?
s反向数值微分 若ode45失败时,可尝 精度中等 试使用 ode23s 刚性 单步法:
2阶Rosebrock算法; 当精度较低时,计算 低精度 时间比ode15s短 ode23tb 刚性 梯形算法;低精度 当精度较低时,计算 时间比ode15s短 要特别说明的是:
ode23,ode45是极其常用的用来求解非刚性的标准形式的一阶常 微分方程的初值问题的解的Matlab的常用程序,其中:
ode23采用龙格-库塔2阶算法,用3阶公式作误差估计来调节步长,具有低等的精度。
ode45则采用龙格-库塔4阶算法,用5阶公式作误差估计来调节步长,具有中等的精度。
5,ezplot(x,y,[tmin,tmax]):
符号函数的作图命令,x,y为关于参数t的符号函数,[tmin,tmax] 为t的取值范围。
6,inline():
建立一个内联函数。
格式:
inline(‘expr’,’var1’,’var2’,…),注意括号里的表达式要 加引号。
如:
Q=dblquad(inline(‘y*sin(x)’),pi,2*pi,0,pi). ?
3?
6例1:
求解微分方程dy?
2xy?
xe?
x,并加以验证。
2dx求解本问题的Matlab程序:
symsxy y=dsolve(‘Dy+2*x*y=x*exp(-x_)’,’x’)simplify(diff(y,x)+2*x*y-x*exp(-x_))说明:
行1是用命令定义x,y为符号变量。
这里可以 不写,但为确保正确性,建议写上; 行2是用命令求出的微分方程的解:
1/2*exp(-x_)*x_+exp(-x_)*C1 行3使用所求得的解。
这里是将解代入原微分 方程,结果应该为0,但这里给出:
-x_*exp(-x_)-2*x*exp(-x_)*C1+2*x*(1/2*exp(-x_)*x_+exp(-x_)*C1) 行4用simplify()函数对上式进行化简,结果为 0,表明y=y(x)的确是微分方程的解。
例2:
求解微分方程xy?
?
y?
ex?
0在初始条件y
(1)?
2e下的特解,并画出解函数的图形。
求解本问题的Matlab程序:
symsxy y=dsolve(‘x*Dy+y-exp(x)=0’,’y
(1)=2*exp
(1)’,’x’)ezplot(y) 171. 第三篇第十章常微分方程求解 微分方程的特解为y=1/x*exp(x)+1/x*exp
(1)(Matlab格式), e?
ex即y?
,解函数图形。
x?
dxt?
5x?
y?
e,?
?
dt例3:
求解微分方程组?
在初始条件xt?
0?
1,yt?
0?
0dy?
?
x?
3y?
0?
?
dt下的特解,并画出解函数的图形。
Matlab程序:
symsxy [x,y]=dsolve(‘Dx+5*x+y=exp(t)’,’Dy-x-3*y=0’,’x(0)=1’,’y(0)=0’,’t’) simple(x);simple(y); ezplot(x,y,[0,]);axisauto2,用ode23,ode45等求解非刚性的标准形式的一阶常微分方程的初值问题的数值解例 ?
dy2?
?
?
2y?
2x?
2x,4:
求解微分方程初值问题?
dx的数值解,求解范围 ?
?
y(0)?
1为[0,]。
Matlab程序:
fun=inline(‘-2*y+2*x_+2*x’,’x’,’y’);[x,y]=ode23(fun,[0,],1);x’;y’; plot(x,y,’o-’)例5:
用Euler 2x?
dy?
?
y?
2,折线法求解微分方程初值问题?
dxy的数 ?
y(0)?
1?
值解,求解范围为[0,2]。
解:
本问题的差分方程为 ?
x0?
0,y0?
1,h?
?
?
xk?
1?
xk?
h,?
2?
yk?
1?
yk?
hf(xk,yk)(其中:
f(x,y)?
y?
2x/y),k?
0,1,2,...,n?
1Matlab程序:
clear f=sym(‘y+2*x/y_‘);a=0; 172. 第三篇第十章常微分方程求解 b=2;h=; n=(b-a)/h+1;x=0;y=1; szj=[x,y];fori=1:
n-1 y=y+h*subs(f,{‘x’,’y’},{x,y}); x=x+h; szj=[szj;x,y];endszj plot(szj(:
1)szj(:
2))练习:
1,求微分方程(x2?
1)y?
?
2xy?
sinx?
0的通解。
2,求微分方程y?
?
?
2y?
?
5y?
exsinx的通解。
?
dx?
x?
y?
0,?
dt3,求微分方程组?
在初值条件xt?
0?
1,yt?
0?
0下的特解,并画?
?
dy?
x?
y?
0?
?
dt出函数y?
f(x)的图形。
4,分别用ode23,ode45求上述第3题中的微分方程初值问题的数值解,求解区间为t?
[0,2]。
利用画图来比较两种求解器之间的差异。
5,用Euler ?
12x2?
y?
?
y?
3,折线法求解微分方程初值问题?
y的数值 ?
y(0)?
1?
解,求解范围为[0,2]。
?
y?
?
y?
excosx,6,用四阶Runge?
Kutta法求解微分方程初值问题?
的数 ?
y(0)?
1值解,求解范围为[0,3]。
四阶Runge?
Kutta法的迭代公式为求解 ?
?
y0?
y(x0),?
?
xk?
1?
xk?
h,?
h?
yk?
1?
yk?
(L1?
2L2?
2L3?
L4),6?
?
k?
0,1,2,...,n?
1法):
?
L1?
f(xk,yk),?
hh?
L2?
f(xk?
yk?
L1),22?
?
hh?
L3?
f(xk?
yk?
L2),22?
?
?
L4?
f(xk?
h,yk?
hL3),试用该方法求解第5题中的初值问题。
Matlab程序:
clear f=sym(‘y-exp(x)*cos(x)’);a=0;b=3;h=; n=(b-a)/h+1;x=0;y=1; szj=[x,y];fori=1:
n-1 l1=subs(f,{‘x’,’y’},{x,y}); l2=subs(f,{‘x’,’y’},{x+h/2,y+l1*h/2});l3=subs(f,{‘x’,’y’},{x+h/2,y+l2*h/2});l4=subs(f,{‘x’,’y’},{x+h,y+l3*h});y=y+h*(l1+2*12+2*13+14)/6; x=x+h; szj=[szj;x,y];endszj plot(szj(:
1)szj(:
2)) 7,用ode45方法求上述第6题的常微分方程初值问题的数值解,从而利用画图来比较两者的差异。
二、常微分方程的MATLAB符号求解 用MATLAB求常微分方程的通解 用Matbal函数dsolve求常微分方程:
F(x,y,y?
y?
?
...,y(n))?
0 174.
第三篇第十章常微分方程求解 holdoffjdwuc1=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] 运行后屏幕显示分别取n?
10,100时,所给的初值问题在[0,1]上的自变量x的值构成的数组Xi(i=1,2),利用向前欧拉公式求出的与Xi(i=1,2)对应的数值解Yi(i=1,2),Yi的绝对误差jwYi(i=1,2)和相对误差xwYi(i=1,2),每个子区间[xk,xk?
1]上的局部截断误差公式Ren2=1/200*dy2,近似解Y与精确解的图形. 例分别用向前欧拉公式()、一、四阶泰勒逼近及用Matlab函数double求解 ?
dy?
?
1?
ycosx,0?
x?
2,初值问题?
dx分别取n?
40,80时,并将计算结果与精确解做比 ?
?
y(0)?
0.较,写出在每个子区间?
xn,xn?
1?
上的局部截断误差公式,画出数值解与精确解在区间[0, 2]上的图形。
解:
(1)建立并保存名为的M文件函数。
functionf=function(x,y) f=1-y*cos(x);
(2)输入程序 >>S1=dsolve(‘Dy=1-y*cos(x)’,’y(0)=0’,’x’) 取一阶泰勒逼近 y?
e?
sinx(1?
x?
cosx). 同理,取四阶泰勒逼近,输入程序>>symsX y=exp(-sin(X))*int(1+sin(u)+((sin(u))_)/2+((sin(u))_)/6+((sin(u))_)/24,u,0,X) 根据运行后,将屏幕显示结果整理得四阶泰勒逼近为 ?
108111?
1017?
?
y?
e?
sinx?
?
x?
cosx?
?
sinx?
sin2x?
sin3x?
?
. 1896?
964?
?
?
964取n?
40时,输入程序 >>symsu h=2/40,X=0:
h:
2; S=exp(-sin(X))*int(exp(sin(u)),u,0,X);S1=double(S);subplot(4,1,1) plot(X,S1,’bo’),grid legend(‘用double计算dy/dx=1-ycos(x),y(0)=0在[0,2]上的精确解的数值解’) subplot(4,1,2) y=exp(-sin(X)).*(1+X-cos(X));plot(X,y,’g-’),grid legend(‘y=exp(-sin(X)).*(1+X-cos(X))在[0,2]上的图形’)subplot(4,1,3) y1=exp(-sin(X)).*81/64*X-10/9*cos(X)-17/64*cos(X).*sin(X)-1/18*sin(X)._.*cos(X)-1/96*sin(X)._.*cos(X)+10/9);plot(X,y1,’m-’),grid legend(‘四阶泰勒逼近的解在[0,2]上的图形’)subplot(4,1,4) 180. 第三篇第十章常微分方程求解 x0=0;y0=0;b=;n=40;tol=;[h,k,X,Y,p,Ren]=QEuler1(@funfcn,x0,y0,b,n,tol)legend(‘n=40时,用向前欧拉公式计算dy/dx=1-ycos(x),y(0)=0在[0,2]上的数值解’) symsu X=0:
:
2;S=exp(-sin(X))*int(exp(sin(u)),u,0,X);S1=double(S);jwS1y1=y1’-S1’, xwS1y1=jwS1y1./S1’,jwYy1=y1’-Y,xwYy1=jwYy1./Y,k1=1:
n;k=[0,k1];p=[k’,X’,Y,y1’,jwYy1,xwYy1],pS=[k’,X’,S1’,jwS1y1],xwS1y1 运行后屏幕显示取n=40时,分别用向前欧拉公式()和一、四阶泰勒多项式逼近法及double求解此初值问题在[0,2]上的自变量X对应的数值解Y,y,y1,S1及其图形,Y与y1,S1与y1的绝对误差向量jwYy1,xwS1y1和相对误差向量xwYy1,xwS1y1.向前欧拉公式在每个子区间 ?
xk,xk?
1?
上的局部截断误差公式为 Ren?
?
?
(?
k),?
k?
?
xk,xk?
1?
k?
1,2,...,40。
取n?
40时,输入程序>>symsu X=0:
:
2; S=exp(-sin(X))*int(exp(sin(u)),u,0,X);S1=double(S);y=exp(-sin(X)).*(1+X-cos(X));plot(X,y,’g-’),gridx0=0;y0=0;b=;n=40;tol=; [h,k,X,Y,p,Ren]=QEuler1(@funfcn,x0,y0,b,n,tol) y4=exp(-sin(X)).*81/64*X-10/9*cos(X)-17/64*cos(X).*sin(X)-1/18*sin(X)._.*cos(X)-1/96*sin(X)._.*cos(X)+10/9);plot(X,S1,’mo’,X,y,’g-.’,X,Y,’rp’,X,y4,’b-’),grid legend(‘用double计算dy/dx=1-ycos(x),y(0)=0在[0,2]上的精确解的数值解’,’y=exp(-sin(X)).*(1+X-cos(X))在[0,2]上的图形’,’n=40时,用向前欧拉公式计算dy/dx=1-ycos(x),y(0)=0在[0,2]上的数值解 ’,’y4=exp(-sin(X))(10/9+81/64X-cos(X)(10/9+17/64sin(X)+1/18sin(X)_+1/96sin(X)_))在[0,2]上的图形’) 运行后屏幕显示取n=40时,分别用向前欧拉公式()和一、四阶泰勒多项式逼近法及double求解此初值问题在[0,2]上的自变量X对应的数值解Y,y,y1,S1及其图形。
取n=40时,向前欧拉公式和四阶泰勒多项式逼近法所求初值问题的数值解的曲线几乎重合,二者的绝对误差和相对误差最小。
一阶泰勒多项式逼近法所求初值问题的数值解的曲线与前两条曲线较近。
而用函数double所求初值问题的数值解的曲线与前三条曲线很远,且绝对误差和相对误差很大,尤其是相对误差,当x=0时,Y,y1,和y与真值0相等,而S1=不等于真值0,此推测,用向前欧拉公式和一、四阶泰勒多项式逼近法所求初值问题的值都接近与真值,是有实际意义的,而用函数double计算的值严重失真,没有实际意义。
下面再取n=80验证。
取n=80时,输入程序>>symsu X=0:
:
2; S=exp(-sin(X))*int(exp(sin(u)),u,0,X);S1=double(S);subplot(4,1,1) plot(X,S1,’bo’),grid legend(‘用double计算dy/dx=1-ycos(x),y(0)=0在[0,2]上的精确解的数值解’) subplot(4,1,2) y=exp(-sin(X)).*(1+X-cos(X));plot(X,y,’g-’),grid legend(‘y=exp(-sin(X)).*(1+X-cos(X))在[0,2]上的图形’)subplot(4,1,3) y1=exp(-sin(X)).*81/64*X-10/9*cos(X)-17/64*cos(X).*sin(X)-1/18*sin(X)._.*cos(X)-1/96*sin(X)._.*cos(X)+10/9);plot(X,y1,’m-’),grid legend(‘四阶泰勒逼近的解在[0,2]上的图形’)subplot(4,1,4) x0=0;y0=0;b=;n=80;tol=; [h,k,X,Y,p,Ren]=QEuler1(@funfcn,x0,y0,b,n,tol) 181. 第三篇第十章常微分方程求解 legend(‘n=80时,用向前欧拉公式计算dy/dx=1-ycos(x),y(0)=0在[0,2]上的数值解’) symsu X=0:
:
2;S=exp(-sin(X))*int(exp(sin(u)),u,0,X);S1=double(S);jwS1y=S1’-y’, xwS1y=jwS1y./S1’,jwYy=Y-y’,xwYy=jwYy./y’,k1=1:
n;k=[0,k1];p=[k’,X’,Y,y’,jwYy,xwYy],pS=[k’,X’,S1’,jwS1y],xwS1y 运行后屏幕显示取n=80时,分别用向前欧拉公式()和一、四阶泰勒多项式逼近法及double求解此初值问题在[0,2]上的自变量X对应的数值解Y,y,y1,S1及其图形,Y与y,S1与y的绝对误差向量jwYy,xwS1y和相对误差向量xwYy,xwS1y.向前欧拉公式在每个子区间 ?
xk,xk?
1?
上的局部截断误差公式为 Ren?
?
?
(?
k),?
k?
?
xk,xk?
1?
k?
1,2,...,80。
?
xk,xk?
1?
上的局部截断误差公式为 运行的结果可见,当x=0时,Y,y1,和y与真值0相等,但是S1却随着n的改变而发生变化,但是不等于真值0,此推测,用向前欧拉公式和一、四阶泰勒多项式逼近法所求初值问题的值都接近与真值,是有实际意义的,而用函数double计算的值严重失真,没有实际意义。
下面再取n=100时,向前欧拉公式在每个子区间 Ren?
?
004y?
?
(?
k),?
k?
?
xk,xk?
1?
k?
1,2,...,100。
这说明,节点数n越大,向前欧拉公式在每个子区间?
xk,xk?
1?
上的局部截断误差越小,所以函数double计算的值严重失真。
(二)向前欧拉公式的MATLAB主程序2 用向前欧拉公式求解常微分方程初值问题的数值解及其截断误差公式的MATLAB主程序2 输入量:
funfcn是函数f(x,y),x0和y0是初值问题y(x0)?
y0,b是自变量x的最 大值,h是步长。
输出量:
k是自变量x取值的个数,向量X的元素是自变量x的取值,数组Y的元素是利用向前欧拉公式()求出常微分方程初值问题()在X的元素处的数值解,REn是在每个子区间[xx,xk?
1]上的局部截断误差公式,h是步长,其中dy2是 y(x)的2阶导数y?
?
(x)。
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(‘用向前欧拉公式计算dy/dx=f(x,y),y(x0)=y0在[x0,b] 上的数值解’) end k1=1:
n;k=[0,k1];P=[k’,X,Y];symsdy2, REn=*dy2*h_; 例用向前欧拉公式求解初值问题 dyx?
y?
y(0)?
1,取dx3y 182. 第三篇第十章常微分方程求解 n?
10,100,并将计算结果与精确解作比较,写出在每个子区间[xk,xk?
1]上的局部截断误差公式,画出数值解与精确解在区间[0,2]上的图形. 解
(1)建立并保存名为的M文件函数。
functionf=function(x,y) f=y-x/(3*y);
(2)建立并保存名为的M文件函数。
(3)输入程序 >>S1=dsolve(‘Dy=y-x/(3*y)’,’y(0)=1’,’x’)(4)输入程序 >>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*e=1在[0,2]上的精确解’) holdoff jwY1=S2-Y1;xwY1=jwY1./Y1;k1=1:
n1;k=[0,k1];P2=[k’,X1,Y1,S2,jwY1,xwY1] 运行后屏幕显示取n?
10,100时,此初值问题在[0,2]上的自变量x的向量X,X1,利用向前欧拉公式求出的与X,X1对应的数值解Y,Y1及其绝对误差向量jwY,jwY1和相对误差向量xwY,xwY1,每个子区间[xk,xk?
1]上的局部截断误差公式Ren2=1/200*dy2,数值解Y与精确解的图形. (三)自适应向前欧拉公式的MATLAB主程序 用自适应向前欧拉公式求解常微分方程初值问题的数值解
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 第十章 常微分方程组求解 第十 微分 方程组 求解