代碼資源精通matlab科學(xué)計(jì)算12-3rd

12。常微分方程問(wèn)題來(lái)近似求解中提供了求解函數(shù)和求解器可以實(shí)現(xiàn)常微分方程。通過(guò)本章學(xué)習(xí)，讀者能夠熟練使用的求解函數(shù)和求解器，以及通過(guò)編程進(jìn)中的求解中的求解函的影響。當(dāng)常微分方程式能夠解析求解時(shí)，可用的符號(hào)工具箱中的功能找到精確解在常微分方程難以獲得解析解的情況下使用的常微分方程求解器solver，符號(hào)解函數(shù)y'g(x,xyxy=f(x,y'=f'=g(x,y)y(x0y0下，可以得到唯一解。常微分方程符號(hào)解的語(yǔ)法是：表二階微分項(xiàng)y''，condition則為初始條件。dsolvedsolve('equation',dsolve('equation',vdsolve('equation','condition',v12-1】1。計(jì)算微分方程dy3xyxex>>dsolve('Dy+3*x*y=x*exp(-ans=(1/3*exp(-x*(x-3*t))+C1)*exp(->>dsolve('Dy+3*x*y=x*exp(-ans=exp(-x^2)+exp(-y

3e

C1，其中C1【例12-2】 常微分方程符號(hào)解求解實(shí)例2。計(jì)算微分方程xy'2yex0在初始條件y|x12e下的特解。>>dsolve('x*Dy+2*y-ans=(exp(x)*x-y

xexex2e12-3】3y''2y'ex0>>ans=-1/3*exp(x)-1/2*exp(-可知通解為y1ex1e2xC1C2，其中CC2為常數(shù) 求解器tspan=[t0tf]y0yf(tytodefunyf(ty中的f(ty，tspant他指定點(diǎn)t0t1t

上的解，則令tspant0t1,t2

，y0solverode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb之一，其中ode45、ode23、ode113ODE12-1所示；ode15s、ode23s、ode23t、ode23tb屬于剛性O(shè)DE類(lèi)型，這些命令的特點(diǎn)如表12-2所示。12-1非剛性O(shè)DE求解命求解器 一步算法；4，5階龍格-庫(kù)塔方程；累計(jì)截?cái)嗾`差達(dá)大部分場(chǎng)合的首選算一步算法；2，3階龍格-庫(kù)塔方程；累計(jì)截?cái)嗾`差達(dá)使用于精度較低的情多步法；Adams算法；高低精度均可到10-3～10-計(jì)算時(shí)間比ode4512-2剛性O(shè)DE求解求解器 梯形算適度剛性情多步法；Gear’s反向數(shù)值微分；精度中ode45失效時(shí)，可嘗試使一步法；2Rosebrock算法；低精當(dāng)精度較低時(shí)，計(jì)算時(shí)間比ode15s梯形算法；低精當(dāng)精度較低時(shí)，計(jì)算時(shí)間比ode15sode45ode23ode45ode23完全一樣。兩個(gè)函數(shù)的差別在于必須與所用的內(nèi)部算法相關(guān)。兩個(gè)函數(shù)都運(yùn)用了基本的龍格-庫(kù)塔ode232/3階龍格-庫(kù)塔-芬爾格(Runge-Kutta-Fehlerg)ode45運(yùn)用組合的4/5階龍格-庫(kù)塔-芬爾格算法。通常，ode45ode23t0tf之間可取較少的時(shí)間步。然而，在同一時(shí)間內(nèi)，ode233ode45每時(shí)間步至少調(diào)用6次。正如使用高階多項(xiàng)式內(nèi)插常常得不到最好的結(jié)果一樣，ode45也不總是比ode23好。如果ode45產(chǎn)生的結(jié)果，對(duì)作圖間隔太大，則必須在更細(xì)的時(shí)間區(qū)間內(nèi)對(duì)數(shù)據(jù)進(jìn)行內(nèi)插，比如用函數(shù)interp1。這個(gè)附加時(shí)間點(diǎn)會(huì)使ode23更有效。作為一條普遍規(guī)則，在所計(jì)算的solverODEF(y,y',...y(n1),t)y(0)y0,y'(0)y1,...y(n1)(0)寫(xiě)為向量的形式為yyy(1y(2y(m1，nmy1

f1(t,y)y2 f2(t,y fn(t,

y1(0) y0 y2 y0 yn odefilesolver12-4】solverode45y'y

y(0) 2

y'y ，其初始條件為y(0)1 1 y'0.51y

y3(0) 1 解：編寫(xiě)函數(shù)文件functiondy=rigid(t,y)dy=zeros(3,1); 行向量dy(1)=y(2)*y(3);dy(2)=-y(1)*dy(3)=-0.51*y(1)*建立文件ex1204.moptions=odeset('RelTol',1e-4,'AbsTol',[1e-41e-41e-5]);[T,Y]=ode45(@rigid,[012],[011],options);12-1100-0-1

12-1非剛性體運(yùn)動(dòng)的微分方程12.2

y'f(x,y(x0)y'(x)y(xh)hyn1ynhf(xn,yn 調(diào)用格式：[outx,outyMyEuler(fun,x0,xty0,PointNum)。其中，fun為函數(shù)名；x0為自變量區(qū)間初值；xt為自變量區(qū)間終值；y0x0PointNum為自變量在[x0,xt]上所取的點(diǎn)數(shù)；outx為輸出自變量區(qū)間上所取點(diǎn)的x值；outy為對(duì)應(yīng)點(diǎn)上的函數(shù)值。歐拉法 function%用前向差分的歐拉方法解微分方程%函數(shù)%y0表示函數(shù)在x0%outx:所取的點(diǎn)的xifnargin<5|

ifnargin<4 %y0默認(rèn)值為0 %計(jì)算步長(zhǎng) y(1，：)=y0（：）’; fork=1:PointNum %計(jì)算f(x,y)在每個(gè)迭代點(diǎn)的y(k1,y(kh*f;%對(duì)于所取的點(diǎn)x迭代計(jì)算y 12-5】y'sinxy(x)1,x h0.2h0.4解：首先建立函數(shù)文件functionf=myfun01(x,y)ex1201.mfunctionex1205() h1=2*pi/15 h2=pi/15%計(jì)算步長(zhǎng) for 12-212-2歐拉法所得方程解與符號(hào)解（精確解）對(duì)12-2中可以看出，用歐拉法得到的解和用符號(hào)法得到的解之間存在一定的誤差，差商和中心差商。有的讀者可以參考相關(guān)資料，編寫(xiě)相應(yīng)的程序，其基本方法，又稱(chēng)為Henu算法。對(duì)方程y'f(x,y)兩邊由xn到xn1積分，并利用梯形，可得yn1ynh[f(xn,yn)f(xn1,yn1)]y(x0)具體迭代如下

yn1ynhf(xn,ynyn1ynh[f(xn,yn)f(xn1,yn1)]在中編程實(shí)現(xiàn)的改進(jìn)的歐拉法函數(shù)為：MyEulerPro。調(diào)用格式：[Xout,Yout]=MyEulerPro(fun,x0,xty0,PointNumber)。其中，fun為函數(shù)名；x0為自變量區(qū)間初值；xt為自變量區(qū)間終值；y0x0Xoutx改進(jìn)的歐拉法的程序代碼如下function用改進(jìn)的歐拉法解微分方程%函數(shù)%y0表示函數(shù)在x0%Xout:所取的點(diǎn)的xifnargin<5|PointNumer<=0 %如果函數(shù)僅輸4個(gè)參數(shù)值，則PointNumer默認(rèn)值為100ifnargin<4%y0默認(rèn)值為0 %表示出離散的自變量xfori=1:PointNumber 【例12-6】 y'sinxy(x)1,x 即對(duì)【例解：編寫(xiě)程序function%ex1206比較改進(jìn)的歐拉法，簡(jiǎn)單歐拉方法以及微分方程符號(hào)解 for12-312-312-3改進(jìn)的歐拉法解、歐拉解、

由此可見(jiàn)，改進(jìn)的歐拉法較之簡(jiǎn)單歐拉法更為龍格-庫(kù)塔yn1ynh(k12k22k3k46k1f(xn,ynk2f(xnk3f(xn

1h,yn 1h,yn

hk2k4f(xnh,ynhk3在中編程實(shí)現(xiàn)的四階龍格-庫(kù)塔算法函數(shù)為：MyRunge_Kutta。調(diào)用格式：[xy]=MyRunge_Kutta(fun,x0,xty0,PointNumvarargin)。其中，fun為函數(shù)名；x0為自變量區(qū)間初值；xt為自變量區(qū)間終值；y0x0PointNum為自變量在[x0,xt]xxy四階龍格-庫(kù)塔法的程序代碼如下function[x,y]=%Runge-Kutta方法解微分方程形為%函數(shù)%y0表示函數(shù)在x0%varargin為可選輸入項(xiàng)可傳適當(dāng)參數(shù)給函數(shù)%x:所取的點(diǎn)的xifnargin<4|PointNum<=0PointNum=100;ifnargin<3y0=0;y(1 hxt-x0)/(PointNum- x= %得x向量fork= f1=f1= %得中f2=h*feval(fun,x(k)+h/2,y(k,:)+f2= %得中f3=h*feval(fun,x(k)+h/2,y(k,:)+f3= %得中f4=h*feval(fun,x(k)+h,y(k,:)+f4= %得中y(k1y(kf12*(f2f3)+ %得12-7】龍格--y'yy(x0)0,x0y1

clear x=x0+[0:Num]*h;a=1;yt1-exp(- funinline('-y %用inline構(gòu)造函數(shù)y0 PointNum [x1,ye]=MyEuler(fun,x0,xt,y0,PointNum);[x2,yh]=MyEulerPro(fun,x0,xt,y0,PointNum);[x3,yr]=MyRunge_Kutta(fun,x0,xt,y0,PointNum);plot(x,yt,'k',x1,ye,'b:',x2,yh,'g:',x3,yr,'r:')holdonplot(x1,ye,'bo',x2,yh,'b+',PointNum= %估計(jì)各算法運(yùn)行PointNum步迭代的時(shí)tic%[x1,ye]=time_Euler=toc [x1,yh]=time_EulerPro=toc [x1,yr]=time_RK4 %得到Runge_Kutta>>time_Euler=time_EulerPro=time_RK412-412-4三種不同方法所得結(jié)程序運(yùn)行時(shí)間較長(zhǎng)，幾乎是簡(jiǎn)單歐拉法的五倍（0.2544s，而龍格-庫(kù)塔法運(yùn)行時(shí)間是1.0195s前面講到，求解器solver中的ode23采用二階、三階龍格-庫(kù)塔法；ode45采用四階、12-8】solver中的龍格-1。分別采用二、三、四、y'3y22x23x,且0剟 y(x)1,x 編寫(xiě)程序文件ex1208a.m：%ex1208a.m用ode23得到微分方程解并計(jì)算出該算法運(yùn)行時(shí)fun=inline('-3*y^2+2*x.^2+3*x','x','y'); %用inline構(gòu)造函數(shù)f(x,y) %可得到x,y輸出向量值ode23(fun,[0,1],1),hold 12-5a圖12- 二、三階龍格-庫(kù)塔法解微分方編寫(xiě)程序文件ex1208b.m：%ex1208b.m用ode45得到微分方程解,fun=inline('-3*y^2+2*x.^2+3*x','x','y'); %用inline構(gòu)造函數(shù)f(x,y)ode45(fun,[0,1],1),holdon t1、t2>>ex1208bt1=0.0256t212-5b12-5b四階、五階龍格-庫(kù)塔算法解微分12-9】solver中的龍格-2ode45求解如下二y'(2)y(2)[0.50.02x0y(125y(22解：ode45functiondx=myfun02(x,y)ex1209.m%[x,y]=ode45('myfun02',[015],[25 %畫(huà)出y(1)，y(2)12-612-612-9方程的12.4預(yù)估- 用的形式，它們分別是ABM（Adams-Bashforth-Moulton）法和Hamming法。ABMABM法與上面幾種方法的主要不同是：它屬于多步算法，yn1的值是由(xn,yn)，(xn1,yn1)(xn2,yn2)(xn3,yn3這幾個(gè)點(diǎn)通過(guò)預(yù)估-校正過(guò)程求得的，要分以下兩步f(xyf值表示。從而可得yn1的預(yù)估表達(dá)式pn1：pn1yn

h(9fn337第二步：對(duì)pn1的預(yù)估值進(jìn)行校正，得到y(tǒng)n1的校正值cn1cn1yn

h(fn25fn1通常使用的是改進(jìn)的ABM算法。即在上述迭代式的基礎(chǔ)上，再加入更正，因此pn1yn

h(9fn337mn1pn1251(cnpncn1yn

h(fn25 cn119(cn1pn1)o(h5，其計(jì)算精度優(yōu)于前面介ode113函數(shù)就是該算法的程序?qū)崿F(xiàn)，由于迭代過(guò)程較復(fù)雜，此處不給出具體程序有的讀者可以參考o(jì)de113的實(shí)現(xiàn)程序其具體命令格式與上一節(jié)介紹的ode23和ode45相同。HammingHamming法是另一種形式的多步預(yù)估-校，其截?cái)嗾`差也是o(h5)，運(yùn)算效率Hamming預(yù)估-校正為pn1yn34h(2fn23mn1pn1112(cn 1[9ynyn23h(

yn1cn1

9(cn1在中編程實(shí)現(xiàn)的Hamming算法函數(shù)為：MyHamming。功能：以Hamming算法求解常微分方程。調(diào)用格式：[xy]=MyHamming(fx0,xty0,NKC)。其中，f為函數(shù)名；x0xty0x0N為自變量在[x0,xt]上所取的點(diǎn)數(shù)；KC為確定是否使用改正；xxyHamming算法 function用hamming方法解向量微分方程y’(t)%函數(shù)f(x,%函數(shù)在x0%x:所取的點(diǎn)的xifnargin<KC ifnargin<5|N<=N %最大迭代數(shù)默認(rèn)為ifnargin<y0 %初值默認(rèn)為y0 %使y0為行向hxt-x0)/(N- xt0=[x,y 用xx(1:3 %重新調(diào)整xfork=F(k-1feval(f,x(k),y(k計(jì)算f2f3p= 預(yù)估量初c= h34= %預(yù)估中系KC11= %更正中系KC91= %最后計(jì)算y值的中更正項(xiàng)的系h312=3*h*[-12 %校正中各f項(xiàng)系fork=4:N-p1y(k-3h34*(2*(F(1F(3F(2p(n+1)m1=p1+KC11*(c-p); c1=(-y(k-2,:)+9*y(k,:)+...h312*[F(2:3,:);feval(f,x(k+1),m1)])/8; y(k+1,:)=c1-KC91*(c1- p c F=[F(2:3,:);feval(f,x(k+1),y(k+【例12-10】預(yù)估-校求解常微分方程應(yīng)用實(shí)例。采用ABM法與Hamming法解y'y1y(x)0,x 解：在中編寫(xiě)ex1210.m來(lái)實(shí)現(xiàn)求解%ex1210.m用ode113(),MyHamming(),解微分方程y'=-y+1clear,clfx0 xt y0 N= fun66inline('- f66=inline('1-exp(-t)','t'); [x113,y113]=ode113(fun66,[x0,xt],y0);t113=toc[x1,yH]=MyHamming(fun66,x0,xt,y0,N);tH=tocyt1= 作legend('精確解','Hamming解')title('Hamming法所得的解')subplot(122plot(x113,yt113,'*rx113,y113,'og')legend('精確解','ABM解')title('ode113即ABM算法所得的解12-712-7ABM法與Hamming法解微分12.5常微分方程求解綜合實(shí) 3d3 d2

dx32dx23dx

dy1 dy2dy3

x3y3x3y2 y1

0

d dxdxy2 1Y

3e2x

x3

x3方程變成上述形式后，其求解過(guò)程與一階微分方程相同，不過(guò)值得注意的是，y的初function %初始化變量 %dy(1)表示y的一階導(dǎo)數(shù),其等于y %dy(2)表示y的二階導(dǎo)dy(3)=2*y(3)/x^3+3*y(2)/x^3+3*exp(2*x)/x^3dy(3)表示y的三階導(dǎo)MyEulerProMyRunge_Kutta%ex1211.m用自編函數(shù)MyEulerPro(),MyRunge_Kutta()求多階微分方程 y0=[11030]; 12-84

x

改進(jìn)歐拉法解R-K法解3210

12-8用兩種方法求多階微分方程的【例12-12】多階常微分方程求解實(shí)例2。 功能函數(shù)ode23、ode453d3 d2

dx32dx23dx

解： %ex1212用ode23ode45ode113解多階微分方[x23,y23]=ode23('myfun03',[1,10],[110[x45,y45]=ode45('myfun03',[1,10],[110[x113,y113]=ode113('myfun03',[1,10],[110 legend('ode23解','ode45解','ode113解')title('ODE函數(shù)求解結(jié)果') %以ode45xODEx3210

圖12-9a中ODE功能函數(shù)求解結(jié)xy,y一階導(dǎo)數(shù),y二階導(dǎo)數(shù)函數(shù)xyyyy一階導(dǎo)數(shù)y兩二導(dǎo)數(shù)876543210 12-9by及其一階、二階導(dǎo)數(shù)12.6差分方程用filter N y(n)brx(nr)aky(nr k在中，可以用函數(shù)filter()求解差分方程。函數(shù)filter有多種使用形式，常用y=filter(B,A,其中，A[a0, ,aN]，B[b0 ,bM]，x為輸入信號(hào)序列】y(n)0.7y(n1)0.6y(n2)y(n3)解：x(n(ny(n，利用filter函數(shù)編程求解，在中編寫(xiě)程序ex1213來(lái)實(shí)現(xiàn)。 n=[- %輸入脈沖等價(jià) 12-10脈沖響43210 12-10單位脈沖響應(yīng)12-14】2y(n)0.7y(n1)0.6y(n2)y(n3)x(n)0.5x(nx0,0,0,0,0,1,1,1,1,1,1,1,1,0,0,0,0,0,0，對(duì)應(yīng)于n的值為從-9到9，并繪制輸入解： 中編寫(xiě)程序文件ex1214.m來(lái)實(shí)現(xiàn)求解%ex1214.m求差分方程的解n=[-yColumns1through