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

第九章方程与微分方程一、代数方程的求根多项式求根（1）matlab中的多项式是用行向量表示，如：p=[2,3,-1,4,5]

表示多项式：（2）求多项式在某点的值的函数：polyval(p,x)p：多项式对应的行向量,x：待求值的点如：求在x=1,2,-3,5点的值p=[2,3,-1,4,5]x=[1,2,-3,5]y=polyval(p,x)y=1365651625

（3）polyder(p)

求多项式p的导数。（4）poly(x)

构造以x为根的多项式。（5）m=roots(p)

求多项式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.求方程根符号解sovle函数格式：solve(f,x)功能：对变量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解一个方程：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.求方程根的数值解大多数非线性方程的求根问题非常困难，通常需要求方程的数值根。常用求数值解的方法有：二分法、不动点迭代法、牛顿迭代法、牛顿下山法、割线法等。（1）求方程数值解的函数fzero格式：fzero(f,x0)功能：求方程f=0在x0附近的一个根例9.1求方程函数sin(x+1)-0.2x=0的根。x=-5:0.05:5;y=sin(x+1)-0.2*x;plot(x,y)gridonpausen=input('输入根的个数n=')

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

end例9.2求解非线性方程组：（2）求非线性方程（组）数值解的函数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将求得的根标注在图形上，观察是否正确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’,‘初始条件,‘自变量’）初始条件省缺时，是求微分方程的通解。

Dy代表y的导数，D2y代表y的二阶导数，D3y代表y的三阶导数……例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')

结果：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.求微分方程数值解多数微分方程是没有解析解的，即使一些看上去形式非常简单的微分方程，如这个微分方程不能用初等函数及其积分来表达它的解!

从实际问题中抽象出来的微分方程，通常主要依靠数值解法来解决！微分方程有解的一个充分条件：微分方程初值问题：（1）若f(x,y)在区域上连续，且满足：则初值问题（1）在[a,b]上有唯一解y=y(x)简单的说微分方程数值解法，就是将连续问题离散化，用数值方法给出它的解在某些离散节点上的近似值。求数值解的方法有欧拉法、龙格-库塔法等。

基于龙格-库塔法，matlab求微分方程数值解函数：[x,y]=ode23(‘函数’,‘求解区间’,‘初始值’)其它求解命令：ode45,ode15s,ode23s等

（采用2,3阶龙格-库塔法）x:自变量值，y:函数值注意：matlab求微分方程数值解时，高阶微分方程必须等价的变换成一阶微分方程组。例9.5求解微分方程初值问题：f=@(x,y)1./(1+x.^2)-2*y.^2ode23(f,[0,6],0)figure(2)ode45(f,[0,6],0)例9.6求解微分方程初值问题：解:设y1=y,问题化为：建立函数文件：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))求微分方程(组)数的值解时，需先定义表示方程(组)的函数文件，如：dy=ffff(x,y)然后调用：[x,y]=ode23(‘ffff’,[a,b],x0)

在求解n个方程的微分方程组时，y与x0均为n维向量，函数文件ffff中的待解方程组应以y的分量形式写出。

例9.7导弹追踪问题设位于坐标原点的甲舰向位于x轴上点(1,0)处的乙舰发射导弹，导弹头始终对准乙舰，若乙舰以最大速度v0（常数）沿平行于y轴的直线行驶，导弹速度是5v0

，求导弹运行的曲线；又乙舰行驶多远时会被击中？解：设在t时刻导弹位于P(x,y)处，乙舰位于Q(1,v0t)处，导弹运动曲线：y=f(x)由导弹头始终对准乙舰：由导弹速度是5v0：弧OP长=5|AQ|解:设y1=y,问题化为：

定义函数文件

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放射性废料的处理有一个时期美国原子能机构处理核废料时，一直采用将它们装入密封的圆桶里，然后扔到水深约91米海底的方法。对此科学家们表示担心，怕圆桶下沉到海底时，与海底碰撞发生破裂而造成核污染。为此工程师们进行碰撞试验，发现当圆桶下沉速度超过12.2米/秒与海底碰撞时，圆桶就可能发生破裂。按此试验结论，分析这种处理核废料的方式是否安全。已知圆桶重量239.456千克，体积为0.208立方米，海水浮力1025.94千克/立方米。（大量试验表明圆桶下沉时的阻力与下沉速度成正比，比例系数为0.12）分析：求出当圆桶下沉到海底时，下沉速度是否大于12.2米/秒重量W=239.456千克，体积V=0.208立方米，海水浮力B=1025.94V=213.396千克，阻力系数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求圆桶沉到海底时的速度：B=1025.94*0.208;k=0.12;W=239.456;m=W/9.8;tt=13.2042;symstdyy=di