MatLab讲义 2002年9月版
end d=d-1; end end
此时,调用char(p)将的如下结果: char(p)
ans=x^3-2*x-5
5.定义显示输出
定义显示输出的文件名为display。很多类的显示输出定义非常简单,仅仅是输出变量名,并将变量的值转换为字符(调用函数char)输出。现为类polynom定义一个简单的显示输出文件@polynom/display.m. function display(p) disp(' ');
disp([inputname(1),' =']); disp(' ');
disp([' ' char(p)]) disp(' ');
4.15 重载与继承
重载包括三个方面的内容:算法的重载,操作符的重载和函数的重载。重载,实际是建立一系列特定的M文件,以指定MATLAB在操作自定义类对象时的具体操作。如polynom对象加法的重载例子。有polynom对象p,q,当遇到p+q时,将自动在@polynom路径下寻找plus.m文件,如果文件存在,执行文件。要重载这个文件,实质就是建立这个文件。 function r=plus(p,q) p=polynom(p); q=polynom(q);
k=length(q,c)-length(p.c);
r=polynom([zeros(1,k) p.c]+[zeros(1,-k) q.c]);
第五章MATLAB基本应用领域
5.1 线性代数
主要进行线性方程的求解、矩阵求逆、LU分解、QR分解、矩阵求幂、矩阵指数函数、求特征值及奇异值分解。
1.MATLIB中的矩阵
矩阵定义为二维实或复阵列。矩阵的加、减、乘、除、转置运算是最基本的运算,它们应符合维数一致的限制。
2.线性代数方程求解 一般线性方程可表示成: AX=B XA=B
当矩阵A为方阵时,它的解为X=A\\B (X=inv(A)B)或X=B/A。
当矩阵A为非奇异时,线性方程的解唯一;当矩阵A为奇异时,线性方程的解要么不存在,要么不唯一。
a=[1 -2 -3 -4;2 1 -1 1;-1 0 -1 2;3 -3 4 -5]; b=[1 1 1 1]' x=a\\b
当矩阵A为(m*n)维矩阵,且m>n时,方程AX=B中,方程个数多于变量个数,因此应采用最小二乘
36
MatLab讲义 2002年9月版
法来求解。
例:(ex32.m)对一组测量数据:
t=[0 0.3 0.8 1.1 1.6 2.3]‘
y=[0.82 0.72 0.63 0.6 0.55 0.5]‘
拟用延迟指数函数来拟合这组数据:
y(t)≈c1+c2e-t
得6个方程,而未知数仅有c1、c2两个。应采用最小二乘法来求解。
A=[ones(size(t)) exp(-t)] c=A\\y
T=[0:.1:2.5]‘;
Y=[ones(size(T))exp(-T)]*c; plot(T,Y,‘-?,t,y,‘o‘)
3.矩阵求逆
det(A):函数可求矩阵A的行列式值。 inv(A): 函数可求矩阵A的逆矩阵值。 pinv(A):用于计算非方阵的伪逆。 4.LU、QR分解
通过高斯对消或LU分解法,可将任意方阵表示成一个下三角阵与一个上三角阵的乘积:A=LU。
a=[1:3;4:6;4 2 6 ] [l,u]=lu(a)
得到det(a)=det(l)*det(u)
inv(a)=inv(u)*inv(l)
5.矩阵求幂和矩阵指数
23
矩阵求幂如A,B可以用A^2,B^3 sqrtm(A):计算A的开方
A
expm(A):计算矩阵A的指数e 6.特征值与特征向量
矩阵A的特征值λ和特征矢量ν,满足Aν=λν。 eig(A)求矩阵A的特征值λ。 7.一元多项式的运算
(1)一元多项式的表示和创建
任意一个多项式都可以用一个行向量来表示,它的系数是按降序排列: 例:x^3-6*x^2+11*x-6 表示为 p=[1 –6 11 –6]
poly2sym(p,‘x‘) 则 ans=x^3-6*x^2+11*x-6 poly2sym(p,‘x‘)函数:多项式的符号表示。
poly()函数:若A是n*n 矩阵,则p=poly(A),p为A的特征多项式,是n+1维向量,第一个元素的值一定是1。而特征多项式的根就是A的特征值。
u=[r1,r2,r3,......rm]
poly(u)=(x-r1)(x-r2).....(x-rm)
例:A=[6 –11 6;1 0 0;0 1 0];p=poly(A);
p=[1.0 –6.0 11.0 –6.0]
roots(p) %%求多项式的根与 eig(A)相同。
例:求1,2,-3为根的多项式 u=[1 2 –3] p=poly(u)
poly2sym(p) %% x^3-7*x+6 (b)多项式的基本运算 ? 多项式的加减运算:
37
MatLab讲义 2002年9月版
MATLIB中没有多项式加减的函数,可以自编函数。 例:(qq.m)
function yy=qq(x,y)
nx=length(x); % x=reshape(x,1,nx) ny=length(y); %y=reshape(y,1,ny); n=max(nx,ny) cc=zeros(1,n); if nx>ny
cc(1,(nx-ny+1):nx)=y;yy=x+cc; elseif nx cc(1,(ny-nx+1):ny)=x;yy=y+cc; else yy=x+y; end 函数的调用为:p1=[3 4 6],p2=[5 2 –4 7] 则:求多项式的和为:p3=qq(p1,p2) 求多项式的差为:p3=qq(p1,-p2) ? 多项式的乘法运算 w=conv(u,v) 即求多项式u和v的乘积,即求向量u和v的卷积,公式如下: 若m=length(u);n=length(v);则w的长度为 m+n+1。 且w(k)=u(1)v(k)+u(2)v(k+1)+....u(k)v(2k-1) 如:p1=[3 4 6];p2=[5 2 –4 7]; p3=conv(p1,p2); poly2sym(p3,‘x‘) %显示多项式15*x^5+26*x^4+26*x^3+17*x^2+4*x+42 两个以上的多项式的乘法需要重复使用conv. ? 多项式的除法运算 [q,r]=deconv(v,u)表示多项式v除以多项式u得到的商q和余数r,q,r均为多项式,进行解卷积的运算。公式如下: 例:[q,r]=deconv(p3,p1); poly2sym(q,’x’) %%q=5*x^3+2*x^2-4*x+7 [q,r]=deconv(p3,p2); poly2sym(q,’x’) %%q=3*x^2+4*x+6 ? 求多项式的根 roots(c)返回多项式的根组成的向量,也是多项式的友元阵的特征值向量,系数为实数的多项式的根,若根为复数,则复数是成对出现的。 如roots(p); 求p=x^3-6x^2+15x-4的根 p=[1 –6 15 –4],x=roots(p) x= 2.8494 + 2.2726i 38 MatLab讲义 2002年9月版 2.8494 - 2.2726i 0.3011 ? 多项式的数组运算 按数组运算规则计算多项式的值polyval。 y=polyval(p,x)计算多项式在x处的值,x可以是矩阵或向量。 x=[1:3;4:6];p=[1 2 3 4]; polyval(p,x) C.多项式的因式分解 n阶多项式有n个根r1,r2,...rn, 分解因式为(x-r1)(x-r2)...(x-rn)有可能有复根,复根则合并,因此在实数范围内的分解后的因式的最高次数不会超过2,分解因式后返回n*3的矩阵,其中n为分解后的因子个数,每一行就代表一个因子。 例:多项式的因式分解的自编函数(expuly.m) function y=expuly(p) s=p(1);r=roots(p);y=[]; while ~isempty(r) c=r(1);r(1)=[]; if imag(c)>eps cc=conv([1 -c],[1 -conj(c)]); y=[y;cc]; r(find(abs(r-conj(c)) cc=real(c);cc=[0 1 -c];y=[y;cc]; end end if abs(s-1)>eps cc=[0 0 s];y=[y;cc]; end D.最大公约式和最小公倍式 最大公约式和最小公倍式是整系数多项式的两个比较重要的概念,MATLAB本身并没有提供来求多项式的最大公因式和最小公倍式的函数。自编函数gcd.m,利用循环相除法求最大公约式。最小公倍式利用两个多项式的乘积除以二者的最大公因式即可以得到最小公约式。 例:最大公约式和最小公倍式的自编函数(gcd.m) function y=gcd(f,g) while ~isempty(f) if abs(f(1))<2*eps f(1)=[]; else break; end end while ~isempty(g) if abs(g(1))<2*eps g(1)=[]; else 39 MatLab讲义 2002年9月版 break; end end if isempty(f)|isempty(g) error('Arguments should not be null'); return; end f=f/f(1);g=g/g(1); while 1 [q,r]=deconv(f,g); while ~isempty(r) if abs(r(1))<2*eps r(1)=[]; else break; end end if isempty(r) y=g/g(1); return else f=g; g=r; end if length(r)==1 y=1; return; end end E.其它: ▲poly函数 可以把向量转化为多项式,向量的元素为多项式的根。 u=[r1 r2 r3 ... rm] poly(u)为多项式 (x-r1)(x-r2)....(x-rm) 例:求1 2 –3为根的多项式 则u=[1 2 –3];p=poly(u) p=[1 0 –7 6] ▲residue函数 [r,l,k]=residue 用于多项式的部分分式展开。 其中p,q分别为分子,分母多项式,r为留数,l为极点,k为余项。 例:num=10*[1 4 5 6 7]; %分子多项式 den=poly([-2;-1;-0.5]); %分子多项式 [res,poles,ki]=residue(num,den) res = -6.6667 -60.0000 40