代碼資源精通matlab科學計算12-3rd

12。常微分方程問題來近似求解中提供了求解函數(shù)和求解器可以實現(xiàn)常微分方程。通過本章學習，讀者能夠熟練使用的求解函數(shù)和求解器，以及通過編程進中的求解中的求解函的影響。當常微分方程式能夠解析求解時，可用的符號工具箱中的功能找到精確解在常微分方程難以獲得解析解的情況下使用的常微分方程求解器solver，符號解函數(shù)y'g(x,xyxy=f(x,y'=f'=g(x,y)y(x0y0下，可以得到唯一解。常微分方程符號解的語法是：表二階微分項y''，condition則為初始條件。dsolvedsolve('equation',dsolve('equation',vdsolve('equation','condition',v12-1】1。計算微分方程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】 常微分方程符號解求解實例2。計算微分方程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他指定點t0t1t

上的解，則令tspant0t1,t2

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

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

ifnargin<4 %y0默認值為0 %計算步長 y(1，：)=y0（：）’; fork=1:PointNum %計算f(x,y)在每個迭代點的y(k1,y(kh*f;%對于所取的點x迭代計算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%計算步長 for 12-212-2歐拉法所得方程解與符號解（精確解）對12-2中可以看出，用歐拉法得到的解和用符號法得到的解之間存在一定的誤差，差商和中心差商。有的讀者可以參考相關資料，編寫相應的程序，其基本方法，又稱為Henu算法。對方程y'f(x,y)兩邊由xn到xn1積分，并利用梯形，可得yn1ynh[f(xn,yn)f(xn1,yn1)]y(x0)具體迭代如下

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

由此可見，改進的歐拉法較之簡單歐拉法更為龍格-庫塔yn1ynh(k12k22k3k46k1f(xn,ynk2f(xnk3f(xn

1h,yn 1h,yn

hk2k4f(xnh,ynhk3在中編程實現(xiàn)的四階龍格-庫塔算法函數(shù)為：MyRunge_Kutta。調用格式：[xy]=MyRunge_Kutta(fun,x0,xty0,PointNumvarargin)。其中，fun為函數(shù)名；x0為自變量區(qū)間初值；xt為自變量區(qū)間終值；y0x0PointNum為自變量在[x0,xt]xxy四階龍格-庫塔法的程序代碼如下function[x,y]=%Runge-Kutta方法解微分方程形為%函數(shù)%y0表示函數(shù)在x0%varargin為可選輸入項可傳適當參數(shù)給函數(shù)%x:所取的點的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構造函數(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= %估計各算法運行PointNum步迭代的時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三種不同方法所得結程序運行時間較長，幾乎是簡單歐拉法的五倍（0.2544s，而龍格-庫塔法運行時間是1.0195s前面講到，求解器solver中的ode23采用二階、三階龍格-庫塔法；ode45采用四階、12-8】solver中的龍格-1。分別采用二、三、四、y'3y22x23x,且0剟 y(x)1,x 編寫程序文件ex1208a.m：%ex1208a.m用ode23得到微分方程解并計算出該算法運行時fun=inline('-3*y^2+2*x.^2+3*x','x','y'); %用inline構造函數(shù)f(x,y) %可得到x,y輸出向量值ode23(fun,[0,1],1),hold 12-5a圖12- 二、三階龍格-庫塔法解微分方編寫程序文件ex1208b.m：%ex1208b.m用ode45得到微分方程解,fun=inline('-3*y^2+2*x.^2+3*x','x','y'); %用inline構造函數(shù)f(x,y)ode45(fun,[0,1],1),holdon t1、t2>>ex1208bt1=0.0256t212-5b12-5b四階、五階龍格-庫塔算法解微分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 %畫出y(1)，y(2)12-612-612-9方程的12.4預估- 用的形式，它們分別是ABM（Adams-Bashforth-Moulton）法和Hamming法。ABMABM法與上面幾種方法的主要不同是：它屬于多步算法，yn1的值是由(xn,yn)，(xn1,yn1)(xn2,yn2)(xn3,yn3這幾個點通過預估-校正過程求得的，要分以下兩步f(xyf值表示。從而可得yn1的預估表達式pn1：pn1yn

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

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

h(9fn337mn1pn1251(cnpncn1yn

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

yn1cn1

9(cn1在中編程實現(xiàn)的Hamming算法函數(shù)為：MyHamming。功能：以Hamming算法求解常微分方程。調用格式：[xy]=MyHamming(fx0,xty0,NKC)。其中，f為函數(shù)名；x0xty0x0N為自變量在[x0,xt]上所取的點數(shù)；KC為確定是否使用改正；xxyHamming算法 function用hamming方法解向量微分方程y’(t)%函數(shù)f(x,%函數(shù)在x0%x:所取的點的xifnargin<KC ifnargin<5|N<=N %最大迭代數(shù)默認為ifnargin<y0 %初值默認為y0 %使y0為行向hxt-x0)/(N- xt0=[x,y 用xx(1:3 %重新調整xfork=F(k-1feval(f,x(k),y(k計算f2f3p= 預估量初c= h34= %預估中系KC11= %更正中系KC91= %最后計算y值的中更正項的系h312=3*h*[-12 %校正中各f項系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】預估-校求解常微分方程應用實例。采用ABM法與Hamming法解y'y1y(x)0,x 解：在中編寫ex1210.m來實現(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常微分方程求解綜合實 3d3 d2

dx32dx23dx

dy1 dy2dy3

x3y3x3y2 y1

0

d dxdxy2 1Y

3e2x

x3

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

x

改進歐拉法解R-K法解3210

12-8用兩種方法求多階微分方程的【例12-12】多階常微分方程求解實例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ù)求解結果') %以ode45xODEx3210

圖12-9a中ODE功能函數(shù)求解結xy,y一階導數(shù),y二階導數(shù)函數(shù)xyyyy一階導數(shù)y兩二導數(shù)876543210 12-9by及其一階、二階導數(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為輸入信號序列】y(n)0.7y(n1)0.6y(n2)y(n3)解：x(n(ny(n，利用filter函數(shù)編程求解，在中編寫程序ex1213來實現(xiàn)。 n=[- %輸入脈沖等價 12-10脈沖響43210 12-10單位脈沖響應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，對應于n的值為從-9到9，并繪制輸入解： 中編寫程序文件ex1214.m來實現(xiàn)求解%ex1214.m求差分方程的解n=[-yColumns1through