微分方程問題的解法

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

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

y=dsolve(f1,f2,…,fm,’x’)fi即可以描述微分方程，又可描述初始條件或邊界條件。如：描述微分方程時(shí)描述條件時(shí)第二頁，共八十五頁，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日分別處理系數(shù)，如：>>[n,d]=rat(double(vpa(-445/26*cos(1)-51/13*sin(1)-69/2)))]ans=-8704185%rat()最接近有理數(shù)的分?jǐn)?shù)判斷誤差：>>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');

如果用推導(dǎo)的方法求Ci的值，每個(gè)系數(shù)的解析解至少要寫出10數(shù)行,故可采用有理式近似的方式表示.>>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微分方程問題的數(shù)值解法

6.2.1微分方程問題算法概述第九頁，共八十五頁，2022年，8月28日微分方程求解的誤差與步長(zhǎng)問題：第十頁，共八十五頁，2022年，8月28日第十一頁，共八十五頁，2022年，8月28日第十二頁，共八十五頁，2022年，8月28日6.2.2四階定步長(zhǎng)Runge-Kutta算法

及MATLAB實(shí)現(xiàn)第十三頁，共八十五頁，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％等間距數(shù)組tout=[t0:h:th]';yout=[];fort=tout'k1=h*eval([odefile‘(t,y0)’]);%odefile是一個(gè)字符串變量，為表示微分方程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％由效果看，該算法不是一個(gè)較好的方法。第十四頁，共八十五頁，2022年，8月28日6.2.3一階微分方程組的數(shù)值解

6.2.3.1四階五級(jí)Runge-Kutta-Felhberg算法通過誤差向量調(diào)節(jié)步長(zhǎng)，此為自動(dòng)變步長(zhǎng)方法。

四階五級(jí)RKF算法有參量系數(shù)表。第十五頁，共八十五頁，2022年，8月28日6.2.3.2基于MATLAB的微分方程

求解函數(shù)

格式1：直接求解

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

格式2：帶有控制參數(shù)

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

格式3：帶有附加參數(shù)

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

[t0,tf]求解區(qū)間，x0初值問題的初始狀態(tài)變量。

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

functionxd=funname(t,x)

可以使用附加變量

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

％t是時(shí)間變量或自變量（必須給），x為狀態(tài)向量，

xd為返回狀態(tài)向量的導(dǎo)數(shù)。flag用來控制求解過程,指定初值，即使初值不用指定，也必須有該變量占位。修改變量：options唯一結(jié)構(gòu)體變量，用odeset()修改。

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

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為設(shè)定的仿真終止時(shí)間>>[t,x]=ode45('lorenzeq',[0,t_final],x0);plot(t,x),>>figure;%打開新圖形窗口>>plot3(x(:,1),x(:,2),x(:,3));>>axis([1042-2020-2025]);%根據(jù)實(shí)際數(shù)值手動(dòng)設(shè)置坐標(biāo)系第十九頁，共八十五頁，2022年，8月28日可采用comet3()函數(shù)繪制動(dòng)畫式的軌跡。>>comet3(x(:,1),x(:,2),x(:,3))第二十頁，共八十五頁，2022年，8月28日描述微分方程是常微分方程初值問題數(shù)值求解的關(guān)鍵。>>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]);得出完全一致的結(jié)果。第二十一頁，共八十五頁，2022年，8月28日6.2.3.3MATLAB下帶有附加參數(shù)的微分方程求解例：第二十二頁，共八十五頁，2022年，8月28日編寫函數(shù)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位置為[]，表示不需修改控制選項(xiàng)>>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微分方程轉(zhuǎn)換

6.2.4.1單個(gè)高階常微分方程處理方法第二十五頁，共八十五頁，2022年，8月28日第二十六頁，共八十五頁，2022年，8月28日例：函數(shù)描述為：

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);

由于變步長(zhǎng)所采用的步長(zhǎng)過小，所需時(shí)間較長(zhǎng)，導(dǎo)致輸出的y矩陣過大，超出計(jì)算機(jī)存儲(chǔ)空間容量。所以不適合采用ode45()來求解，可用剛性方程求解算法ode15s()。第二十八頁，共八十五頁，2022年，8月28日6.2.4.2高階常微分方程組的變換方法第二十九頁，共八十五頁，2022年，8月28日例：第三十頁，共八十五頁，2022年，8月28日第三十一頁，共八十五頁，2022年，8月28日描述函數(shù)：

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得出的軌道不正確，默認(rèn)精度RelTol設(shè)置得太大，從而導(dǎo)致的誤差傳遞，可減小該值。第三十三頁，共八十五頁，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))%繪制出軌跡曲線顯而易見，這樣求解是錯(cuò)誤的，應(yīng)該采用更小的步長(zhǎng)。第三十七頁，共八十五頁，2022年，8月28日>>tic,[t2,y2]=rk_4('apolloeq',[0,20,0.001],x0);tocelapsed_time=124.4990％計(jì)算時(shí)間過長(zhǎng)>>plot(y2(:,1),y2(:,3))%繪制出軌跡曲線嚴(yán)格說來某些點(diǎn)仍不滿足10－6的誤差限，所以求解常微分方程組時(shí)建議采用變步長(zhǎng)算法，而不是定步長(zhǎng)算法。第三十八頁，共八十五頁，2022年，8月28日例：第三十九頁，共八十五頁，2022年，8月28日用MATLAB符號(hào)工具箱求解，令％>>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)

對(duì)于更復(fù)雜的問題來說，手工變換的難度將很大，所以如有可能，可采用計(jì)算機(jī)去求解有關(guān)方程，獲得解析解。如不能得到解析解，也需要在描寫一階常微分方程組時(shí)列寫出式子，得出問題的數(shù)值解。第四十頁，共八十五頁，2022年，8月28日6.3特殊微分方程的數(shù)值解

6.3.1剛性微分方程的求解剛性微分方程一類特殊的常微分方程，其中一些解變化緩慢，另一些變化快，且相差懸殊，這類方程常常稱為剛性方程。MATLAB采用求解函數(shù)ode15s()，該函數(shù)的調(diào)用格式和ode45()完全一致。[t,x]=ode15s(Fun,[t0,tf],x0,options,p1,p2,…)

第四十一頁，共八十五頁，2022年，8月28日例：％計(jì)算>>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)變化在某些點(diǎn)上較快。第四十三頁，共八十五頁，2022年，8月28日例：定義函數(shù)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日步長(zhǎng)分析：>>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隱式微分方程求解隱式微分方程為不能轉(zhuǎn)化為顯式常微分方程組的方程例：第四十八頁，共八十五頁，2022年，8月28日編寫函數(shù)：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微分代數(shù)方程求解例：第五十頁，共八十五頁，2022年，8月28日編寫函數(shù)

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微分代數(shù)方程中的質(zhì)量矩陣（控制參數(shù)）>>x0=[0.8;0.1;0.1];

>>[t,x]=ode15s(@c7eqdae,[0,20],x0,options);plot(t,x)第五十一頁，共八十五頁，2022年，8月28日編寫函數(shù)：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:結(jié)構(gòu)體數(shù)據(jù)，sol.x:時(shí)間向量t,sol.y:狀態(tài)向量。

第五十四頁，共八十五頁，2022年，8月28日例：第五十五頁，共八十五頁，2022年，8月28日編寫函數(shù)：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)];歷史數(shù)據(jù)函數(shù)：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邊值問題的計(jì)算機(jī)求解第五十八頁，共八十五頁，2022年，8月28日6.4.1邊值問題的打靶算法數(shù)學(xué)方法描述：以二階方程為例

第五十九頁，共八十五頁，2022年，8月28日編寫函數(shù)：線性的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日例：編寫函數(shù)：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日原方程的解析解為解的檢驗(yàn)>>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)%整個(gè)解函數(shù)檢驗(yàn)ans=4.4790e-008>>norm(y(end,1)-2)%終點(diǎn)條件檢驗(yàn)ans=2.2620e-008第六十二頁，共八十五頁，2022年，8月28日非線性方程邊值問題的打靶算法：

用Newton迭代法處理第六十三頁，共八十五頁，2022年，8月28日編寫函數(shù)：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日例：編寫兩個(gè)函數(shù)：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]);精確解：檢驗(yàn)：>>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日編寫函數(shù)：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日例：編寫函數(shù)：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偏微分方程組求解函數(shù)描述：第七十一頁，共八十五頁，2022年，8月28日邊界條件的函數(shù)描述：初值條件的函數(shù)描述：

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

第七十四頁，共八十五頁，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日描述邊界條件的函數(shù)function[pa,qa,pb,qb]=c7mpbc(xa,ua,xb,ub,t)pa=[0;ua(2)];qa=[1;0];