《MATLAB解微分方程》課件

MATLABODE初值問題的數(shù)值解PDE問題的數(shù)值解1精選課件ppt問題提出倒葫蘆形狀容器壁上的刻度問題.對于如圖所示圓柱形狀容器壁上的容積刻度,可以利用圓柱體體積公式其中直徑D為常數(shù).對于幾何形狀不是規(guī)則的容器,比如倒葫蘆形狀容器壁上如何標(biāo)出刻度呢?下表是經(jīng)過測量得到部分容器高度與直徑的關(guān)系.2精選課件pptx1o根據(jù)上表的數(shù)據(jù),可以擬合出倒葫蘆形狀容器的圖,建立如圖所示的坐標(biāo)軸,問題為如何根據(jù)任意高度x標(biāo)出容器體積V的刻度,由微元思想分析可知H00.20.40.60.81.0D0.040.110.260.561.041.17其中x表示高度,直徑D是高度x的函數(shù),記為D(x).x1o3精選課件ppt只要求解上述方程,就可求出體積V與高度x之間的函數(shù)關(guān)系,從而可標(biāo)出容器壁上容積的刻度,但問題是函數(shù)D(x)無解析表達(dá)式,無法求出其解析解.

因此,得到如下微分方程初值問題4精選課件ppt包含自變量、未知函數(shù)及未知函數(shù)的導(dǎo)數(shù)或微分的方程稱為微分方程。在微分方程中,自變量的個數(shù)只有一個,稱為常微分方程。自變量的個數(shù)為兩個或兩個以上的微分方程叫偏微分方程。微分方程中出現(xiàn)的未知函數(shù)最高階導(dǎo)數(shù)的階數(shù)稱為微分方程的階數(shù)。微分方程分類5精選課件ppt常微分方程：(1)，(2)式稱為初值問題.在實際應(yīng)用中還經(jīng)常需要求解常微分方程組：(3)式稱為邊值問題。6精選課件ppt但能求解析解的常微分方程是有限的，大多數(shù)的常微分方程是給不出解析解的.這個一階微分方程就不能用初等函數(shù)及其積分來表達(dá)它的解。

例例

的解的值仍需插值方法來計算.7精選課件ppt8精選課件ppt事實上，從實際問題當(dāng)中抽象出來的微分方程，通常主要依靠數(shù)值解法來解決?？梢宰C明:如果函數(shù)在帶形區(qū)域R=a≤x≤b,-∞＜y＜∞}內(nèi)連續(xù)，且關(guān)于y滿足李普希茲(Lipschitz)條件，即存在常數(shù)L(它與x,y無關(guān))使

對R內(nèi)任意兩個都成立,則方程的解在

a,b

上存在且唯一。

在區(qū)間a≤x≤b上的數(shù)值解法。

主要討論一階常微分方程初值問題9精選課件ppt微分方程數(shù)值方法的基本思想對常微分方程初值問題的數(shù)值解法，就是要算出精確解y(x)在區(qū)間

a,b

上的一系列離散節(jié)點

處的函數(shù)值相鄰兩個節(jié)點的間距稱為步長，步長可以相等，也可以不等。假定h為定數(shù)，稱為定步長，這時節(jié)點可表示為10精選課件ppt數(shù)值解法需要把連續(xù)性的問題加以離散化，從而求出離散節(jié)點的數(shù)值解。離散化對常微分方程數(shù)值解法的基本出發(fā)點就是離散化。其數(shù)值解法的基本特點：采用“步進(jìn)式”，即求解過程順著節(jié)點排列的次序一步一步地向前推進(jìn)，11精選課件ppt描述這類算法，要求給出用已知信息計算的遞推公式。建立這類遞推公式的基本方法是在這些節(jié)點上用數(shù)值積分、數(shù)值微分、泰勒展開等離散化方法，對初值問題中的導(dǎo)數(shù)進(jìn)行不同的離散化處理。12精選課件ppt數(shù)值解和精確解用數(shù)值方法求解初值問題，不是求出它的解析解或其近似解析式，而是給出它的解在某些離散節(jié)點上的近似值。用y(x)表示問題的準(zhǔn)確解y(x0),y(x1),y(xN)表示解y(x)在節(jié)點x0,x1,…,xN處的準(zhǔn)確值y0,y1,…,yN表示數(shù)值解，即問題的解y(x)在相應(yīng)節(jié)點處的近似值。13精選課件ppt單步法和多步法單步法：在計算yi+1時只利用yi多步法：在計算yi+1時不僅利用yi,還要利用yi?1,yi?2,…,k步法：在計算yi+1時要用到y(tǒng)i,yi?1,…,yi?k+1顯式格式可寫成:yk+1=yk+hΦf(xk，yk;h）隱式格式:yk+1=yk+hΦf（xk，yk,yk+1;h）它每步求解yk+1需要解一個隱式方程。14精選課件ppt歐拉（Euler）方法在x=x0處，用差商代替導(dǎo)數(shù)：由得15精選課件ppt同理，在x=xn

處，用差商代替導(dǎo)數(shù)：由得若記則上式可記為此即為求解初值問題的Euler方法，又稱顯式Euler方法。16精選課件ppt例:用Euler方法求解常微分方程初值問題并將數(shù)值解和該問題的解析解比較。解：Euler方法的具體格式：17精選課件ppth=0.2;y(1)=0.2;x=0.2:h:3;forn=1:14xn=x(n);yn=y(n);y(n+1)=yn+h*(yn/xn-2*yn*yn);endx0=0.2:h:3;y0=x0./(1+x0.^2);plot(x0,y0,x,y,x,y,'o')程序?qū)崿F(xiàn)18精選課件ppt

xn

y(xn) yn

yn-y(xn) 0.0 0 0 00.2 0.1923 0.2000 0.00770.4 0.3448 0.3840 0.03920.6 0.4412 0.5170 0.07580.8 0.4878 0.5824 0.0946 1.0 0.5000 0.5924 0.09241.2 0.4918 0.5705 0.07871.4 0.4730 0.5354 0.0624h=0.2,xn=nh,(n=0,1,2…,15),f(x,y)=y/x–2y2計算中取f(0,0)=1.計算結(jié)果如下：19精選課件pptxn

y(xn) yn

yn-y(xn)1.6 0.4494 0.4972 0.0478 1.8 0.4245 0.4605 0.03592.0 0.4000 0.4268 0.02682.2 0.3767 0.3966 0.0199 2.4 0.3550 0.3698 0.01472.6 0.3351 0.3459 0.01082.8 0.3167 0.3246 0.0079 3.0 0.3000 0.3057 0.0057由表中數(shù)據(jù)可以看到，微分方程初值問題的數(shù)值解和解析解的誤差一般在小數(shù)點后第二位或第三位小數(shù)上，說明Euler方法的精度是比較差的。20精選課件ppt二階Runge-Kutta方法其中c1,c2,

2,

21待定。上式的局部截斷誤差為：21精選課件ppt即c1=1-a，

2=

21=1/(2a)方程組解不唯一，可令c2=a

0

，則

滿足上述條件的公式都為2階R-K公式。22精選課件ppt例

蛇形曲線的初值問題令f(x,y)=y/x–2y2,取f(0,0)=1,h=0.2,xn=hn,(n=1,2,…,15)2階龍格-庫塔公式計算格式：k1=yn/xn–2yn

2,k2=(yn+hk1)/(xn+h)–2(yn+hk1)2yn+1=yn+0.5h[k1+k2]23精選課件pptx0=0;y0=0;h=.2;x=.2:h:3;k1=1;k2=(y0+h*k1)/x(1)-2*(y0+h*k1)^2;y(1)=y0+.5*h*(k1+k2);forn=1:14k1=y(n)/x(n)-2*y(n)^2;k2=(y(n)+h*k1)/x(n+1)-2*(y(n)+h*k1)^2;y(n+1)=y(n)+0.5*h*(k1+k2);endy1=x./(1+x.^2);plot(x,y,'o',x,y1)24精選課件ppt常用的一個公式為四階Runge-Kutta方法25精選課件pptfunctionydot=harmonic(t,y)ydot=[y(2);-y(1)]y=inline(‘[01;-10]*y’,’t’,’y’);SystemofEquations26精選課件pptfunctionydot=twobody(t,y)r=sqrt(y(1)^2+y(2)^2);ydot=[y(3);y(4);-y(1)/r^3;-y(2)/r^3];

TwoBodyProblem27精選課件pptLinearizedDifferentialEquations28精選課件pptJ的特征值是解增長解衰減解振蕩29精選課件ppt基于龍格－庫塔法，MATLAB求常微分方程數(shù)值解的函數(shù)，一般調(diào)用格式為：[t,y]=ode23('fname',tspan,y0)[t,y]=ode45('fname',tspan,y0)其中fname是定義f(t,y)的函數(shù)文件名，該函數(shù)文件必須返回一個列向量。tspan形式為[t0,tf],表示求解區(qū)間。y0是初始狀態(tài)列向量。t和y分別給出時間向量和相應(yīng)的狀態(tài)向量。MATLAB求常微分方程數(shù)值解的函數(shù)30精選課件pptode23:Bogacki,Shampine(1989)和Shampine(1994),”23”表示用兩算法：一個2階，一個3階Bogacki,P.andLFShampine,"A3(2)pairofRunge-Kuttaformulas,"Appl.

Math.Letters,Vol.2,1989,pp1-9.

BS23algorithm31精選課件pptF=inline('[y(2);-y(1)]','t','y')ode23(F,[02*pi],[1;0])opts=odeset(‘reltol’,1.e-4,’abstol’,1.e-6,’outputfcn’)Examples32精選課件pptode23(@twobody,[02*pi],[1;0;0;1]);01234567-1.5-1-0.500.511.5Examples33精選課件ppty0=[1;0;0;3];ode23(@twobody,[02*pi],y0);01234567-202468101214161834精選課件ppty0=[1;0;0;3];[t,y]=ode23(@twobody,[02*pi],y0);plot(y(:,1),y(:,2));axisequal35精選課件ppty0=[1;0;0;3];[t,y]=ode23(@twobody,[02*pi],y0);plot(y(:,1))plot(y(:,2))36精選課件pptAproblemisstiffifthesolutionbeingsoughtisvaryingslowly,buttherearenearbysolutionsthatveryrapidly,sothenumericalmethodmusttakesmallstepstoobtainsatisfactoryresults.Amodelofflamepropagation（火焰燃燒）:y是球的半徑，y^2和y^3與球的表面積和體積有關(guān)想一下，點燃一根火柴，光球迅速增長，到達(dá)一個關(guān)鍵的大小，然后維持它的大?。ㄓ捎谶M(jìn)入球內(nèi)氧氣和消耗的氧氣平衡）StiffProblem(剛性問題)37精選課件ppt010203040506070809010000.20.40.60.811.21.4eta=0.02;symy;F=inline(‘y^2-y^3’,’t’,’y’);ode23(F,[02/eta],eta);38精選課件ppteta=0.00002;symy;F=inline(‘y^2-y^3’,’t’,’y’);ode23(F,[02/eta],eta);39精選課件ppt012345678910x10400.20.40.60.811.21.4ode23seta=0.00002;ode23s(inline('y^2-y^3','t','y'),[02/eta],eta);40精選課件ppt例蛇形曲線的常微分方程初值問題MATLAB數(shù)值求解命令F=inline('1./(1+x.^2)-2*y.^2')；ode23(F,[0,6],0)輸出結(jié)果為圖形41精選課件ppt[T,y]=ode23(f,[0,6],0)將得到自變量和函數(shù)的離散數(shù)據(jù)T=00.00010.00050.00250.01250.04960.10850.18630.28370.40910.59910.85131.05671.26801.51101.80502.17882.68423.28423.88424.48425.08425.68426.0000y=00.00010.00050.00250.01250.04950.10730.18000.26260.35050.44110.49440.50000.48680.46070.42420.37930.32700.27830.24110.21220.18910.17050.162042精選課件ppt例洛倫茲模型由如下常微分方程組描述取

=8/3，

=10，

=28。初值：x(0)=0，y(0)=0，z(0)=0.01。利用MATLAB求解常微分方程數(shù)值解命令計算出t∈[0，80]內(nèi)，三個未知函數(shù)的數(shù)據(jù)值，并繪出相空間在Y-X平面的投影曲線43精選課件ppt氣象學(xué)家Lorenz提出一篇論文，名叫「一只蝴蝶拍一下翅膀會不會在Taxas州引起龍卷風(fēng)？」論述某系統(tǒng)如果初期條件差一點點，結(jié)果會很不穩(wěn)定，他把這種現(xiàn)象戲稱做「蝴蝶效應(yīng)」。Lorenz為何要寫這篇論文呢？這故事發(fā)生在1961年的某個冬天，他如往常一般在辦公室操作氣象電腦。平時，他只需要將溫度、濕度、壓力等氣象數(shù)據(jù)輸入，電腦就會依據(jù)三個內(nèi)建的微分方程式，計算出下一刻可能的氣象數(shù)據(jù)，因此模擬出氣象變化圖。這一天，Lorenz想更進(jìn)一步了解某段紀(jì)錄的後續(xù)變化，他把某時刻的氣象數(shù)據(jù)重新輸入電腦，讓電腦計算出更多的後續(xù)結(jié)果。當(dāng)時，電腦處理數(shù)據(jù)資料的數(shù)度不快，在結(jié)果出來之前，足夠他喝杯咖啡并和友人閑聊一陣。在一小時後，結(jié)果出來了，不過令他目瞪口呆。結(jié)果和原資訊兩相比較，初期數(shù)據(jù)還差不多，越到後期，數(shù)據(jù)差異就越大了，就像是不同的兩筆資訊。而問題并不出在電腦，問題是他輸入的數(shù)據(jù)差了0.000127，而這些微的差異卻造成天壤之別。所以長期的準(zhǔn)確預(yù)測天氣是不可能的。44精選課件ppt天氣預(yù)報的準(zhǔn)確性：.tw/senior/seniorteach/stfastnews/stnfastnews/stnfastnews4/stnfastnews4_2.htmLorenz現(xiàn)象的數(shù)學(xué)：/news_detail.php?news_id=225分形藝術(shù)電子版：/personal/huajie/fractalart/html/notice.htm混沌理論：/~yzhao/doc/chaos.html45精選課件ppt記向量[y1，y2，y3]=[x，y，z]，創(chuàng)建MATLAB函數(shù)文件如下functionz=flo(t,y)z(1,:)=-8*y(1)/3+y(2).*y(3);z(2,:)=-10*(y(2)-y(3));z(3,:)=-y(1).*y(2)+28*y(2)-y(3);用MATLAB命令求解并繪出Y-X平面的投影圖

y0=[0;0;0.01];[x,y]=ode23(@flo,[0,80],y0)plot(y(:,2),y(:,1))46精選課件pptplot(y(:,3),y(:,1))47精選課件pptplot3(y(:,1),y(:,2),y(:,3))48精選課件ppt

y0=[30;0;-40];plot(y(:,i))49精選課件ppt50精選課件ppt51精選課件ppt非剛性系統(tǒng)：ode45(Runge-Kutta45)ode23(Runge-Kuatta23)ode113(Adams-Bashforth-MoultonPECE)多步方法剛性系統(tǒng)：ode15s(Gear方法),多步方法ode23s(二階modifiedRosenbroackformula)，單步ode23t(trapezoidalrule),solveDAEsode23tb(TR-BDF2)lowordermethodMatlab'sODESolvers52精選課件pptLaplacian算子：Poisson方程（elliptic）：Laplacian算子的特征值問題:Heatequation(parabolic):Waveequation(hyperbolic):PDEModel53精選課件ppt五點離散Poisson方程離散：特征值問題:FineteDifferenceMethods54精選課件ppt熱方程：波動方程：55精選課件ppt-41100000000000001-40100000000000010-41100000000000011-40100000000000010-41100000000000011-40100000000000010-41000100000000011-41000100000000001-41000100000000001-41000100000000001-40000100000000000-41000000000000001-41000000000000001-4100000