(精品)第6章微分方程问题的解法

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,第六章,微分方程问题的解法,微分方程的解析解方法,常微分方程问题的数值解法,微分方程问题算法概述,四阶定步长,Runge-Kutta,算法及,MATLAB,实现,一阶微分方程组的数值解,微分方程转换,特殊微分方程的数值解,边值问题的计算机求解,偏微分方程的解,6.1,微分方程的解析解方法,格式：,y=dsolve(f,1, f,2, , f,m,),格式：指明自变量,y=dsolve(f,1, f,2, , f,m,x),f,i,即可以描述微分方程，又可描述初始条件或边界条件。如：,描述微分方程时,描述条件时,例：,syms,t; u=exp(-5*t)*cos(2*t+1)+5;,uu,=5*diff(u,t,2)+4*diff(u,t)+2*u,uu,=,87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+10,syms,t y;, y=dsolve(D4y+10*D3y+35*D2y+50*Dy+24*y=,.,87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+10), y=dsolve(D4y+10*D3y+35*D2y+50*Dy+24*y=,.,87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1) . +10, y(0)=3, Dy(0)=2, D2y(0)=0, D3y(0)=0),分别处理系数，如：, n,d=rat(double(vpa(-445/26*cos(1)-51/13*sin(1)-69/2),ans,=,-8704 185 % rat(),最接近有理数的分数,判断误差：, vpa(-445/26*cos(sym(1)-51/13*sin(1)-69/2+8704/185),ans,=,.114731975864790922564144636e-4,y=dsolve(D4y+10*D3y+35*D2y+50*Dy+24*y=,.,87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1) + . 10,y(0)=1/2,Dy(pi)=1,D2y(2*pi)=0,Dy(2*pi)=1/5);,如果用推导的方法求,Ci,的值，每个系数的解析解至少要写出,10,数行,故可采用有理式近似 的方式表示,., vpa(y,10) %,有理近似值,ans,=,1.196361839*exp(-5.*t)+.4166666667-.4785447354*sin(t)*,cos(t,)*exp(-5.*t)-.4519262218e-1*cos(2.*t)*exp(-5.*t)-2.392723677*cos(t)2*exp(-5.*t)+.2259631109*sin(2.*t)*exp(-5.*t)-473690.0893*exp(-3.*t)+31319.63786*exp(-2.*t)-219.1293619*exp(-1.*t)+442590.9059*exp(-4.*t),例：求解, x,y=dsolve(D2x+2*,Dx,=x+2*y-exp(-t), ,Dy,=4*x+3*y+4*exp(-t),例：,syms,t x, x=,dsolve(Dx,=x*(1-x2),x =, 1/(1+exp(-2*t)*C1)(1/2), -1/(1+exp(-2*t)*C1)(1/2),syms,t x; x=,dsolve(Dx,=x*(1-x2)+1),Warning: Explicit solution could not be found; implicit solution returned., In D:MATLAB6p5toolboxsymbolicdsolve.m at line 292,x =,t-Int(1/(a-a3+1),a=.x)+C1=0,故只有部分非线性微分方程有解析解。,6.2,微分方程问题的数值解法,6.2.1,微分方程问题算法概述,微分方程求解的误差与步长问题：,6.2.2,四阶定步长,Runge-Kutta,算法 及,MATLAB,实现,function ,tout,yout,=rk_4(odefile,tspan,y0),y0,初值列向量,t0=tspan(1);,th,=tspan(2);,if,length(tspan,) t_final=100; x0=0;0;1e-10;,% t_final,为设定的仿真终止时间, t,x=ode45(lorenzeq,0,t_final,x0); plot(t,x), figure; %,打开新图形窗口, plot3(x(:,1),x(:,2),x(:,3);, axis(10 42 -20 20 -20 25); %,根据实际数值手动设置坐标系,可采用,comet3( ),函数绘制动画式的轨迹。, comet3(x(:,1),x(:,2),x(:,3),描述微分方程是常微分方程初值问题数值求解的关键。, f1=inline(-8/3*x(1)+x(2)*x(3); -10*x(2)+10*x(3);,.,-x(1)*x(2)+28*x(2)-x(3),t,x);, t_final=100; x0=0;0;1e-10;, t,x=ode45(f1,0,t_final,x0);, plot(t,x), figure;, plot3(x(:,1),x(:,2),x(:,3); axis(10 42 -20 20 -20 25);,得出完全一致的结果。,6.2.3.3 MATLAB,下带有附加参数的微分方程求解,例：,编写函数,function,xdot,=lorenz1(t,x,flag,beta,rho,sigma),flag,变量是不能省略的,xdot,=-beta*x(1)+x(2)*x(3);,-,rho,*x(2)+rho*x(3);,-x(1)*x(2)+sigma*x(2)-x(3);,求微分方程：, t_final=100; x0=0;0;1e-10;, b2=2; r2=5; s2=20;, t2,x2=ode45(lorenz1,0,t_final,x0,b2,r2,s2);, plot(t2,x2),options,位置为,，,表示不需修改控制选项, figure; plot3(x2(:,1),x2(:,2),x2(:,3); axis(0 72 -20 22 -35 40);,f2=inline(-beta*x(1)+x(2)*x(3); -,rho,*x(2)+rho*x(3);,.,-x(1)*x(2)+sigma*x(2)-x(3), ,t,x,flag,beta,rho,sigma,);,flag,变量是不能省略的,6.2.4,微分方程转换,6.2.4.1,单个高阶常微分方程处理方法,例：,函数描述为：,function y=,vdp_eq(t,x,flag,mu,),y=x(2); -,mu,*(x(1)2-1)*x(2)-x(1);, x0=-0.2; -0.7; t_final=20;,mu,=1; t1,y1=ode45(vdp_eq,0,t_final,x0,mu);,mu,=2; t2,y2=ode45(vdp_eq,0,t_final,x0,mu);, plot(t1,y1,t2,y2,:), figure; plot(y1(:,1),y1(:,2),y2(:,1),y2(:,2),:), x0=2;0; t_final=3000;,mu,=1000; t,y=ode45(vdp_eq,0,t_final,x0,mu);,由于变步长所采用的步长过小，所需时间较长，导致输出的,y,矩阵过大，超出计算机存储空间容量。所以不适合采用,ode45(),来求解，可用刚性方程求解算法,ode15s( ),。,6.2.4.2,高阶常微分方程组的变换方法,例：,描述函数：,function,dx,=,apolloeq(t,x,),mu,=1/82.45; mu1=1-mu;,r1=sqrt(x(1)+mu)2+x(3)2);,r2=sqrt(x(1)-mu1)2+x(3)2);,dx,=x(2);,2*x(4)+x(1)-mu1*(x(1)+mu)/r13-mu*(x(1)-mu1)/r23;,x(4);,-2*x(2)+x(3)-mu1*x(3)/r13-mu*x(3)/r23;,求解：, x0=1.2; 0; 0; -1.04935751;, tic, t,y=ode45(apolloeq,0,20,x0);,toc,elapsed_time =,0.8310, length(t), plot(y(:,1),y(:,3),ans,=,689,得出的轨道不正确，,默认精度,RelTol,设置,得太大，从而导致的,误差传递，可减小该,值。,改变精度：, options=,odeset,;,options.RelTol,=1e-6;, tic, t1,y1=ode45(apolloeq,0,20,x0,options);,toc,elapsed_time =,0.8110, length(t1), plot(y1(:,1),y1(:,3),ans,=,1873, min(diff(t1),ans,=,1.8927e-004, plot(t1(1:end-1),diff(t1),例：, x0=1.2; 0; 0; -1.04935751;, tic, t1,y1=rk_4(apolloeq,0,20,0.01,x0);,toc,elapsed_time =,4.2570, plot(y1(:,1),y1(:,3),%,绘制出轨迹曲线,显而易见，这样求解,是错误的，应该采用,更小的步长。, tic, t2,y2=rk_4(apolloeq,0,20,0.001,x0);,toc,elapsed_time =,124.4990,计算时间过长, plot(y2(:,1),y2(:,3),%,绘制出轨迹曲线,严格说来某些点仍不,满足,10,6,的误差限，,所以求解常微分方程,组时建议采用变步长,算法，而不是定步长,算法。,例：,用,MATLAB,符号工具箱求解，,令 ,syms,x1 x2 x3 x4, ,dx,dy,=solve(dx+2*x4*x1=2*,dy, ,dx,*x4+ 3*x2*dy+x1*x4-x3=5,dx,dy) %,dx,dy,为指定变量,dx,=,-2*(3*x4*x1*x2+x4*x1-x3-5)/(2*x4+3*x2),dy,=,(2*x42*x1-x4*x1+x3+5)/(2*x4+3*x2),对于更复杂的问题来说，手工变换的难度将很大，所以如有可能，可采用计算机去求解有关方程，获得解析解。如不能得到解析解，也需要在描写一阶常微分方程组时列写出式子，得出问题的数值解。,6.3,特殊微分方程的数值解,6.3.1,刚性微分方程的求解,刚性微分方程,一类特殊的常微分方程，其中一些解变化缓慢，另一些变化快，且相差悬殊，这类方程常常称为刚性方程。,MATLAB,采用求解函数,ode15s(),，,该函数的调用格式和,ode45(),完全一致。,t,x=ode15s(Fun,t,0,t,f,x,0,options,p,1,p,2,),例：,计算, h_opt=,odeset,;,h_opt.RelTol,=1e-6;, x0=2;0; t_final=3000;, tic,mu,=1000; t,y=ode15s(vdp_eq,0,t_final,x0,h_opt,mu);,toc,elapsed_time =,2.5240,作图, plot(t,y(:,1); figure; plot(t,y(:,2),y(:,1),曲线变化较平滑, y(:,2),变化在某些点上较快。,例：,定义函数,function,dy,=c7exstf2(t,y),dy,=0.04*(1-y(1)-(1-y(2)*y(1)+0.0001*(1-y(2)2;,-104*y(1)+3000*(1-y(2)2;,方法一, tic,t2,y2=ode45(c7exstf2,0,100,0;1);,toc,elapsed_time =,229.4700, length(t2), plot(t2,y2),ans,=,356941,步长分析：, format long, min(diff(t2), max(diff(t2),ans,=,0.00022220693884 0.00214971787184, plot(t2(1:end-1),diff(t2),方法二，用,ode15s(),代替,ode45(), opt=,odeset,;,opt.RelTol,=1e-6;, tic,t1,y1=ode15s(c7exstf2,0,100,0;1,opt);,toc,elapsed_time =,0.49100000000000, length(t1), plot(t1,y1),ans,=,169,6.3.2,隐式微分方程求解,隐式微分方程为不能转化为显式常微分方程组的方程,例：,编写函数：,function,dx,=c7ximp(t,x),A=sin(x(1) cos(x(2); -cos(x(2) sin(x(1);,B=1-x(1); -x(2);,dx,=inv(A)*B;,求解：, opt=,odeset,;,opt.RelTol,=1e-6;, t,x=ode45(c7ximp,0,10,0; 0,opt); plot(t,x),6.3.3,微分代数方程求解,例：,编写函数,function,dx,=c7eqdae(t,x),dx,=-0.2*x(1)+x(2)*x(3)+0.3*x(1)*x(2);,2*x(1)*x(2)-5*x(2)*x(3)-2*x(2)*x(2);,x(1)+x(2)+x(3)-1;, M=1,0,0; 0,1,0; 0,0,0;, options=,odeset,;, options.Mass=M;,Mass,微分代数方程中,的质量矩阵（控制参数）, x0=0.8; 0.1; 0.1;,t,x=,ode15s,(c7eqdae,0,20,x0,options); plot(t,x),编写函数：,function,dx,=c7eqdae1(t,x),dx,=-0.2*x(1)+x(2)*(1-x(1)-x(2)+0.3*x(1)*x(2);,2*x(1)*x(2)-5*x(2)*(1-x(1)-x(2)-2*x(2)*x(2);, x0=0.8; 0.1;,fDae,=inline(-0.2*x(1)+x(2)*(1-x(1)-x(2)+0.3*x(1)*x(2);,.,2*x(1)*x(2)-5*x(2)*(1-x(1)-x(2)-2*x(2)*x(2),t,x);, t1,x1=ode45(fDae,0,20,x0); plot(t1,x1,t1,1-sum(x1),6.3.3,延迟微分方程求解,sol:,结构体数据，,sol.x:,时间向量,t, sol.y:,状态向量。,例：,编写函数：,function,dx,=c7exdde(t,x,z),xlag1=z(:,1); %,第一列表示提取,xlag2=z(:,2);,dx,=1-3*x(1)-xlag1(2)-0.2*xlag2(1)3-xlag2(1);,x(3); 4*x(1)-2*x(2)-3*x(3);,历史数据函数：,function S=c7exhist(t),S=zeros(3,1);,求解：, lags=1 0.5;,tx,=dde23(c7exdde,lags,zeros(3,1),0,10);, plot(tx.x,tx.y(2,:),与,ode45(),等返回的,x,矩阵不一样，它是按行排列的。,6.4,边值问题的计算机求解,6.4.1,边值问题的打靶算法,数学方法描述：,以二阶方程为例,编写函数： 线性的,function t,y=shooting(f1,f2,tspan,x0f,varargin),t0=tspan(1);,tfinal,=tspan(2);,ga,=x0f(1);,gb,=x0f(2);,t,y1=ode45(f1,tspan,1;0,varargin);,t,y2=ode45(f1,tspan,0;1,varargin);,t,yp,=ode45(f2,tspan,0;0,varargin);,m=(,gb-ga,*y1(end,1)-yp(end,1)/y2(end,1);,t,y=ode45(f2,tspan,ga;m,varargin);,例：,编写函数：,function,xdot,=c7fun1(t,x),xdot,=x(2); -2*x(1)+3*x(2);,function,xdot,=c7fun2(t,x),xdot,=x(2); t-2*x(1)+3*x(2);, t,y=shooting(c7fun1, ,c7fun2,0,1,1;2); plot(t,y),原方程的解析解为,解的检验, y0=(exp(2)-3)*exp(t)+(3-exp(1)*exp(2*t)/(4*exp(1)*(exp(1)-1)+3/4+t/2;, norm(y(:,1)-y0) %,整个解函数检验,ans,=,4.4790e-008, norm(y(end,1)-2) %,终点条件检验,ans,=,2.2620e-008,非线性方程边值问题的打靶算法：,用,Newton,迭代法处理,编写函数：,function t,y=nlbound(funcn,funcv,tspan,x0f,tol,varargin),t0=tspan(1);tfinal=tspan(2);,ga,=x0f(1);,gb,=x0f(2); m=1; m0=0;,while (norm(m-m0),tol,), m0=m;,t,v=ode45(funcv,tspan,ga;m;0;1,varargin);,m=m0-(v(end,1)-gb)/(v(end,3);,end,t,y=ode45(funcn,tspan,ga;m,varargin);,例：,编写两个函数：,function,xdot,=c7fun3(t,x),xdot,=x(2); 2*x(1)*x(2); x(4); 2*x(2)*x(3)+2*x(1)*x(4);,function,xdot,=c7fun4(t,x),xdot,=x(2); 2*x(1)*x(2);, t,y=nlbound(c7fun4,c7fun3,0,pi/2,-1,1,1e-8);, plot(t,y); set(gca,xlim,0,pi/2);,精确解：,检验：, y0=tan(t-pi/4);, norm(y(:,1)-y0),ans,=,1.6629e-005, norm(y(end,1)-1),ans,=,5.2815e-006,6.4.2,线性微分方程的有限差分算法,把等式左边用差商表示。,编写函数：,function x,y=fdiff(funcs,tspan,x0f,n),t0=tspan(1);tfinal=tspan(2);,ga,=x0f(1);,gb,=x0f(2);,h=(tfinal-t0)/n;,for i=1:n, x(i)=t0+h*(i-1); end,x0=x(1:n-1); t=-2+h2*feval(funcs,x0,2);,tmp,=feval(funcs,x0,1);,v=1+h*tmp/2; w=1-h*tmp/2; b=h2*feval(funcs,x0,3);,b(1)=b(1)-w(1)*,ga,; b(n-1)=b(n-1)-v(n-1)*,gb,; b=b; A=,diag(t,);,for i=1:n-2, A(i,i+1)=v(i); A(i+1,i)=w(i+1); end,y=inv(A)*b; x=x,tfinal,; y=,ga,; y;,gb,;,例：,编写函数：,function y=c7fun5(x,key),switch key,case 1, y=1+x;,case 2, y=1-x;,otherwise, y=1+x.2;,end, t,y=fdiff(c7fun5,0,1,1,4,50); plot(t,y),6.5,偏微分方程求解入门,6.5.1,偏微分方程组求解,函数描述：,边界条件的函数描述：,初值条件的函数描述：,u0=,pdeic(x,),例：,函数描述：,function c,f,s=c7mpde(x,t,u,du),c=1;1; y=u(1)-u(2); F=exp(5.73*y)-exp(-11.46*y); s=F*-1; 1;,f=0.024*du(1); 0.17*du(2);,描述边界条件的函数,function ,pa,qa,pb,qb,=c7mpbc(xa,ua,xb,ub,t),pa=0; ua(2);,qa,=1;0;,pb,=ub(1)-1; 0;,qb,=0;1;,描述初值：,function u0=c7mpic(x),u0=1; 0;,求解：, x=0:.05:1; t=0:0.05:2; m=0;, sol=pdepe(m,c7mpde,c7mpic,c7mpbc,x,t);, surf(x,t,sol(:,:,1),，,figure; surf(x,t,sol(:,:,2),6.5.2,二阶偏微分方程的数学描述,椭圆型偏微分方程：,抛物线型偏微分方程,：,双曲型偏微分方程：,特征值型偏微分方程：,6.5.3,偏微分方程的求解界面应用简介,6.5.3.1,偏微分方程求解程序概述,启动偏微分方程求解界面,在,MATLAB,下键入,pdetool,该界面分为四个部分,菜单系统,工具栏,集合编辑,求解区域,6.5.3.2,偏微分方程求解区域绘制,1,）用工具栏中的椭圆、矩形等绘制一些区域。,2,）在集合编辑栏中修改其内容。,如（,R1,E1,E2,）,E3,3,）,单击工具栏中 按纽可得求解边界。,4,）选择,Boundary-Remove All,Subdomain,Borders,菜单项，消除相邻区域中间的分隔线。,5,）单击 按纽可将求解区域用三角形划分成网格。可用 按纽加密。,6.5.3.3,偏微分方程边界条件描述,选择,Boundary-Specify Boundary Conditions,菜单,6.5.3.4,偏微分方程求解举例,例：,求解：,1,）绘制求解区域。,2,）描述边界条件,（,Boundary-Specify Boundary Conditions,）,。,3,）选择偏微分方程的类型。单击工具栏中的,PDE,图标，在打开的新窗口选择,Hyperbolic,选项，输入参数,c,a,f,d.,4,）,求解。单击工具栏中的等号按钮。,显示：,1,）图形颜色表示,t=0,时,u(x,y),的函数值。,2,）单击工具栏中的三维图标将打开一新的对话框，若再选择,Contour,可绘制等值线，若选择,Arrows,选项可绘制引力线。若单独选择,Height,（,3d-plot),，,则在另一窗口绘制出三维图形。,3,）可在单击三维图标打开的新对话框中，对,Property,栏目的各个项目重新选择。,4,）可修改微分方程的边界条件，重新求解。,动画：,1,）,Solve-Parameters,对话框时间向量改为,0:0.1:2,。,2,）三维图标打开的对话框中选择,Animation,选项，单击,Options,按纽设置播放速度。,Plot-Export Movie,菜单。,6.5.3.5,函数参数的偏微分方程求解,例：,(,椭圆型）,求解：,1,）求解区域不变。,2,）描述边界条件，,u=0,。,3,）选择偏微分方程的类型。单击工具栏中的,PDE,图标，在打开的新窗口选择,Elliptic,选项，输入参数,c=1./sqrt(1+ux.2+uy.2), a=x.2+y.2 , f=exp(-x.2-y.2).,4),再,打开,Solve-Parameters,对话框，选定,Use nonlinear solve,属性（该属性只适于椭圆性偏微分方程）,5,）,求解。单击工具栏中的等号按钮。,