第9章微分方程

1第9章

MATLAB微分方程求解

29.1常微分方程初值問題常微分方程的初值問題是科學(xué)計(jì)算中的常見問題，MATLAB提供了求解常微分方程（ODE）初值問題的一系列ode函數(shù)。ode函數(shù)主要采用Runge-Kutta法求解常微分方程。常用ode函數(shù)函數(shù)名含義函數(shù)名含義ode23普通2-3階法解ODEode23s變階法解剛性odeode45普通4-5階法解ODEode23t解適度剛性odeode113普通變階法解ODEode23tb低階法解剛性ode39.1.1ODE的機(jī)理將高階常微分方程的初值問題，轉(zhuǎn)變?yōu)?階常微分方程組的初值問題。9.1.2普通2-3階法解ODE例9-1編制常微分方程的函數(shù)文件dfun.mfunctiondy=dfun(x,y)dy=y-2*x/y編制如下計(jì)算程序clear;clc;X12=linspace(0,5,50);[X,Y]=ode23('dfun',X12,2);plot(X,Y)運(yùn)行結(jié)果如圖9-1所示圖9-149.1.3普通4-5階法解ODE例9-2編制微分方程組的函數(shù)文件dfun.mfunctiondy=dfun(x,y)dy(1,1)=-2*y(1)+y(2)+2*sin(x);dy(2,1)=10*y(1)-9*y(2)+9*(cos(x)-sin(x));編制如下程序clear;clc;X12=linspace(0,9,30);ode45('dfun',X12,[2,3]);運(yùn)行結(jié)果如圖9-2所示圖9-259.2常微分方程邊值問題9.2.1邊值問題簡介一般來說，微分方程邊值問題可能有解、也可能無解，可能有唯一解、也可能有無數(shù)解。在假定有唯一解的前提下，邊值問題有三種基本解法：迭加法：假如微分方程和邊界條件是線性的，邊值問題可以轉(zhuǎn)化為初值問

題。試射法：將問題轉(zhuǎn)化為對漏缺初值的搜索，一旦漏缺初值確定，就可以按

初值問題求解。松弛法：先猜測滿足邊界條件的區(qū)間網(wǎng)點(diǎn)上的解值，然后利用微分方程進(jìn)

行迭代改善。69.2.2bvp4c求解思路為求解常微分方程的邊值問題，MATLAB提供了bvp4c函數(shù)文件。該文件是依據(jù)有限元法中的配置法編寫的。其解題的基本思路如下：

把待解問題轉(zhuǎn)化為標(biāo)準(zhǔn)邊值問題為期望解指定初始猜測利用原微分方程，構(gòu)造殘差函數(shù)用原微分方程和邊界條件，借助迭代使殘差不斷減小為了bvp4c使用的方便，MATLAB提供了能產(chǎn)生較好猜測解的bvpinit指令，減輕了用戶進(jìn)行初始猜測的負(fù)擔(dān)。此外解算指令bvp4c也能處理包含未知參數(shù)的微分方程的邊值問題。79.2.3bvp4c配套指令例9-3編制微分方程組的函數(shù)文件dfun.mfunctiondy=dfun(x,y,c)dy(1,1)=y(2);dy(2,1)=-c*abs(y(1));編制殘差方程組的函數(shù)文件bfun.mfunctionres=bfun(ya,yb,c)res(1,1)=ya(1);res(2,1)=yb(1)+2編制如下計(jì)算程序sol=bvp4c('dfun','bfun',sinit,[],c);plot(sol.x,sol.y(1,:),'r--')holdonc=1.5;sol=bvp4c('dfun','bfun',sinit,[],c);plot(sol.x,sol.y(1,:),'b-')legend('c=1','c=1.5')運(yùn)行結(jié)果如圖9-3所示clear;clc;x=linspace(0,4,5);v=[1;0];sinit=bvpinit(x,v);c=1;圖9-389.3常微分方程的解析解9.3.1微分方程的通解例9-4編制如下程序clear;clc;fun='Dx=-a*t'x=dsolve(fun,'t')運(yùn)行結(jié)果為fun=Dx=-a*t

x=

-1/2*a*t^2+C1例9-5編制如下程序clear;clc;fun1='Df=f+g';fun2='Dg=g-f';[f,g]=dsolve(fun1,fun2,'t')運(yùn)行結(jié)果為f=

-exp(t)*(C1*cos(t)-C2*sin(t))

g=

exp(t)*(C1*sin(t)+C2*cos(t))99.3.2微分方程的特解例9-6編制如下程序clear;clc;fun='D2y=4'x=dsolve(fun,'y(0)=1','y(1)=3','x')運(yùn)行結(jié)果為fun=D2y=4

x=

2*x^2+1例9-7編制如下程序clear;clc;fun1='Df=f+g';fun2='Dg=g-f';v1='f(0)=2';v2='g(0)=2';[f,g]=dsolve(fun1,fun2,v1,v2,'t')運(yùn)行結(jié)果為f=

exp(t)*(2*cos(t)+2*sin(t))

g=

exp(t)*(-2*sin(t)+2*cos(t))10例9-8編制如下程序clear;clc;fun='D3y=-y';v1='y(0)=3';v2='Dy(0)=0';v3='D2y(0)=0';y=dsolve(fun,v1,v2,v3,'t')運(yùn)行結(jié)果為y=

exp(-t)+2*exp(1/2*t)*cos(1/2*3^(1/2)*t)119.4偏微分方程求解9.4.1偏微分方程簡介很多科學(xué)計(jì)算問題都可以歸結(jié)為求解含定解條件(初始條件和邊界條件)的偏微分方程問題。MATLAB提供了專門求解偏微分方程問題工具箱(PDEToolbox)。MATLAB可以求解的線性偏微分方程類型包括：

橢圓型方程拋物線方程雙曲型方程特征值方程129.4.2偏微分方程的求解過程利用MATLAB求解偏微分方程的過程如下：1.完成幾何造型2.輸入定解條件3.選擇方程類型4.劃分網(wǎng)格5.數(shù)值求解6.結(jié)果可視化13例9-91415169.4.3PDE圖形用戶界面簡介例9-10179.5綜合實(shí)例9.5.1有阻尼自由振動例9-11編制標(biāo)準(zhǔn)一階常微分方程組的函數(shù)文件dfun.mfunctiondz=dfun(t,z)a=-2;b=-100;dz(1,1)=z(2);dz(2,1)=a*z(2)+b*z(1);編制如下程序clear;clc;t12=linspace(0,2,100);[t,z]=ode45('dfun',t12,[1,0]);plot(t,z(:,1))gridon運(yùn)行后振動規(guī)律如圖9-19所示圖9-19189.5.2點(diǎn)的運(yùn)動軌跡例9-12編制標(biāo)準(zhǔn)微分方程組的函數(shù)文件dfun.mfunctiondu=dfun(t,u)du(1,1)=u(3);du(2,1)=u(4);du(3,1)=0.2*u(1);du(4,1)=0.2*u(2);編制如下程序clear;clc;t12=linspace(0,2,100);[t,u]=ode45('dfun',t12,[2,0,0,1])plot(u(:,1),u(:,2))gridon運(yùn)行后，得到軌跡如圖9-20所示圖9-20199.5.3梁的撓曲線編制標(biāo)準(zhǔn)微分方程組的函數(shù)文件dfun.mfunctiondu=dfun(x,u)du(1,1)=u(2);du(2,1)=-0.02*x+0.01*x^3;根據(jù)定解條件,編制殘差方程組的函數(shù)文件bfun.mfunctionres=bfun(ya,yb)res(1,1)=ya(1);r