MATLAB应用求解非线性方程.docx
- 文档编号:27904140
- 上传时间:2023-07-06
- 格式:DOCX
- 页数:24
- 大小:250.30KB
MATLAB应用求解非线性方程.docx
《MATLAB应用求解非线性方程.docx》由会员分享,可在线阅读,更多相关《MATLAB应用求解非线性方程.docx(24页珍藏版)》请在冰豆网上搜索。
MATLAB应用求解非线性方程
第7章求解非线性方程
7.1多项式运算在MATLAB中的实现
一、多项式的表达
n次多项式表达为:
,是n+1项之和
在MATLAB中,n次多项式可以用n次多项式系数构成的长度为n+1的行向量表示
[a0,a1,……an-1,an]
二、多项式的加减运算
设有两个多项式
和
。
它们的加减运算实际上就是它们的对应系数的加减运算。
当它们的次数相同时,可以直接对多项式的系数向量进行加减运算。
当它们的次数不同时,应该把次数低的多项式无高次项部分用0系数表示。
例2计算
a=[1,-2,5,3];b=[0,0,6,-1];c=a+b
例3设
,
,求f(x)+g(x)
f=[3,-5,2,-7,5,6];g=[3,5,-3];g1=[0,0,0,g];%为了和f的次数找齐
f+g1,f-g1
三、多项式的乘法运算
conv(p1,p2)
例4在上例中,求f(x)*g(x)
f=[3,-5,2,-7,5,6];g=[3,5,-3];
conv(f,g)
四、多项式的除法运算
[Q,r]=deconv(p1,p2)
表示p1除以p2,给出商式Q(x),余式r(x)。
Q,和r仍为多项式系数向量
例4在上例中,求f(x)/g(x)
f=[3,-5,2,-7,5,6];g=[3,5,-3];
[Q,r]=deconv(f,g)
五、多项式的导函数
p=polyder(P):
求多项式P的导函数
p=polyder(P,Q):
求P·Q的导函数
[p,q]=polyder(P,Q):
求P/Q的导函数,导函数的分子存入p,分母存入q。
参数P,Q是多项式的向量表示,p,q也是多项式的向量表示。
例4求有理分式
的导函数
P=[3,5,0,-8,1,-5];%有理分式分子
Q=[10,5,0,0,6,0,0,7,-1,0,-100];%有理分式分母
[p,q]=polyder(P,Q)
六、多项式求根
多项式求根就是求满足多项式p(x)=0的x值。
N次多项式应该有n个根。
这些根可能是实根,也可能是若干对共轭复根。
其调用格式是
x=roots(P)
其中P为多项式的系数向量,求得的根赋给向量x,即x
(1),x
(2),…,x(n)分别代表多项式的n个根。
该命令每次只能求一个一元多项式的根,该指令不能用于求方程组的解,必须把多项式方程变成Pn(x)=0的形式;
例4求方程
的解。
首先将方程变成Pn(x)=0的形式:
roots([1-10-1])
例5求多项式x4+8x3-10的根。
A=[1,8,0,0,-10];
x=roots(A)
若已知多项式的全部根,则可以用poly函数建立起该多项式,其调用格式为:
P=poly(x)
若x为具有n个元素的向量,则poly(x)建立以x为其根的多项式,且将该多项式的系数赋给向量P。
例6已知f(x)=3x5+4x3-5x2-7.2x+5
(1)计算f(x)=0的全部根。
(2)由方程f(x)=0的根构造一个多项式g(x),并与f(x)进行对比。
P=[3,0,4,-5,-7.2,5];
X=roots(P)%求方程f(x)=0的根
G=poly(X)%求多项式g(x)
将这个结果乘以3,就与f(x)一致
7.2求解非线性方程f(x)=0
方程求根的一般形式是求下列方程的根:
f(x)=0(l)
实际上,就是寻找使函数f(x)等于零的变量x,所以求方程(l)的根,也叫求函数f(x)的零点。
如果变量x是列阵,则方程(l)就代表方程组。
当方程(l)中的函数f(x)是有限个指数、对数、三角、反三角或幂函数的组合时,则方程(l)被称为超越方程,例如e-x-sin(πx/2)+lnx=0就是超越方程。
当方程(l)中的函数f(x)是多项式时,即f(x)=Pn(x)=anxn+an-1xn+…+alx+a0,则方程(l)就成为下面的多项式方程,也称代数方程:
Pn(x)=anxn+an-1xn+…+alx+a0=0
(2)
Pn(x)的最高次数n等于2、3时,用代数方法可以求出方程
(2)的解析解,但是,当n≥5时,伽罗瓦(Galois)定理已经证明它是没有代数求根方法的。
至于超越方程,通常很难求出其解析解。
所以,方程(l)的求解经常使用作图法或数值法,而计算机的发展和普及又为这些方法提供了广阔的发展前景,使之成为科学和工程中最实用的方法之一。
本章首先介绍求解f(x)=0的MATLAB符号法指令,然后介绍求方程数值解的基本原理,最后再介绍求解f(x)=0的MATLAB数值法指令。
一、符号方程求解
在MATLAB中,求解用符号表达式表示的代数方程可由函数solve实现,其调用格式为:
solve(s):
求解符号表达式s的代数方程,求解变量为默认变量。
当方程右端为0时,方程可以不标出等号和0,仅标出方程的左端。
solve(s,v):
求解符号表达式s的代数方程,求解变量为v。
solve(s1,s2,…,sn,v1,v2,…,vn):
求解符号表达式s1,s2,…,sn组成的代数方程组,求解变量分别v1,v2,…,vn。
例1.解下列方程。
1.
x=solve('1/(x+2)+4*x/(x^2-4)=1+2/(x-2)','x')
2.
f=sym('x-(x^3-4*x-7)^(1/3)=1')
x=solve(f)
3.
x=solve('2*sin(3*x-pi/4)=1')
4.
x=solve('x+x*exp(x)-10','x')%仅标出方程的左端
二、求方程f(x)=0数值解的基本方法
并非所有的方程f(x)=0都能求出精确解或解析解,不存在这种解的方程就需要用数值解法求出近似解,有几种常见的数值解法基本原理:
二分法。
1求实根的二分法原理
设方程f(x)=0中的函数f(x)为实函数,且满足:
①函数f(x)在[a,b]上单调、连续;
②方程f(x)=0在(a,b)内只有一个实根x*。
则求方程f(x)=0的根,就是在(a,b)内找出使f(x)为零的点x*:
f(x*)=0,即求函数f(x)的零点。
因为f(x)单调连续,由连续函数的性质可知,若任意两点aj,bj[a,b],而且满足条件f(aj)f(bj)<0,则闭区间[aj,bj]上必然存在方程的根x*,即x*[aj,bj]。
据此原理提出求实根的二分法如下图所示,
图1方程求根二分法原理示意图
先用中点
将区间[a,b]平分为两个子区间(a,b1)和(b1,b),方程的根必然在子区间两端点上函数值之积小于零的那一半中,即不在(a,b1)内,就在(b1,b)内,除非f(b1)=0,于是寻根的范围缩小了一半。
图1中的根x*在区间中点左侧,即x*(a,bl)。
再将新的含根区间(a,b1)分成两半,重复上述步骤确定出更新的含根子区间。
如此重复n次,设含根区间缩小为(an,bn),则方程的根x*(an,bn),这一系列含根的子区间满足:
(a,b)D(al,bl)(a2,b2)…(a0,b0)…
由于含根区间范围每次减半,子区间的宽度为
(n=1,2,….),显然当n→时,(bn一an)→0,即子区间收敛于一点x*,这个点就是方程的根。
若n为有限整数,取最后一个子区间的中点
作为方程根的近似值,它满足f(xn)≈0,于是有:
这就是近似值xn的绝对误差限。
假定预先要求的误差为,由
便可以求出满足误差要求的最小等分次数n。
下面是二分法的程序
function[c,err,yc]=bisect(f,a,b,delta)
%Input-fisthefunctioninputasastring‘f’
%-aandbaretheleftandrightendpoints
%.-deltaisthetolerance
%Output-cisthezero
%-yc=f(c)
%-erristheerrorestimateforc
ya=feval(f,a);
yb=feval(f,b);
ifya*yb>0,break,end%表示无解,结束
maxl=l+round((log(b-a)-log(delta))/log
(2));%从误差表达式得到最小等分次数n
fork=1:
max1
c=(a+b)/2;%取区间中点
yc=feval(f,c);
ifyc==0
a=c;
b=c;%这时解已经找到
elseifyb*yc>0
b=c;%区间减半
yb=yc;
else
a=c;
ya=yc;
end
ifb-a end c=(a+b)/2; err=abs(b-a); yc=feval(f,c) 2迭代法 迭代法是计算数学中的一种重要方法,用途很广,求解线性方程组和矩阵特征值时也要用到它。 这里结合非线性方程的迭代法求解,介绍一下它的基本原理。 迭代法基本原理 迭代法的基本原理就是构造一个迭代公式,反复用它得出一个逐次逼近方程根的数列,数列中每个元素都是方程根的近似值,只是精度不同。 迭代法求解方程 f(x)=0 (1) 时,先把方程等价地变换成形式 f(x)=x-g(x)=0, (2) 移项得出: x=g(x)(3) 若函数g(x)连续,则称(3)为迭代函数。 用它构造出迭代公式: xk+1=g(xk),k=0,l,2,…(4) 从初始值x0出发,便可得出迭代序列: {xk}=x0,x1,x2,….xk,…..(5) 如果迭代序列(5)收敛,且收敛于x*,则由式(4)有: 可见x*便是方程(l)的根。 迭代法几何意义: 如下图所示,解方程f(x)=0可以等价地变换成求解x=g(x), 图4-2方程求根迭代法原理示意图 在几何上,就等价求曲线y=x和y=g(x)交点P*的坐标x*。 求迭代序列(5),就等于从图中x0点出发,由函数y=g(x0)得出y=P0,代入函数y=x中得出Q1,再把Q1的x坐标x1代入方程y=g(x)得出P1,如此继续下去,便可在曲线y=g(x)上得到一系列的点P0,P1,…,Pk,…,这些点的x坐标便是迭代数列xl,x2,…,xk,…,它趋向于方程(l)的根x*,数列的元素就是方程根的近似值。 数列的收敛就等价于曲线y=x和y=g(x)能够相交于一点。 迭代公式收敛定理 要想用迭代法求出方程根的近似值,迭代序列(4-5)必须收敛。 下面的定理给出了迭代法的收敛条件,同时也给出了迭代公式的误差。 收敛定理: 方程x=g(x)在(a,b)内有根x*,如果: ①当x[a,b]时,g(x)[a,b]; ②g(x)可导,且存在正数q<1,使得对于任意x[a,b]都有|g’(x)|q<1,则有以下结论。 ①方程x=g(x)在(a,b)内有唯一的根x*。 ②迭代公式xk+1=g(xk)对(a,b)内任意初始近似根x0均收敛于x*。 ③近似根xk的误差估计公式为: (4-6) 3切线法 切线法就是从函数曲线上的一点出发,不断用曲线的切线代替曲线,求得收敛于根的数列。 切线法原理: 解非线性方程f(x)=0的切线法也称牛顿法,它是把方程线性化的一种近似方法,用函数f(x)的切线代替曲线产生一个收敛于方程根的迭代序列,从而得到方程的近似根。 把函数f(x)在某一初始值x。 点附近展开成泰勒级数: (4-7) 取其线性部分,近似地代替函数f(x)可得方程的近似式: 设 ,解该近似方程可得: 把函数f(x)在xl点附近展开成泰勒级数,取其线性部分替代函数f(x),设 ,得: 如此继续做下去,就可以得到牛顿迭代公式: ,k=0,l,2,…(8) 由式(8)得出的迭代序列xl,x2,…,xk…,在一定的条件下收敛于方程的根x*。 2.几何意义 图4一3方程求根切线法原理示意图 选取初值x0后,过 点作曲线 的切线,其方程为 。 设切线与X釉的交点为x1,则 ,再过 作切线,与x轴的交点为 ,如此不断作切线,求与x轴的交点,便可得出的一系列的交点x1,x2,…,xk,…,它们逐渐逼近方程的根x*。 3.切线法的收敛性 理论可以证明,在有根区间[a,b]上,如果 、 连续且不变号,则只要选取的初始近似根x0满足 f(x。 ),切线法必定收敛。 它的收敛速度经推导可得出: (9) 是个常数,式(9)表明用牛顿迭代公式在某次算得的误差,与上次误差的平方成正比,可见牛顿迭代公式的收敛速度很快。 4.2.4割线法(弦截法) 应用切线法的牛顿迭代公式时,每次都得计算导数 ,若将该导数用差商代替,就成为割线法(有时称快速弦截法)的迭代公式: ,k=0,l,2,…(4一10) 割线法的几何意义也很明显。 如图所示, 图4一4方程求根割线法原理示意图 过点(x0,f(x0))和(x1,f(x1))作函数y=f(x)曲线的割线,交X轴于点x2,再过点(x1,f(x1))和(x2,f(x2))作曲线的割线,交X轴于点x3,……一直做下去,则得到一系列割线与X轴的交点,这些交点序列将趋于方程的根x*。 非线性方程的数值解法还有许多,这里仅介绍了几种基本方法的原理。 二分法简单方便,但收敛速度慢; 迭代法虽然收敛速度稍微快点,但需要判断能否收敛; 只要初值选取得当,切线法具有恒收敛且收敛速度快的优点,但需要求出函数的导数; 弦截法不需要求导数,特别是前面介绍的快速弦截法,收敛速度很快,但是需要知道两个近似的初始根值才能作出弦,要求的初始条件较多。 这些方法各有千秋,需根据具体情况选用。 三、方程f(x)=0数值解的MATLAB实现 MATLAB中求方程数值解的办法很多,有的是专用指令,有的是根据方程性质而借用其他专用指令求得的。 4.3.2求函数零点指令fzero 求解方程f(x)=0的实数根也就是求函数f(x)的零点。 MATLAB中设有求函数f(x)零点的指令fzero,可用它来求方程的实数根。 该指令的使用格式为: fzero(fun,x0,options) ①输入参数fun为函数f(x)的字符表达式、内联函数名或M函数文件名。 ②输入参数x0为函数某个零点的大概位置(不要取零)或存在的区间[xi,xj],要求函数f(x)在x0点左右变号,即f(xi)f(xj)<0。 ③输入参数options可有多种选择,若用optimset('disp','iter')代替options时,将输出寻找零点的中间数据。 ④该指令无论对多项式函数还是超越函数都可以使用,但是每次只能求出函数的一个零点,因此在使用前需摸清函数零点数目和存在的大体范围。 为此,一般先用绘图指令plot,fplot或ezplot画出函数f(x)的曲线,从图上估计出函数零点的位置。 例4一5求方程x2+4sin(x)=25的实数根(-2π<x<2π)。 解: 一、fun为函数f(x)的字符表达式 (l)首先要确定方程实数根存在的大致范围。 为此,先将方程变成标准形式f(x)=x2+4sin(x)-25=0。 作f(x)的曲线图: x=-2*pi: 0.1: 2*pi; f=x.^2+4*sin(x)-25; plot(x,f);gridon; 从曲线上可以看出,函数的零点大约在x1≈-4和x2≈5附近。 (2)直接使用指令fzero求出方程在x1≈-4时的根。 x1=fzero('x^2+4*sin(x)-25',-4) 若键入: fzero('x^2+4*sin(x)-25',-4,optimset('disp','iter')),将显示迭代过程。 中间数据表明,求根过程中不断缩小探测范围,最后得出-4附近满足精度的近似根。 (3)求x2≈5的根: x2=fzero('x^2+4*sin(x)-25',5) 二、fun为函数f(x)的M函数文件名 将方程x2+4sin(x)=25编成M函数文件(实用中在函数较为复杂、而又多次重复调用时,才这样做),用fzero求解。 (1)在M文件编辑调试窗中键入: functionyy=li4_5(x) yy=x^2+4*sin(x)-25; 以li4_5为文件名存盘,退出编辑调试窗,回到指令窗。 (2)确定根的大体位置; (3)在指令窗中键入下述指令可求出-4附近的根: x1=fzero('li4_5',-4) 键入下述指令可求出5附近的根: x2=fzero('li4_5',5) 三、fun为函数f(x)的内联函数名 内联函数是MATLAB提供的一个对象(Object)。 它的性状表现和函数文件一样,但内联函数的创建比较容易。 inline('CE') 'CE'是字符串,CE为不包含赋值符号“=”的表达式。 上式把串表达式转化为输入宗量自动生成的内联函数。 上述调用格式将自动地对CE进行辫识,把CE中由字母/数字组成的连续字符认做变量,除“预定义变量名(如i,j,pi)”和“常用函数名(如sin)”以外的由字母/数字组成的连续字符将被认做变量。 但注意: 若连续字符后紧接“左圆括号”,那么将不被当作输入宗量。 如x (1),就不会认做输入宗量处理。 inline('CE',arg1,arg2,…) 上述调用格式把串表达式转化为arg1,arg2等指定输入宗量的内联函数;这种调用格式是创建内联函数的最稳妥、可靠途径。 输入宗量字符可表达得更自如。 将函数f(x)写成内联函数的形式: f=inline('x^2+4*sin(x)-25') 这时内联函数名为f 分别求x1≈-4和x2≈5时的根: x1=fzero(f,-4),x2=fzero(f,5) 例4-6求f(x)=x-10x+2=0在x0=0.5附近的根。 从f(x)的曲线看(x=-2.5: 0.01: 0.5;fx=x-10.^x+2;plot(x,fx)), 曲线的零点有两个,一个在x=-2附近,另一个在x=0.5附近 (1)建立函数文件funx.m。 (function[输出变量列表]=函数名(输入变量列表)) functionfx=funx(x) fx=x-10.^x+2; (2)调用fzero函数求根。 z=fzero('funx',0.5) 例2求 在 和 作为迭代初值时的零点。 从f(x)的曲线看(x=-7: 0.01: 2;f=x-1./x+5;plot(x,f)), 曲线的零点有两个,一个在x=-5附近,另一个在x=1附近 (1)建立函数文件fz.m。 functionf=fz(x) f=x-1./x+5; (2)调用fz.函数求根。 fzero('fz',-5)%以-5作为迭代初值 fzero('fz',1)%以1作为迭代初值 7.3求解非线性方程组数值解的迭代法 一、符号方程组求解 在MATLAB中,求解用符号表达式表示的方程组仍然可由函数solve实现,其调用格式与解用符号表达式表示的方程一样。 例解下列方程组。 1. [xy]=solve('1/x^3+1/y^3=28','1/x+1/y=4','x,y') 2. [uv]=solve('u^3+v^3=98','u+v=2','u,v') 3. [xy]=solve('x+y=98','x^(1/3)+y^(1/3)=2','x,y') 回车后出现下面的提示 Warning: Explicitsolutioncouldnotbefound. >InD: \MATLAB6p5\toolbox\symbolic\solve.matline136 如果做代换: ,方程3就变成方程2,就可解 这个问题说明,符号求解并不是万能的。 如果用MATLAB得出无解或未找到所期望的解时,应该用其它方法试探求解。 4. [xy]=solve('x^2+y^2=5','2*x^2-3*x*y-2*y^2')%变量由默认规则确定 二、求解非线性方程组的基本方法 对于非线性方程组(以二元方程组为例,其他可以类推) (11) 的数值解求法,跟一元非线性方程的切线法(牛顿法)雷同,也是把非线性函数线性化,近似替代原方程得出数值解,所以也叫作牛顿迭代法。 假设方程组(11)的初始估计值为(x0,y0),可以把方程组(11)中的两个函数f1(x,y)和f2(x,y)在(x0,y0)处用二元泰勒级数展开,只取线性部分,移项得出: (12) 若系数矩阵行列式 ,则方程组(12)的解为: , 方程组(11)中的两个函数f1(x,y)和f2(x,y)在(x1,y1)处,再用二元泰勒级数展开,只取线性部分,……如此继续替代下去,直到方程组的根达到所要求的精度,就完成了方程组的求解。 求解方程组(11)还有许多其他办法,如“最速下降法”,它是利用方程组(11)构成所谓模函数 ,通过求模函数极小值的方法得到方程组的数值解,诸如此类在此不再一一列举。 三、求方程组数值解的指令 fsolve是用最小二乘法求解非线性方程组F(X)=0的指令,变量X可以是向量或矩阵,方程组可以由代数方程或者超越方程构成。 它的使用格式为: fsolve('fun',X0,OPTIONS) ①参数fun是编辑并存盘的M函数文件的名称,可以用@代替单引号对它进行标识。 M函数文件主要内容是方程F(X)=0中的函数F(X),即方程左边的函数。 ②参数X0是向量或矩阵,为探索方程组解的起始点。 求解将从X0出发,逐渐趋向,最终得到满足精度要求、最接近X0的近似根X*: F(X*)≈0。 由于X0是向量或矩阵,无法用画图方法进行估计,实际问题中常常是根据专业知识、物理意义等进行估
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- MATLAB 应用 求解 非线性 方程
![提示](https://static.bdocx.com/images/bang_tan.gif)