MATLAB教程：第九章 方程與微分方程

第九章方程與微分方程一、代數(shù)方程的求根多項(xiàng)式求根（1）matlab中的多項(xiàng)式是用行向量表示，如：p=[2,3,-1,4,5]

表示多項(xiàng)式：（2）求多項(xiàng)式在某點(diǎn)的值的函數(shù)：polyval(p,x)p：多項(xiàng)式對(duì)應(yīng)的行向量,x：待求值的點(diǎn)如：求在x=1,2,-3,5點(diǎn)的值p=[2,3,-1,4,5]x=[1,2,-3,5]y=polyval(p,x)y=1365651625

（3）polyder(p)

求多項(xiàng)式p的導(dǎo)數(shù)。（4）poly(x)

構(gòu)造以x為根的多項(xiàng)式。（5）m=roots(p)

求多項(xiàng)式p的的根如：x=[-1,2,3]q=poly(x)如：p=[-1,2,3,1]x=roots(p)q=1-416x=3.0796-0.5398+0.1826i-0.5398-0.1826i2.求方程根符號(hào)解sovle函數(shù)格式：solve(f,x)功能：對(duì)變量x解方程f=0symsxy='2*x^2-3*x-8'm=solve(y,x)m1=vpa(m,7)求方程2x2-3x-8=0的根m=3/4-73^(1/2)/473^(1/2)/4+3/4m1=-1.3860012.886001解一個(gè)方程：y=@(x)2*x^2-3*x-8m=solve(y)解方程組y1=@(x,y)2*x^2-3*y-8;y2=@(x,y)x-2*y+1;[m,n]=solve(y1,y2)m1=vpa(m,7)n1=vpa(n,7)m=313^(1/2)/8+3/83/8-313^(1/2)/8n=313^(1/2)/16+11/1611/16-313^(1/2)/16m1=2.586476-1.836476n1=1.793238-0.4182379

解方程組

symsxyy1='2*x^2-3*y-8';y2='x-2*y+1';f=solve(y1,y2)x=f.xy=f.yf=x:[2x1sym]y:[2x1sym]

x=3/8-313^(1/2)/8313^(1/2)/8+3/8

y=11/16-313^(1/2)/16313^(1/2)/16+11/163.求方程根的數(shù)值解大多數(shù)非線性方程的求根問(wèn)題非常困難，通常需要求方程的數(shù)值根。常用求數(shù)值解的方法有：二分法、不動(dòng)點(diǎn)迭代法、牛頓迭代法、牛頓下山法、割線法等。（1）求方程數(shù)值解的函數(shù)fzero格式：fzero(f,x0)功能：求方程f=0在x0附近的一個(gè)根例9.1求方程函數(shù)sin(x+1)-0.2x=0的根。x=-5:0.05:5;y=sin(x+1)-0.2*x;plot(x,y)gridonpausen=input('輸入根的個(gè)數(shù)n=')

fork=1:nb0=input('輸入初值b0=')a=fzero('sin(x+1)-0.2*x',b0)

end例9.2求解非線性方程組：（2）求非線性方程（組）數(shù)值解的函數(shù)fsolve格式：fsolve(f,x0)（x0為初值）h1=ezplot(@(x,y)5*x.^2+y.^2.*sin(y)-3)set(h1,'color',[1,0,0])holdonh2=ezplot(@(x,y)y.^2+10*cos(x.*y)+4*x+2*y-7);set(h2,'color',[0,0,1])holdoffgridon1,-10.5，-3

-0.5，-3-1，0-0.8，3-1.2,4①②③④⑤⑥

functionf=fun902(x)f1=5*x(1).^2+x(2).^2.*sin(x(2))-3;f2=x(2).^2+10*cos(x(1).*x(2))+4*x(1)+2*x(2)-7;

f=[f1;f2];clear,clcx=[1,0.5,-0.5,-1,-0.8,-1.2];y=[-1,-3,-3,0,3,4];m=length(x);r=[];

fork=1:mr(:,k)=fsolve(@fun902,[x(k);y(k)]);

endrr=0.93600.4483-0.3502-0.7746-0.7121-1.1380-1.2139-3.3232-3.35520.05183.09303.43971,-10.5，-3

-0.5，-3-1，0-0.8，3-1.2，4將求得的根標(biāo)注在圖形上，觀察是否正確h1=ezplot(@(x,y)5*x.^2+y.^2.*sin(y)-3)set(h1,'color',[1,0,0])holdonh2=ezplot(@(x,y)y.^2+10*cos(x.*y)+4*x+2*y-7);set(h2,'color',[0,0,1])r=[0.9360,0.4483,-0.3502,-0.7746,-0.7121,-1.1380;-1.2139,-3.3232,-3.3552,0.0518,3.0930,3.4397];u=r(1,:);v=r(2,:);plot(u,v,'mp')holdoffgridon二、微分方程求解1.求微分方程解析解

dsolve(‘方程1’,‘方程2’,…,‘方程n’,‘初始條件,‘自變量’）初始條件省缺時(shí)，是求微分方程的通解。

Dy代表y的導(dǎo)數(shù)，D2y代表y的二階導(dǎo)數(shù)，D3y代表y的三階導(dǎo)數(shù)……例9.3求解微分方程：y=dsolve('D2y-2*Dy+y-x^2=0','x')y=4*x+C48*exp(x)+x^2+C49*x*exp(x)+6y=dsolve('D2y+4*Dy+29*y','y(0)=0,Dy(0)=15','x')

結(jié)果：y=(3*sin(5*x))/exp(2*x)y=dsolve('D2y-1000*(1-y^2)*Dy-y=0','y(0)=2,Dy(0)=0','x')y=[emptysym]例9.4求微分方程組的通解[x,y,z]=dsolve('Dx=2*x-3*y+3*z','Dy=4*x-5*y+3*z','Dz=4*x-4*y+2*z')2.求微分方程數(shù)值解多數(shù)微分方程是沒(méi)有解析解的，即使一些看上去形式非常簡(jiǎn)單的微分方程，如這個(gè)微分方程不能用初等函數(shù)及其積分來(lái)表達(dá)它的解!

從實(shí)際問(wèn)題中抽象出來(lái)的微分方程，通常主要依靠數(shù)值解法來(lái)解決！微分方程有解的一個(gè)充分條件：微分方程初值問(wèn)題：（1）若f(x,y)在區(qū)域上連續(xù)，且滿足：則初值問(wèn)題（1）在[a,b]上有唯一解y=y(x)簡(jiǎn)單的說(shuō)微分方程數(shù)值解法，就是將連續(xù)問(wèn)題離散化，用數(shù)值方法給出它的解在某些離散節(jié)點(diǎn)上的近似值。求數(shù)值解的方法有歐拉法、龍格-庫(kù)塔法等。

基于龍格-庫(kù)塔法，matlab求微分方程數(shù)值解函數(shù)：[x,y]=ode23(‘函數(shù)’,‘求解區(qū)間’,‘初始值’)其它求解命令：ode45,ode15s,ode23s等

（采用2,3階龍格-庫(kù)塔法）x:自變量值，y:函數(shù)值注意：matlab求微分方程數(shù)值解時(shí)，高階微分方程必須等價(jià)的變換成一階微分方程組。例9.5求解微分方程初值問(wèn)題：f=@(x,y)1./(1+x.^2)-2*y.^2ode23(f,[0,6],0)figure(2)ode45(f,[0,6],0)例9.6求解微分方程初值問(wèn)題：解:設(shè)y1=y,問(wèn)題化為：建立函數(shù)文件：functiondy=fun906(x,y)dy=zeros(2,1);dy(1)=y(2);dy(2)=1000*(1-y(1)^2)*y(2)-y(1);[x,y]=ode15s('fun906',[0,3000],[2,0])plot(x,y(:,1))求微分方程(組)數(shù)的值解時(shí)，需先定義表示方程(組)的函數(shù)文件，如：dy=ffff(x,y)然后調(diào)用：[x,y]=ode23(‘ffff’,[a,b],x0)

在求解n個(gè)方程的微分方程組時(shí)，y與x0均為n維向量，函數(shù)文件ffff中的待解方程組應(yīng)以y的分量形式寫出。

例9.7導(dǎo)彈追蹤問(wèn)題設(shè)位于坐標(biāo)原點(diǎn)的甲艦向位于x軸上點(diǎn)(1,0)處的乙艦發(fā)射導(dǎo)彈，導(dǎo)彈頭始終對(duì)準(zhǔn)乙艦，若乙艦以最大速度v0（常數(shù)）沿平行于y軸的直線行駛，導(dǎo)彈速度是5v0

，求導(dǎo)彈運(yùn)行的曲線；又乙艦行駛多遠(yuǎn)時(shí)會(huì)被擊中？解：設(shè)在t時(shí)刻導(dǎo)彈位于P(x,y)處，乙艦位于Q(1,v0t)處，導(dǎo)彈運(yùn)動(dòng)曲線：y=f(x)由導(dǎo)彈頭始終對(duì)準(zhǔn)乙艦：由導(dǎo)彈速度是5v0：弧OP長(zhǎng)=5|AQ|解:設(shè)y1=y,問(wèn)題化為：

定義函數(shù)文件

functiondy=fun907(x,y)dy=zeros(2,1);dy(1)=y(2);dy(2)=1/5*sqrt(1+y(1)^2)/(1-x);[x,y]=ode15s('fun907',[0,0.99999],[0,0]);yy=y(:,1);plot(x,yy,'*')s=yy(end)得：s=0.2017例9.8放射性廢料的處理有一個(gè)時(shí)期美國(guó)原子能機(jī)構(gòu)處理核廢料時(shí)，一直采用將它們裝入密封的圓桶里，然后扔到水深約91米海底的方法。對(duì)此科學(xué)家們表示擔(dān)心，怕圓桶下沉到海底時(shí)，與海底碰撞發(fā)生破裂而造成核污染。為此工程師們進(jìn)行碰撞試驗(yàn)，發(fā)現(xiàn)當(dāng)圓桶下沉速度超過(guò)12.2米/秒與海底碰撞時(shí)，圓桶就可能發(fā)生破裂。按此試驗(yàn)結(jié)論，分析這種處理核廢料的方式是否安全。已知圓桶重量239.456千克，體積為0.208立方米，海水浮力1025.94千克/立方米。（大量試驗(yàn)表明圓桶下沉?xí)r的阻力與下沉速度成正比，比例系數(shù)為0.12）分析：求出當(dāng)圓桶下沉到海底時(shí)，下沉速度是否大于12.2米/秒重量W=239.456千克，體積V=0.208立方米，海水浮力B=1025.94V=213.396千克，阻力系數(shù)k=0.12，y=dsolve('m*D2y=W-B-k*Dy','y(0)=0','Dy(0)=0')y=(B*m-W*m-B*k*t+W*k*t)/k^2-(B*m-W*m)/(k^2*exp((k*t)/m))得：求解方程：B=1025.94*0.208;k=0.12;W=239.456;m=W/9.8;y=@(t)(B*m-W*m-B*k*t+W*k*t)/k^2-(B*m-W*m)./(k^2*exp((k*t)/m))-91;t=0:0.1:15;yy=y(t);plot(t,yy,t,0,'r.-')t0=input('輸入初值t0=')tt=fzero(y,t0)方程的解在t=13附近，鍵盤輸入t0=13得：tt=13.2042求圓桶沉到海底時(shí)的速度：B=1025.94*0.208;k=0.12;W=239.456;m=W/9.8;tt=13.2042;symstdyy=di