第六讲：Matlab多项式与代数方程求解器.ppt

1、Matlab 多项式运算与代数方程求解器,Matlab 多项式运算,Matlab 中多项式的表示方法,在 Matlab 中，n 次多项式是用一个长度为 n+1的向量来表示，缺少的幂次项系数为 0。,注：系数中的零不能省！ 按降幂顺序,例：,多项式的符号形式：poly2sym， p1=poly2str(p,x), poly2sym(2,-1,0,3),1.系数向量的直接输入法,例：输入多项式：,解：p=1 -5 6 -33; Poly2sym(p),2.由矩阵求其特征多项式，由函数poly实现。,a=1 2 3;2 3 4; 3 4 5; p1=poly(a) poly2sym(p1),由特征多

2、项式生成的多项式的首项系数一定为1。 N阶矩阵一般生成N次多项式。,3.由给定的根求其对应的多项式，也由poly函数实现。,若要生成实系数多项式，则根中的复数必定对应共轭。,root=-5 -3+4i -3-4i p=poly(root),syms a b c d v A=a b;c d; poly(A) %poly(A,v),4.多项式的符号表示,若矩阵A是符号表示其调用格式为：,Poly(A)%返回A的用符号x或t表示的特征多项式 Poly(A,v) )%返回A的用符号v表示的特征多项式 Poly2sym(C)%将向量C表示的多项式转化为用符号表示。其变量用符号x表示。 Poly2sym(

3、C,v) )%将向量C表示的多项式转化为用符号表示。其变量用符号v表示。,C=1 0 -2 -5; %poly2sym(C) poly2sym(C，t),多项式四则运算,多项式加减运算,Matlab 没有提供专门进行多项式加减运算的函数，事实上，多项式的加减就是其所对应的系数向量的加减运算,对于次数相同的多项式，可以直接对其系数向量进行加减运算； 如果两个多项式次数不同，则应该把低次多项式中系数不足的高次项用 0 补足，然后进行加减运算。,例：,多项式四则运算,多项式乘法运算：,k = conv(p,q),其中 k 返回的是多项式 p 除以 q 的商，r 是余式。,k,r=deconv(p,q

4、),p=conv(q,k)+r,p=2 -5 6 -1 9; poly2sym(p) d=3 -90 -18; poly2sym(d) pd=conv(p,d) poly2sym(pd) p1=deconv(pd,d),计算两多项式的乘除法,多项式的求导,polyder,k=polyder(p) : 多项式 p 的导数； k=polyder(p,q): p*q 的导数； k,d=polyder(p,q): p/q 的导数，k 是分子，d 是分母, k1=polyder(2,-1,0,3); k2=polyder(2,-1,0,3,2,1); k2,d=polyder(2,-1,0,3,2,1)

5、;,例：已知 ， ， 求,多项式的值,计算多项式在给定点的值,有两种形式，对应两种算法。,代数多项式求值,y = polyval(p,x): 计算多项式 p 在 x 点的值,注：若 x 是向量或矩阵，则采用数组运算 (点运算)！,多项式的值,矩阵多项式求值,Y=polyvalm(p,X),与代数多项式求值不同之处： 采用的是普通矩阵运算； X 必须是方阵,例：已知 ，则,polyvalm(p,A) = 2*A*A*A - A*A + 3*eye(size(A) polyval(P,A) = 2*A.*A.*A - A.*A + 3*ones(size(A), p=2,-1,0,3; x=-1,

6、 2;-2,1;polyval(p,x) polyvalm(p,x),多项式的零点（根）,x=roots(p)：若 p 是 n 次多项式，则输出是 p=0 的 n 个根组成的 n 维向量。,多项式的拟合,拟合函数：polyfit(X,Y,n) polyfit(X,Y,n) %X,Y为拟合的数据，n为拟合多项式的阶数。 p,s=polyfit(X,Y,n) % p为拟合多项式系数向量，s为拟合多项式系数向量的结构信息。,例：用5阶多项式对0，/2上的正弦函数值进行最小二乘法拟合,x=0:pi/20:pi/2;y=sin(x); a=polyfit(x,y,5); x1=0:pi/30:2*pi;

7、y1=sin(x1); y2=a(1)*x1.5+a(2)*x1.4+a(3)*x1.3+a(4)*x1.2+a(5)*x1+a(6); plot(x1,y1,b-,x1,y2,r*) legend(原曲线,拟合曲线) axis(0, 7,-1.2,4),多项式的部分分式展开,调用格式： r,p,k=residue(b,a) %多项式b(s)/a(s)的部分分式展开，r是留数，p是极点,k是直向量 b,a=residue(r,p,k),a=2 45 82;b=3 5 8 3; r,p,k=residue(b,a) bb,aa=residue(r,p,k),因式分解和展开,MATLAB提供了对符

8、号表达式、符号矩阵等进行因式分解、同类项合并和展开等操作。,因式分解的主要函数有：factor、collet、expand,多项式的简化,简化的主要函数有： pretty(s)%将代数式s转化为手写形式。 simplify%单一化简命令 simple,k = conv(p,q) k,r = deconv(p,q),k = polyder(p) k = polyder(p,q) k,d = polyder(p,q),y = polyval(p,x) Y = polyvalm(p,X),x = roots(p),多项式运算小结,多项式运算中， 使用的是多项式 系数向量， 不涉及符号计算！,poly

9、2sym(p),polyfit(X,Y,n),线性方程组求解,对于方程ax=b，a 为anm矩阵，有三种情况： 当n=m时，此方程成为“恰定”方程,寻求精确解 当nm时，此方程成为“超定”方程，寻求最小二乘解 当nm时，此方程成为“欠定”方程，寻求基本解 matlab定义的除运算可以很方便地解上 述三种方程,1.恰定方程组的解,方程ax=b(a为非奇异) x=a-1 b 矩阵逆 两种解: x=inv(a)b 采用求逆运算解方程 x=ab 采用左除运算解方程,方程ax=b a=1 2;2 3;b=8;13; x=inv(a)*b x=ab x = x = 2.00 2.00 3.00 3.00,

10、=,a x = b,例: x1+2x2=8 2x1+3x2=13,2.超定方程组的解,方程 ax=b ,mn时此时不存在唯一解。 方程解 (a a)x=a b x=(a a)-1 a b 求逆法 x=ab matlab用最小二乘法找一 个准确地基本解。,例: x1+2x2=1 2x1+3x2=2 3x1+4x2=3 a=1 2;2 3;3 4;b=1;2;3; 解1 x=ab 解2 x=inv(aa) a b x = x = 1.0000 1.0000 -0.0000 0.0000,=,a x = b,3.欠定方程组的解,当方程数少于未知量个数时,即不定 情况,有无穷多个解存在。 matlab

11、可求出两个解： 用除法求的解x是具有最多零元素的解 是具有最小长度或范数的解，这个解是基于伪逆pinv求得的。,x1+2x2+3x3=1 2x1+3x2+4x3=2 a=1 2 3;2 3 4;b=1;2; x=ab x=pinv(a)b x = x = 1.00 0.83 0 0.33 0 -0.17,a x = b,内部包含许多自适应算法：对超定方程用最小二乘法，对欠定方程它将给出范数最小的一个解。,求解下列方程组 0.4096x1+0.1234x2+0.3678x3+0.2943x4=0.4043 0.2246x1+0.3872x2+0.4015x3+0.1129x4=0.1550 0.

12、3645x1+0.1920 x2+0.3781x3+0.0643x4=0.4240 0.1784x1+0.4002x2+0.2786x3+0.3927x4=-0.2557,a=0.4096 0.1234 0.3678 0.2943 0.2246 0.3872 0.4015 0.1129 0.3645 0.1920 0.3781 0.0643 0.1784 0.4002 0.2786 0.3927; b=0.4043 0.1550 0.4240 -0.2557; x=ab,在MATLAB命令窗口输入,x = -0.1819 -1.6630 2.2172 -0.4467,线性方程组的解析解,lin

13、solve(A,b)：解线性方程组,例：解方程组, A=1 2 1; 1 0 1; 1 3 0; b=2;3;8; x=linsolve(A,b),b是列向量！,非线性方程的根,非线性方程的解法：,1.二分法,function y=erfen(fun,a,b,esp) if narginesp if feval(fun,a)*feval(fun,c)0 b=c;c=(a+b)/2; elseif feval(fun,c)*feval(fun,b)0 a=c;c=(a+b)/2; else y=c;esp=10000; end n=n+1; end y=c; elseif feval(fun,a

14、)=0 y=a; elseif feval(fun,b)=0 y=b; else disp(These may not be a root on the intercal); end n,用二分法求方程 的正根，要求误差小于0.05。,在命令窗口输入：erfen(fc,0,10,0.05),先编制函数文件fc.m： function y=fc(x) y=x2-x-1;,n = 56 ans = 1.6180,Matlab 非线性方程的数值求解,fzero(f,x0)：求方程 f=0 在 x0 附近的根。,方程可能有多个根，但 fzero 只给出距离 x0 最近的一个,fzero 先找出一个包含

15、 x0 的区间，使得 f 在这个区间两个端点上的函数值异号，然后再在这个区间内寻找方程 f=0 的根；如果找不到这样的区间，则返回 NaN。,x0 是一个标量，不能缺省,由于 fzero 是根据函数是否穿越横轴来决定零点，因此它无法确定函数曲线仅触及横轴但不穿越的零点，如 |sin(x)| 的所有零点。,非线性方程的根,fzero 的另外一种调用方式,fzero(f,a,b),方程在 a,b 内可能有多个根，但 fzero 只给出一个,求方程 f=0 在 a,b 区间内的根。,参数 f 可通过以下方式给出：,fzero(x3-3*x+1,2); f=inline(x3-3*x+1); fzer

16、o(f,2) fzero(x)x3-3*x+1,2);,f 不是方程！也不能使用符号表达式！,例：, fzero(sin(x),10) fzero(sin,10), fzero(x3-3*x+1,1) fzero(x3-3*x+1,1,2), fzero(x3-3*x+1=0,1),X, fzero(x3-3*x+1,-2,0) f=inline(x3-3*x+1); fzero(f,-2,0),Matlab 符号方程求解器,s=solve(f,v)：求方程关于指定自变量的解； s=solve(f)：求方程关于默认自变量的解。,f 可以是用字符串表示的方程，或符号表达式； 若 f 中不含等号，则表示解方程 f=0。,solve,例：解方程 x3-3*x+1=0, syms x; f=x3-3*x+1; s=solve(f,x), s=solve(x3-3*x+1,x), s=solve(x3-3*x+1=0,x),Matlab 符号方程求解器,solve 也可以用来解方程组,solve( f1 , f2 , . , fN , v1 , v2 , . , vN),求解由 f1 , f2 , . , fN 确定的方程组关于 v1 , v2 , . , vN 的解,例：解方程组, x,y,z=solve(x+2*y-z=27,x+z=3, .