matlab在科學(xué)計(jì)算中的應(yīng)用6+-+微分方程問(wèn)題的解法.ppt

第六章微分方程問(wèn)題的解法,微分方程的解析解方法常微分方程問(wèn)題的數(shù)值解法微分方程問(wèn)題算法概述四階定步長(zhǎng)Runge-Kutta算法及MATLAB實(shí)現(xiàn)一階微分方程組的數(shù)值解微分方程轉(zhuǎn)換特殊微分方程的數(shù)值解邊值問(wèn)題的計(jì)算機(jī)求解偏微分方程的解,6.1微分方程的解析解方法,格式：y=dsolve(f1,f2,fm)格式：指明自變量y=dsolve(f1,f2,fm,x)fi即可以描述微分方程，又可描述初始條件或邊界條件。如：描述微分方程時(shí)描述條件時(shí),例：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),分別處理系數(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=.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);如果用推導(dǎo)的方法求Ci的值，每個(gè)系數(shù)的解析解至少要寫(xiě)出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),例：求解x,y=dsolve(D2x+2*Dx=x+2*y-exp(-t),Dy=4*x+3*y+4*exp(-t),例：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.InD:MATLAB6p5toolboxsymbolicdsolve.matline292x=t-Int(1/(a-a3+1),a=.x)+C1=0故只有部分非線性微分方程有解析解。,6.2微分方程問(wèn)題的數(shù)值解法6.2.1微分方程問(wèn)題算法概述,微分方程求解的誤差與步長(zhǎng)問(wèn)題：,functionoutx,outy=MyEuler(fun,x0,xt,y0,PointNum)%fun表示f(x,y);x0,xt:自變量的初值和終值;y0:函數(shù)在x0處的值,其可以為向量形式;PointNum表示自變量在x0,xt上取的點(diǎn)數(shù)ifnargin5|PointNumex0801h1=0.4189h2=0.2094,functionXout,Yout=MyEulerPro(fun,x0,xt,y0,PointNumber)%MyEulerPro用改進(jìn)的歐拉法解微分方程ifnargin5|PointNumberex0802,6.2.2四階定步長(zhǎng)Runge-Kutta算法及MATLAB實(shí)現(xiàn),functiontout,yout=rk_4(odefile,tspan,y0)y0初值列向量t0=tspan(1);th=tspan(2);iflength(tspan)t_final=100;x0=0;0;1e-10;%t_final為設(shè)定的仿真終止時(shí)間t,x=ode45(lorenzeq,0,t_final,x0);plot(t,x),figure;%打開(kāi)新圖形窗口plot3(x(:,1),x(:,2),x(:,3);axis(1042-2020-2025);%根據(jù)實(shí)際數(shù)值手動(dòng)設(shè)置坐標(biāo)系,可采用comet3()函數(shù)繪制動(dòng)畫(huà)式的軌跡。comet3(x(:,1),x(:,2),x(:,3),描述微分方程是常微分方程初值問(wèn)題數(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é)果。,6.2.3.3MATLAB下帶有附加參數(shù)的微分方程求解,例：,編寫(xiě)函數(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);,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微分方程轉(zhuǎn)換6.2.4.1單個(gè)高階常微分方程處理方法,例：函數(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),:),x0=2;0;t_final=3000;mu=1000;t,y=ode45(vdp_eq,0,t_final,x0,mu);由于變步長(zhǎng)所采用的步長(zhǎng)過(guò)小，所需時(shí)間較長(zhǎng)，導(dǎo)致輸出的y矩陣過(guò)大，超出計(jì)算機(jī)存儲(chǔ)空間容量。所以不適合采用ode45()來(lái)求解，可用剛性方程求解算法ode15s()。,6.2.4.2高階常微分方程組的變換方法,例：,描述函數(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)/r13-mu*(x(1)-mu1)/r23;x(4);-2*x(2)+x(3)-mu1*x(3)/r13-mu*x(3)/r23;,不用end結(jié)尾也可以定義函數(shù),求解：x0=1.2;0;0;-1.04935751;tic,t,y=ode45(apolloeq,0,20,x0);tocelapsed_time=0.8310length(t),plot(y(:,1),y(:,3)ans=689得出的軌道不正確，默認(rèn)精度RelTol設(shè)置得太大，從而導(dǎo)致的誤差傳遞，可減小該值。,改變精度：options=odeset;options.RelTol=1e-6;tic,t1,y1=ode45(apolloeq,0,20,x0,options);tocelapsed_time=0.8110length(t1),plot(y1(:,1),y1(:,3),ans=1873,min(diff(t1)ans=1.8927e-004plot(t1(1:end-1),diff(t1),例：,x0=1.2;0;0;-1.04935751;tic,t1,y1=rk_4(apolloeq,0,20,0.01,x0);tocelapsed_time=4.2570plot(y1(:,1),y1(:,3)%繪制出軌跡曲線顯而易見(jiàn)，這樣求解是錯(cuò)誤的，應(yīng)該采用更小的步長(zhǎng)。,tic,t2,y2=rk_4(apolloeq,0,20,0.001,x0);tocelapsed_time=124.4990計(jì)算時(shí)間過(guò)長(zhǎng)plot(y2(:,1),y2(:,3)%繪制出軌跡曲線嚴(yán)格說(shuō)來(lái)某些點(diǎn)仍不滿足106的誤差限，所以求解常微分方程組時(shí)建議采用變步長(zhǎng)算法，而不是定步長(zhǎng)算法。,例：,用MATLAB符號(hào)工具箱求解，令symsx1x2x3x4dx,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)對(duì)于更復(fù)雜的問(wèn)題來(lái)說(shuō)，手工變換的難度將很大，所以如有可能，可采用計(jì)算機(jī)去求解有關(guān)方程，獲得解析解。如不能得到解析解，也需要在描寫(xiě)一階常微分方程組時(shí)列寫(xiě)出式子，得出問(wèn)題的數(shù)值解。,6.3特殊微分方程的數(shù)值解6.3.1剛性微分方程的求解,剛性微分方程一類特殊的常微分方程，其中一些解變化緩慢，另一些變化快，且相差懸殊，這類方程常常稱為剛性方程。MATLAB采用求解函數(shù)ode15s()，該函數(shù)的調(diào)用格式和ode45()完全一致。t,x=ode15s(Fun,t0,tf,x0,options,p1,p2,),例：計(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,作圖plot(t,y(:,1);figure;plot(t,y(:,2)y(:,1)曲線變化較平滑,y(:,2)變化在某些點(diǎn)上較快。,例：定義函數(shù)functiondy=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);tocelapsed_time=229.4700length(t2),plot(t2,y2)ans=356941,步長(zhǎng)分析：formatlong,min(diff(t2),max(diff(t2)ans=0.000222206938840.00214971787184plot(t2(1:end-1),diff(t2),方法二，用ode15s()代替ode45()opt=odeset;opt.RelTol=1e-6;tic,t1,y1=ode15s(c7exstf2,0,100,0;1,opt);tocelapsed_time=0.49100000000000length(t1),plot(t1,y1)ans=169,6.3.2隱式微分方程求解,隱式微分方程為不能轉(zhuǎn)化為顯式常微分方程組的方程例：,編寫(xiě)函數(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),6.3.3微分代數(shù)方程求解,例：,編寫(xiě)函數(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),編寫(xiě)函數(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);,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:結(jié)構(gòu)體數(shù)據(jù)，sol.x:時(shí)間向量t,sol.y:狀態(tài)向量。,例：,編寫(xiě)函數(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);,求解：lags=10.5;tx=dde23(c7exdde,lags,zeros(3,1),0,10);plot(tx.x,tx.y(2,:)與ode45()等返回的x矩陣不一樣，它是按行排列的。,6.4邊值問(wèn)題的計(jì)算機(jī)求解,6.4.1邊值問(wèn)題的打靶算法,數(shù)學(xué)方法描述：以二階方程為例,編寫(xiě)函數(shù)：線性的functiont,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);,例：編寫(xiě)函數(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),原方程的解析解為解的檢驗(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-008norm(y(end,1)-2)%終點(diǎn)條件檢驗(yàn)ans=2.2620e-008,非線性方程邊值問(wèn)題的打靶算法：用Newton迭代法處理m,編寫(xiě)函數(shù)：functiont,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);endt,y=ode45(funcn,tspan,ga;m,varargin);,例：編寫(xiě)兩個(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);,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-005norm(y(end,1)-1)ans=5.2815e-006,6.4.2線性微分方程的有限差分算法,把等式左邊用差商表示。,編寫(xiě)函數(shù)：functionx,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+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);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;,例：編寫(xiě)函數(shù)：functiony=c7fun5(x,key)switchkeycase1,y=1+x;case2,y=1-x;otherwise,y=1+x.2;endt,y=fdiff(c7fun5,0,1,1,4,50);plot(t,y),6.5偏微分方程求解入門(mén)6.5.1偏微分方程組求解,函數(shù)描述：,邊界條件的函數(shù)描述：初值條件的函數(shù)描述：u0=pdeic(x),例：,函數(shù)描述：,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);,描述邊界條件的函數(shù)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),6.5.2二階偏微分方程的數(shù)學(xué)描述,橢圓型偏微分方程：,拋物線型偏微分方程：雙曲型偏微分方程：特征值型偏微分方程：,6.5.3偏微分方程的求解界面應(yīng)用簡(jiǎn)介6.5.3.1偏微分方程求解程序概述,啟動(dòng)偏微分方程求解界面在MATLAB下鍵入pdetool該界面分為四個(gè)部分菜單系統(tǒng)工具欄集合編輯求解區(qū)域,6.5.3.2偏微分方程求解區(qū)域繪制,1）用工具欄中的橢圓、矩形等繪制一些區(qū)域。2）在集合編輯欄中修改其內(nèi)容。如（R1E1E2）E33）單擊工具欄中按紐可得求解邊界。4）選擇Boundary-RemoveAllSubdomainBorders菜單項(xiàng)，消除相鄰區(qū)域