微分方程问题的解法

微分方程问题的解法第一页，共八十五页，2022年，8月28日6.1微分方程的解析解方法格式：

y=dsolve(f1,f2,…,fm)格式：指明自变量

y=dsolve(f1,f2,…,fm,’x’)fi即可以描述微分方程，又可描述初始条件或边界条件。如：描述微分方程时描述条件时第二页，共八十五页，2022年，8月28日例：>>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)+10>>symsty;>>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'])第三页，共八十五页，2022年，8月28日>>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')第四页，共八十五页，2022年，8月28日分别处理系数，如：>>[n,d]=rat(double(vpa(-445/26*cos(1)-51/13*sin(1)-69/2)))]ans=-8704185%rat()最接近有理数的分数判断误差：>>vpa(-445/26*cos(sym(1))-51/13*sin(1)-69/2+8704/185)ans=第五页，共八十五页，2022年，8月28日>>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)第六页，共八十五页，2022年，8月28日例：求解>>[x,y]=dsolve('D2x+2*Dx=x+2*y-exp(-t)',…'Dy=4*x+3*y+4*exp(-t)')第七页，共八十五页，2022年，8月28日例：>>symstx>>x=dsolve('Dx=x*(1-x^2)')x=[1/(1+exp(-2*t)*C1)^(1/2)][-1/(1+exp(-2*t)*C1)^(1/2)]>>symstx;x=dsolve('Dx=x*(1-x^2)+1')Warning:Explicitsolutioncouldnotbefound;implicitsolutionreturned.>InD:\MATLAB6p5\toolbox\symbolic\dsolve.matline292x=t-Int(1/(a-a^3+1),a=``..x)+C1=0故只有部分非线性微分方程有解析解。第八页，共八十五页，2022年，8月28日6.2微分方程问题的数值解法

6.2.1微分方程问题算法概述第九页，共八十五页，2022年，8月28日微分方程求解的误差与步长问题：第十页，共八十五页，2022年，8月28日第十一页，共八十五页，2022年，8月28日第十二页，共八十五页，2022年，8月28日6.2.2四阶定步长Runge-Kutta算法

及MATLAB实现第十三页，共八十五页，2022年，8月28日function[tout,yout]=rk_4(odefile,tspan,y0)％y0初值列向量

t0=tspan(1);th=tspan(2);iflength(tspan)<=3,h=tspan(3);%tspan=[t0,th,h]else,h=tspan(2)-tspan(1);th=tspan(end);end％等间距数组tout=[t0:h:th]';yout=[];fort=tout'k1=h*eval([odefile‘(t,y0)’]);%odefile是一个字符串变量，为表示微分方程f()的文件名。

k2=h*eval([odefile'(t+h/2,y0+0.5*k1)']);k3=h*eval([odefile'(t+h/2,y0+0.5*k2)']);k4=h*eval([odefile'(t+h,y0+k3)']);y0=y0+(k1+2*k2+2*k3+k4)/6;yout=[yout;y0'];end％由效果看，该算法不是一个较好的方法。第十四页，共八十五页，2022年，8月28日6.2.3一阶微分方程组的数值解

6.2.3.1四阶五级Runge-Kutta-Felhberg算法通过误差向量调节步长，此为自动变步长方法。

四阶五级RKF算法有参量系数表。第十五页，共八十五页，2022年，8月28日6.2.3.2基于MATLAB的微分方程

求解函数

格式1：直接求解

[t,x]=ode45(Fun,[t0,tf],x0)

格式2：带有控制参数

[t,x]=ode45(Fun,[t0,tf],x0，options)

格式3：带有附加参数

[t,x]=ode45(Fun,[t0,tf],x0,options,p1,p2,…)

[t0,tf]求解区间，x0初值问题的初始状态变量。

第十六页，共八十五页，2022年，8月28日描述需要求解的微分方程组：不需附加变量的格式

functionxd=funname(t,x)

可以使用附加变量

functionxd=funname(t,x,flag,p1,p2,…)

％t是时间变量或自变量（必须给），x为状态向量，

xd为返回状态向量的导数。flag用来控制求解过程,指定初值，即使初值不用指定，也必须有该变量占位。修改变量：options唯一结构体变量，用odeset()修改。

options=odeset(‘RelTol’,1e-7);options=odeset;options.RelTol=1e-7;第十七页，共八十五页，2022年，8月28日例：自变函数

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)];第十八页，共八十五页，2022年，8月28日>>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]);%根据实际数值手动设置坐标系第十九页，共八十五页，2022年，8月28日可采用comet3()函数绘制动画式的轨迹。>>comet3(x(:,1),x(:,2),x(:,3))第二十页，共八十五页，2022年，8月28日描述微分方程是常微分方程初值问题数值求解的关键。>>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]);得出完全一致的结果。第二十一页，共八十五页，2022年，8月28日6.2.3.3MATLAB下带有附加参数的微分方程求解例：第二十二页，共八十五页，2022年，8月28日编写函数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]);第二十三页，共八十五页，2022年，8月28日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变量是不能省略的第二十四页，共八十五页，2022年，8月28日6.2.4微分方程转换

6.2.4.1单个高阶常微分方程处理方法第二十五页，共八十五页，2022年，8月28日第二十六页，共八十五页，2022年，8月28日例：函数描述为：

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),':')第二十七页，共八十五页，2022年，8月28日>>x0=[2;0];t_final=3000;>>mu=1000;[t,y]=ode45('vdp_eq',[0,t_final],x0,[],mu);

由于变步长所采用的步长过小，所需时间较长，导致输出的y矩阵过大，超出计算机存储空间容量。所以不适合采用ode45()来求解，可用刚性方程求解算法ode15s()。第二十八页，共八十五页，2022年，8月28日6.2.4.2高阶常微分方程组的变换方法第二十九页，共八十五页，2022年，8月28日例：第三十页，共八十五页，2022年，8月28日第三十一页，共八十五页，2022年，8月28日描述函数：

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)/r1^3-mu*(x(1)-mu1)/r2^3;x(4);-2*x(2)+x(3)-mu1*x(3)/r1^3-mu*x(3)/r2^3];第三十二页，共八十五页，2022年，8月28日求解：>>x0=[1.2;0;0;-1.04935751];>>tic,[t,y]=ode45('apolloeq',[0,20],x0);tocelapsed_time=0.8310>>length(t),>>plot(y(:,1),y(:,3))ans=689得出的轨道不正确，默认精度RelTol设置得太大，从而导致的误差传递，可减小该值。第三十三页，共八十五页，2022年，8月28日改变精度：>>options=odeset;options.RelTol=1e-6;>>tic,[t1,y1]=ode45('apolloeq',[0,20],x0,options);tocelapsed_time=0.8110>>length(t1),>>plot(y1(:,1),y1(:,3)),ans=1873第三十四页，共八十五页，2022年，8月28日>>min(diff(t1))ans=1.8927e-004>>plot(t1(1:end-1),…diff(t1))第三十五页，共八十五页，2022年，8月28日例：第三十六页，共八十五页，2022年，8月28日>>x0=[1.2;0;0;-1.04935751];>>tic,[t1,y1]=rk_4('apolloeq',[0,20,0.01],x0);tocelapsed_time=4.2570>>plot(y1(:,1),y1(:,3))%绘制出轨迹曲线显而易见，这样求解是错误的，应该采用更小的步长。第三十七页，共八十五页，2022年，8月28日>>tic,[t2,y2]=rk_4('apolloeq',[0,20,0.001],x0);tocelapsed_time=124.4990％计算时间过长>>plot(y2(:,1),y2(:,3))%绘制出轨迹曲线严格说来某些点仍不满足10－6的误差限，所以求解常微分方程组时建议采用变步长算法，而不是定步长算法。第三十八页，共八十五页，2022年，8月28日例：第三十九页，共八十五页，2022年，8月28日用MATLAB符号工具箱求解，令％>>symsx1x2x3x4>>[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*x4^2*x1-x4*x1+x3+5)/(2*x4+3*x2)

对于更复杂的问题来说，手工变换的难度将很大，所以如有可能，可采用计算机去求解有关方程，获得解析解。如不能得到解析解，也需要在描写一阶常微分方程组时列写出式子，得出问题的数值解。第四十页，共八十五页，2022年，8月28日6.3特殊微分方程的数值解

6.3.1刚性微分方程的求解刚性微分方程一类特殊的常微分方程，其中一些解变化缓慢，另一些变化快，且相差悬殊，这类方程常常称为刚性方程。MATLAB采用求解函数ode15s()，该函数的调用格式和ode45()完全一致。[t,x]=ode15s(Fun,[t0,tf],x0,options,p1,p2,…)

第四十一页，共八十五页，2022年，8月28日例：％计算>>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);tocelapsed_time=2.5240第四十二页，共八十五页，2022年，8月28日％作图>>plot(t,y(:,1));figure;plot(t,y(:,2))y(:,1)曲线变化较平滑,y(:,2)变化在某些点上较快。第四十三页，共八十五页，2022年，8月28日例：定义函数functiondy=c7exstf2(t,y)dy=[0.04*(1-y(1))-(1-y(2))*y(1)+0.0001*(1-y(2))^2;-10^4*y(1)+3000*(1-y(2))^2];第四十四页，共八十五页，2022年，8月28日方法一>>tic,[t2,y2]=ode45('c7exstf2',[0,100],[0;1]);tocelapsed_time=229.4700>>length(t2),plot(t2,y2)ans=356941第四十五页，共八十五页，2022年，8月28日步长分析：>>formatlong,[min(diff(t2)),max(diff(t2))]ans=>>plot(t2(1:end-1),diff(t2))第四十六页，共八十五页，2022年，8月28日方法二，用ode15s()代替ode45()>>opt=odeset;opt.RelTol=1e-6;>>tic,[t1,y1]=ode15s('c7exstf2',[0,100],[0;1],opt);tocelapsed_time=>>length(t1),>>plot(t1,y1)ans=169第四十七页，共八十五页，2022年，8月28日6.3.2隐式微分方程求解隐式微分方程为不能转化为显式常微分方程组的方程例：第四十八页，共八十五页，2022年，8月28日编写函数：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)第四十九页，共八十五页，2022年，8月28日6.3.3微分代数方程求解例：第五十页，共八十五页，2022年，8月28日编写函数

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)第五十一页，共八十五页，2022年，8月28日编写函数：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)];第五十二页，共八十五页，2022年，8月28日>>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'))第五十三页，共八十五页，2022年，8月28日延迟微分方程求解sol:结构体数据，sol.x:时间向量t,sol.y:状态向量。

第五十四页，共八十五页，2022年，8月28日例：第五十五页，共八十五页，2022年，8月28日编写函数：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)];历史数据函数：functionS=c7exhist(t)S=zeros(3,1);第五十六页，共八十五页，2022年，8月28日求解：>>lags=[10.5];tx=dde23('c7exdde',lags,zeros(3,1),[0,10]);>>plot(tx.x,tx.y(2,:))％与ode45()等返回的x矩阵不一样，它是按行排列的。第五十七页，共八十五页，2022年，8月28日6.4边值问题的计算机求解第五十八页，共八十五页，2022年，8月28日6.4.1边值问题的打靶算法数学方法描述：以二阶方程为例

第五十九页，共八十五页，2022年，8月28日编写函数：线性的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);第六十页，共八十五页，2022年，8月28日例：编写函数：functionxdot=c7fun1(t,x)xdot=[x(2);-2*x(1)+3*x(2)];functionxdot=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)第六十一页，共八十五页，2022年，8月28日原方程的解析解为解的检验>>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第六十二页，共八十五页，2022年，8月28日非线性方程边值问题的打靶算法：

用Newton迭代法处理第六十三页，共八十五页，2022年，8月28日编写函数：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);第六十四页，共八十五页，2022年，8月28日例：编写两个函数：functionxdot=c7fun3(t,x)xdot=[x(2);2*x(1)*x(2);x(4);2*x(2)*x(3)+2*x(1)*x(4)];functionxdot=c7fun4(t,x)xdot=[x(2);2*x(1)*x(2)];第六十五页，共八十五页，2022年，8月28日>>[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第六十六页，共八十五页，2022年，8月28日6.4.2线性微分方程的有限差分算法把等式左边用差商表示。第六十七页，共八十五页，2022年，8月28日第六十八页，共八十五页，2022年，8月28日编写函数：function[x,y]=fdiff(funcs,tspan,x0f,n)t0=tspan(1);tfinal=tspan(2);ga=x0f(1);gb=x0f(2);h=(tfinal-t0)/n;fori=1:n,x(i)=t0+h*(i-1);end,x0=x(1:n-1);t=-2+h^2*feval(funcs,x0,2);tmp=feval(funcs,x0,1);v=1+h*tmp/2;w=1-h*tmp/2;b=h^2*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);fori=1:n-2,A(i,i+1)=v(i);A(i+1,i)=w(i+1);endy=inv(A)*b;x=[xtfinal];y=[ga;y;gb]';第六十九页，共八十五页，2022年，8月28日例：编写函数：functiony=c7fun5(x,key)switchkeycase1,y=1+x;case2,y=1-x;otherwise,y=1+x.^2;end>>[t,y]=fdiff('c7fun5',[0,1],[1,4],50);plot(t,y)第七十页，共八十五页，2022年，8月28日6.5偏微分方程求解入门

6.5.1偏微分方程组求解函数描述：第七十一页，共八十五页，2022年，8月28日边界条件的函数描述：初值条件的函数描述：

u0=pdeic(x)第七十二页，共八十五页，2022年，8月28日例：第七十三页，共八十五页，2022年，8月28日函数描述：

第七十四页，共八十五页，2022年，8月28日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)];

第七十五页，共八十五页，2022年，8月28日描述边界条件的函数function[pa,qa,pb,qb]=c7mpbc(xa,ua,xb,ub,t)pa=[0;ua(2)];qa=[1;0];