MATLAB-微分方程问题解法14

微分方程（组）数值解,常微分方程（组）数值解偏微分方程（组）数值解,1.微分方程的解析解方法,格式：y=dsolve(f1,f2,fm)格式：指明自变量y=dsolve(f1,f2,fm,x)fi即可以描述微分方程，又可描述初始条件或边界条件。如：描述微分方程时描述条件时,例1：symst;u=exp(-5*t)*cos(2*t+1)+5;uu=5*diff(u,t,2)+4*diff(u,t)+2*uuu=87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+10symsty;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=-8704185%rat()最接近有理数的分数判断误差：vpa(-445/26*cos(1)-51/13*sin(1)-69/2+8704/185)ans=.11473197588429684401489794254303e-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=-.6596153846*exp(-5.*t)*cos(2.*t+1.)-1.051923077*exp(-5.*t)*sin(2.*t+1.)+.4166666667+31319.63806*exp(-2.*t)-219.1293625*exp(-1.*t)-473690.0926*exp(-3.*t)+442590.9086*exp(-4.*t),例1：求解x,y=dsolve(D2x+2*Dx=x+2*y-exp(-t),Dy=4*x+3*y+4*exp(-t),例2：symstxx=dsolve(Dx=x*(1-x2)x=1/(1+exp(-2*t)*C1)(1/2)-1/(1+exp(-2*t)*C1)(1/2)symstx;x=dsolve(Dx=x*(1-x2)+1)Warning:Explicitsolutioncouldnotbefound;implicitsolutionreturned.Indsolveat315x=t+Int(-1/(_a-_a3+1),_a=.x)+C1=0故只有部分非线性微分方程有解析解。,2.微分方程问题的数值解法2.1微分方程问题算法概述,常微分方程数值解/*NumericalMethodsforOrdinaryDifferentialEquations*/,待求解的问题：一阶常微分方程的初值问题/*Initial-ValueProblem*/:,解的存在唯一性（“常微分方程”理论）：只要f(x,y)在a,bR1上连续，且关于y满足Lipschitz条件，即存在与x,y无关的常数L使对任意定义在a,b上的y1(x)和y2(x)都成立，则上述IVP存在唯一解。,解析解法：（常微分方程理论）只能求解极少一类常微分方程；实际中给定的问题不一定是解析表达式，而是函数表，无法用解析解法。,如何求解,步进式：根据已知的或已求出的节点上的函数值计算当前节点上的函数值，一步一步向前推进。因此只需建立由已知的或已求出的节点上的函数值求当前节点函数值的递推公式即可。,-EulersMethod,欧拉方法/*EulersMethod*/,EulersMethod,向前差商代替导数,说明,EulersMethod,截断误差:实际上，y(xn)yn，yn也有误差，它对yn+1的误差也有影响，见下图。但这里不考虑此误差的影响，仅考虑方法或公式本身带来的误差，因此称为方法误差或截断误差。局部截断误差的分析：由于假设yn=y(xn)，即yn准确，因此分析局部截断误差时将y(xn+1)和yn+1都用点xn上的信息来表示，工具：Taylor展开。,欧拉法的局部截断误差：,Rn+1的主项/*leadingterm*/,EulersMethod,欧拉法具有1阶精度。,前面介绍了差商的概念和性质，各阶差商可以近似各阶导数，具有不同的精度，差商近似导数,EulersMethod,EulersMethod,欧拉公式的改进：,隐式欧拉法或后退欧拉法/*implicitEulermethodorbackwardEulermethod*/,隐式或后退欧拉公式,隐式欧拉法的局部截断误差：,即隐式欧拉公式具有1阶精度。,EulersMethod,比较欧拉显式公式和隐式公式及其局部截断误差,显式公式,隐式公式,EulersMethod,若将这两种方法进行算术平均，即可消除误差的主要部分/*leadingterm*/而获得更高的精度,称为梯形法,EulersMethod,梯形公式/*trapezoidformula*/,显、隐式两种算法的平均,注：的确有局部截断误差，即梯形公式具有2阶精度，比欧拉方法有了进步。但注意到该公式是隐式公式，计算时不得不用到迭代法，其迭代收敛性与欧拉公式相似。,中点欧拉公式/*midpointformula*/,假设，则可以导出即中点公式也具有2阶精度，且是显式的。,需要2个初值y0和y1来启动递推过程，这样的算法称为双步法/*double-stepmethod*/，而前面的三种算法都是单步法/*single-stepmethod*/。,EulersMethod,预测-校正-改进系统中点法具有二阶精度，且是显式的，与梯形公式精度相匹配，用中点公式作预测，梯形公式作校正，得到如下预测校正系统：,校正误差约为预测误差的1/4,EulersMethod,summary,龙格-库塔法/*Runge-KuttaMethod*/,Runge-KuttaMethod,单步递推法的基本思想是从(xn,yn)点出发，以某一斜率沿直线达到(xn+1,yn+1)点。欧拉法及其各种变形所能达到的最高精度为2阶。,由微分中值定理，有,k*称为区间xn,xn+1上的平均斜率，只要知道平均斜率，就可计算y(xn+1).因此只要对平均斜率提供一种近似算法，则由(4)式可导出一种相应的求解公式。,(4),例,Runge-KuttaMethod,由此看出，改进的欧拉公式用xn与xn+1两个节点的斜率的算术平均作为平均斜率，xn+1点的斜率通过已知信息yn来预测。,考察改进的欧拉法，可以将其改写为：,斜率一定取K1，K2的平均值吗？,步长一定是h吗？即第二个节点一定是xn+1吗？,Runge-KuttaMethod,Runge-KuttaMethod,首先希望能确定系数1、2、p，使得到的算法格式有2阶精度，即在的前提假设下，使得,Step1:将K2在(xi,yi)点作Taylor展开,Step2:将K2代入第1式，得到,2阶RungeKuttaMethod,Runge-KuttaMethod,Step3:将yi+1与y(xi+1)在xi点的泰勒展开作比较,要求，则必须有：,这里有个未知数，个方程。,3,2,存在无穷多个解。所有满足上式的格式统称为2阶龙格-库塔格式。,注意到，就是改进的欧拉法；,p=1/2,1=0,2=1,变形欧拉公式。,Q:为获得更高的精度，应该如何进一步推广？改进的Euler公式推广为二阶Runge-Kutta公式带来这样的启示：若在xn,xn+1上多预测几个点的斜率值，然后将它们的算术平均作为平均斜率，则有可能构造出具有更高精度的计算公式。-Runge-Kutta方法的基本思想。,注：二阶Runge-Kutta公式用多算一次函数值f的办法避开了二阶Taylor级数法所要计算的f的导数。在这种意义上，可以说Runge-Kutta方法实质上是Taylor级数法的变形。,Runge-KuttaMethod,其中i(i=1,m)，i(i=2,m)和ij(i=2,m;j=1,i1)均为待定系数，确定这些系数的步骤与前面相似。,Runge-KuttaMethod,高阶RungeKuttaMethod,Gill公式：4阶经典龙格-库塔公式的一种改进,Runge-KuttaMethod,最常用为4阶经典龙格-库塔法/*ClassicalRunge-KuttaMethod*/：,Runge-KuttaMethod,由于龙格-库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好采用低阶算法而将步长h取小。,Matlab编程,例3：用Euler法解初值问题,解：,取步长h=(2-0)/n=2/n，得差分方程,a=0;%ab为积分区间b=2;n=5;%中间点数h=(b-a)/(n-1);%步长x=linspace(a,b,n);y=zeros(size(x);y(1)=1;fork=1:n-1y(k+1)=y(k)+h*(y(k)+x(k)/(y(k).2);endplot(x,y,o-),Matlab编程,解析解：,Matlab编程,为了减小误差，可采用以下方法：,让步长h取得更小一些；,改用具有较高精度的数值方法：,Runge-Kutta(龙格-库塔)方法,步长减小，精度会提高。,Matlab编程,四阶R-K方法,其中,Matlab编程,a=0;b=2;n=50;h=(b-a)/(n-1);x=linspace(a,b,n);y=zeros(size(x);y(1)=1;f=(u,v)u+v/u2;fork=1:n-1L1=f(y(k),x(k);L2=f(y(k)+h*L1/2,x(k)+h/2);L3=f(y(k)+h*L2/2,x(k)+h/2);L4=f(y(k)+h*L3,x(k)+h);y(k+1)=y(k)+h*(L1+2*L2+2*L3+L4)/6;end,Matlab编程,Matlab编程,用Maltab自带函数解初值问题,求解析解：dsolve,求数值解：ode45、ode23、ode113、ode23t、ode15s、ode23s、ode23tb,Matlab编程,T,Y=solver(odefun,tspan,y0),其中y0为初值条件，tspan为求解区间；Matlab在数值求解时自动对求解区间进行分割，T(向量)中返回的是分割点的值(自变量)，Y(向量)中返回的是解函数在这些分割点上的函数值。solver为Matlab的ODE求解器（可以是ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb）,没有一种算法可以有效地解决所有的ODE问题，因此MATLAB提供了多种ODE求解器，对于不同的ODE，可以调用不同的求解器。,Matlab编程,Matlab编程,例4：自变函数functionxdot=lorenzeq(t,x)xdot=-8/3*x(1)+x(2)*x(3);-10*x(2)+10*x(3);-x(1)*x(2)+28*x(2)-x(3);,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(1042-2020-2025);%根据实际数值手动设置坐标系,可采用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(1042-2020-2025);得出完全一致的结果。,2.3.3MATLAB下带有附加参数的微分方程求解,例5：,编写函数functionxdot=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(072-2022-3540);,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变量是不能省略的,2.4微分方程转换2.4.1单个高阶常微分方程处理方法,Vanderpol方程概述,电子科技大学软件学院,BalthazarvanderPol(1889-1959),例6：函数描述为：functiony=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()。,2.4.2高阶常微分方程组的变换方法,例7：,描述函数：functiondx=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);tocElapsedtimeis0.211905seconds.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);tocElapsedtimeis0.232374seconds.length(t1),plot(y1(:,1),y1(:,3),ans=1873,3.特殊微分方程的数值解3.1刚性微分方程的求解,刚性微分方程一类特殊的常微分方程，其中一些解变化缓慢，另一些变化快，且相差悬殊，这类方程常常称为刚性方程。MATLAB采用求解函数ode15s()，该函数的调用格式和ode45()完全一致。t,x=ode15s(Fun,t0,tf,x0,options,p1,p2,),3.2隐式微分方程求解,隐式微分方程为不能转化为显式常微分方程组的方程例8：,编写函数：functiondx=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),3.3微分代数方程求解,例9：,编写函数functiondx=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),编写函数：functiondx=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),3.3延迟微分方程求解,sol:结构体数据，sol.x:时间向量t,sol.y:状态向量。,例10：,编写函数：functiondx=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);,求解：lags=10.5;tx=dde23(c7exdde,lags,zeros(3,1),0,10);plot(tx.x,tx.y(2,:)与ode45()等返回的x矩阵不一样，它是按行排列的。,4.偏微分方程求解入门4.1偏微分方程组求解,函数描述：,边界条件的函数描述：初值条件的函数描述：u0=pdeic(x),例11：,函数描述：,functionc,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);,描述边界条件的函数functionpa,qa,pb,qb=c7mpbc(xa,ua,xb,ub,t)pa=0;ua(2);qa=1;0;pb=ub(1)-1;0;qb=0;1;,描述初值：functionu0=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),4.2二阶偏微分方程的数学描述,椭圆型偏微分方程：,抛物线型偏微分方程：双曲型偏微分方程：特征值型偏微分方程：,例12求解Poisson方程：,%（2）产生初始的三角形网格p,e,t=initmesh(g);%（3）迭代直至得到误差允许范围内的合格解error=;err=1;whileerr0.01,p,e,t=refinemesh(g,p,e,t);u=assempde(b,p,e,t,c,a,f);%求得数值解exact=(1-p(1,:).2-p(2,:).2)/4;err=norm(u-exact,inf);error=errorerr;end%结果显示subplot(1,3,1),pdemesh(p,e,t);subplot(1,3,2),pdesurf(p,t,u)subplot(1,3,3),pdesurf(p,t,u-exact),例13考虑最小表面问题：,例14：,%（2）产生初始的三角形网格p,e,t=initmesh(g);%(3)定义初始条件u0=zeros(size(p,2),1);ix=find(sqrt(p(1,:).2+p(2,:).2)0.4);u0(ix)=1%(4)在时间段为0到0.1的20个点上求解nframe=20;tlist=linspace(0,0.1,nframe);u1=parabolic(u0,tlist,b,p,e,t,c,a,f,d);%(5)动画图示结果forj=1:nframepdesurf(p,t,u1(:,j);mv(j)=getframe;endmovie(mv,10),例15：,%(3)定义初始条件x=p(1,:);y=p(2,:);u0=atan(cos(pi*x);ut0=3*sin(pi*x).*exp(cos(pi*y);%(4)在时间段为0到5的31个点上求解n=31;tlist=linspace(0,5,n);uu=hyperbolic(u0,ut0,tlist,b,p,e,t,c,a,f,d);%(5)动画图示结果forj=1:npdesurf(p,t,uu(:,j);mv(j)=getframe;endmovie(mv,10),4.3偏微分方程的求解界面应用简介4.3.1偏微分方程求解程序概述,启动偏微分方程求解界面在MATLAB下键入pdetool该界面分为四个部分菜单系统工具栏集合编辑求解区域,4.3.2偏微分方程求解区域绘制,1）用工具栏中的椭圆、矩形等绘制一些区域。2）在集合编辑栏中修改其内容。如（R1E1E2）E33）单击工具栏中按纽可得求解边界。4）选择Boundary-RemoveAllSubdomainBorders菜单项，消除相邻区域中间的分隔线。5）单击按纽可将求解区域用三角形划分成网格。可用按纽加密。,4.3.3偏微分方程边界条件描述,选择Boundary-SpecifyBoundaryConditions菜单,4.3.4偏微分方程求解举例,例：求解：1）绘制求解区域。2）描述边界条件（Boundary-SpecifyBoundaryConditions）。3）选择偏微分方程的类型。单击工具栏中的PDE图标，在打开的新窗口选择Hyperbolic选项