第四章 MATLAB程序设计数值运算基础.ppt

1、1,第四章 数值计算基础,2,第四章 数值计算基础,4.1 多项式 4.2求解线性方程组的 4.3 差分、梯度 4.4 插值和拟合等。 4.5 基本数学函数,3,4.1 多项式,Matlab用行向量表示多项式，行向量由多项式系数按降幂排列组成。例如，多项式 可以用它的长度为n+1的系数行向量表示： 注意：系数行向量中元素的排列顺序必须是从高次幂系数到低次幂系数，多项式中缺少的幂次要用0补齐。,4,1、系数矢量的直接输入法: 由于在MATLAB中的多项式是以向量形式储存的，因此，创建多项式的最简单的方法即为直接输入多项式的系数行向量，MATLAB自动将向量元素按降幂顺序分配给各系数值。为了查看方

2、便，可利用转换函数poly2sym，将多项式由系数行向量形式，转换为符号形式。,4.1.1 多项式的创建,5,4.1.1 多项式的创建,例4-1 输入系数矢量 ，创建多项式 poly2sym(1, -5 ,6, -33) ans= x3-5*x2+6*x-33,6,4.1.1 多项式的创建,2、通过矩阵的特征多项式来创建多项式 可以通过求矩阵的特征多项式，来创建多项式，由函数poly实现。调用格式为： p=poly(A)：求矩阵A的特征多项式系数，要求输入参数A是nn的方阵，输出参数p是包含n+1个元素的行向量，是A的特征多项式系数向量。,7,4.1.1 多项式的创建,例4-2 求矩阵的特征多

3、项式。 A=1 2 3;4,5,6;7,8,0; p=poly(A) p = 1.0000 -6.0000 -72.0000 -27.0000 poly2sym(p) ans = x3-6*x2-72*x-27,8,4.1.1 多项式的创建,3、由根矢量创建多项式 由给定的多项式方程的根矢量创建该多项式，用poly函数实现。调用格式为： p=poly(r)：返回一个行向量，该行向量是以r为根的多项式系数向量。,9,4.1.1 多项式的创建,例4-3 由根矢量-5 -3 4创建多项式。 r=-5 -3 4 p=poly(r) p = 1 4 -17 -60 poly2sym(p) ans = x

4、3+4*x2-17*x-60,10,4.1.1 多项式的创建,由给定根矢量创建多项式时应注意: （1）如果希望生成实系数多项式，则根矢量的复数根必须共轭成对。 （2）有时生成的多项式向量包含很小的虚部，可用real命令将其滤掉。,11,4.1.1 多项式的创建,例4-4 根据根矢量r=-1+2i,-1-2i,0.2创建多项式 r=-1+2i,-1-2i,0.2 r = -1.0000 + 2.0000i -1.0000 - 2.0000i 0.2000 p=poly(r) % 求多项式系数矢量 p = 1.0000 1.8000 4.6000 -1.0000 pr=real(p) % 取实部

5、pr = 1.0000 1.8000 4.6000 -1.0000 poly2sym(pr) ans = x3+9/5*x2+23/5*x-1,12,4.1.2 多项式运算,1、求多项式的值 求多项式的值可以有两种形式，对应两种算法： （1）在输入变量值代入多项式计算时，以数量或矩阵中每个元素为计算单元的，对应函数为polyval; （2）以矩阵为计算单元的，进行矩阵式运算来求得多项式的值，对应函数为polyvalm。,13,4.1.2 多项式运算,调用格式为： polyval(p,x)：求多项式p在x点的值，x可以是数量或矩阵，x是矩阵时，表示求多项式p在x中各元素的值。 polyvalm(

6、p,x)：求多项式p对于矩阵x的值，x可以是数量或矩阵。x如果是数量，求得的值与函数polyval相同，如果x是矩阵则必须是方阵。,14,4.1.2 多项式运算,例4-5 求多项式 在2，4，6，8处的值，对于矩阵 的值，及在矩阵中各元素处的值。,15,4.1.2 多项式运算,通过上例可以得出：设A为方阵，P代表多项式，则， polyval(P,A): A.3-5*A.2+8*ones(size(A) polyvalm(P,A): A3-5*A2+8*eye(size(A),16,4.1.2 多项式运算,2、求多项式的根 可以直接调用MATLAB的函数roots，求解多项式的所有根。 调用形式

7、为： R=roots(C) 其中输入参数C是多项式系数行向量，输出参数R是多项式的根，一般用列向量表示。,17,4.1.2 多项式运算,例4-6 求出多项式 的根 a=1 2 -5 -6; r=roots(a) r = 2.0000 -3.0000 -1.0000 poly(r) ans = 1.0000 2.0000 -5.0000 -6.0000,18,4.1.2 多项式运算,3、多项式的乘除法运算 多项式的乘法和除法实质就是多项式系数向量的卷积和解卷运算。 乘法：c=conv(a,b) 求多项式a和b的乘法，如果向量a的长度为m，b的长度为n，则c的长度为m+n-1。 多项式的除法用函数

8、deconv实现，此函数也是向量的卷积函数的逆函数。 b,r=deconv(c,a) 向量a对向量c进行解卷，得到商向量b和余量r。,19,4.1.2 多项式运算,例4-7 （1）求两多项式的乘积 （2）求上述结果被 除所得的结果。 a=3 -3 4 -1 2; b=2 -1 1 2; c=conv(a,b) c = 6 -9 14 -3 3 5 0 4 d,r=deconv(c,b) d = 3 -3 4 -1 2 r = 0 0 0 0 0 0 0 0,20,4.1.2 多项式运算,4、多项式的微积分 Matlab中多项式的微分函数为polyder，多项式的积分函数为polyint。两个函

9、数的调用格式为： polyder(a) 求系数行向量为a的多项式的微分。 polyint(a) 求系数行向量为a的多项式的积分。,21,4.1.2 多项式运算,例4-8 （1）求多项式 的微分。 （2）将上述结果求积分。 a=1 -5 3 -6 4 -10; d=polyder(a) d = 5 -20 9 -12 4 poly2sym(d) ans = 5*x4-20*x3+9*x2-12*x+4 polyint(d) ans = 1 -5 3 -6 4 0,22,5、多项式的部分分式展开: 对于多项式b(x)和不含重根的n阶多项式a(x)之比，有如下展开： 式中 称为极点(Poles),

10、称为留数(Residue)，k(x)称为直项(Direct term)。,4.1.2 多项式运算,23,假如a(x)有m重根p(j)=p(j+m-1)，则相应部分写成：,4.1.2 多项式运算,24,在MALAB中，两个多项式之比用部分分式展开的函数为residue，有两种调用方法： r,p,k=residue(b,a) 求多项式之比b(x)/a(x)的部分分式展开，输出参数r为留数，p为极点和k为直项。 b,a=residue(r,p,k) 从部分分式得出多项式表达式b(x)和a(x)的系数向量，结果为对于表达式分母的归一形式。,4.1.2 多项式运算,25,例4-9 两个多项式的比为， 求

11、部分分式展开，再用展开的结果转换回原来的两个多项式, a=2,3,-4 b=1,2 r,p,k=residue(b,a) r = 0.0548 0.4452 p = -2.3508 0.8508 k = , b,a=residue(r,p,k) b = 0.5000 1.0000 a = 1.0000 1.5000 -2.0000,26,4.1.2 多项式运算,6、多项式拟合 matlab中多项式拟合的函数为polyfit ，其采用最小二乘法对给定的数据进行多项式拟合，给出拟合的多项式系数。函数的调用方式为： p=polyfit(x,y,n) 应用最小二乘法求出n阶拟合多项式p(x)。即用p(

12、x)拟合y(x)。x、y为数据的横纵坐标向量，n为拟合多项式的阶数，输出参数p为多项式p(x)的系数向量。,27,4.1.2 多项式运算,例4-10 x=1:10,求y（x）的5阶拟合多项式。 x=1:10; y=sqrt(x)+3*cos(x); p=polyfit(x,y,5) p = 0.0074 -0.1737 1.3312 -3.3680 0.3459 4.5606 %把原始数据点和拟合曲线绘制在同一个坐标系中，原始数%据用“。”表%示(如图4-1)。 plot(x,y,o,x,polyval(p,x),-) % 绘图函数用法在第6章介绍,28,4.2 求解线性方程组,4.2.1 齐

13、次线性方程组的解法 对于齐次线性方程组AX=0而言，可以通过求系数矩阵A的秩来判断解的情况： 1、如果系数矩阵的秩=n（方程组中未知数的个数），则方程组只有零解。 2、如果系数矩阵的秩n，则方程组有无穷多解。 可以利用MATLAB函数null(A)，求它的一个基本解。,29,4.2.1 齐次线性方程组的解法,例4-11 用matlab 求解方程组 A=1 1 1 1 -3 -1 1;1 0 0 0 1 1 0;-2 0 0 -1 0 -1 -2;r=rank(A); % 求矩阵A的秩x=null(A, r )得到解为:x = -0.2555 0.0565 -0.3961 -0.3138 -0.

14、0215 0.7040 0.5428 0.0967 0.2218 -0.1603 -0.2941 0.7991 0.8915 0.0717 -0.0151 -0.2386 0.1752 0.4429 -0.2353 0.2039 0.0803 -0.4994 0.6314 0.1099 -0.2304 0.1573 0.0879 0.3781 x的列向量为Ax=0的一个基本解。,30,4.2.2 非齐次线性方程组的解法,对于非齐次线性方程组AX=b而言，则要根据系数矩阵A的秩和增广矩阵B=A b的秩和未知数个数n的关系，才能判断方程组AX=b的解的情况。 （1）如果系数矩阵的秩=增广矩阵的秩=

15、n，则方程组有唯一解。 （2）如果系数矩阵的秩=增广矩阵的秩n，则方程组有无穷多解。 （3）如果系数矩阵的秩增广矩阵的秩，则方程组无解。,31,4.2.2 非齐次线性方程组的解法,求非齐次线性方程组（A*X=b）的通解时，需要先判断方程组是否有解，若有解，再去求通解。 求非齐次线性方程组（A*X=b）的通解的步骤为： 第一步：判断AX=b是否有解，若有解则进行第二步； 第二步：求AX=b的一个特解； 第三步：求AX=0的通解； 第四步：AX=b的通解为：AX=0的通解加上AX=b的一个特解。,32,4.2.2 非齐次线性方程组的解法,用matlab求解时，求Ax=b对应的齐次方程组Ax=0的通

16、解，可以利用函数null； 求Ax=b的特解，根据方程组中方程的个数m和未知数的个数n，可以把方程组Ax=b分为：恰定方程组（m=n），超定方程组（mn），欠定方程组（mn, 超定方程组，可以尝试计算最小二乘解； （3）mn，欠定方程组，可以尝试计算含有最少m的基解。,33,1、恰定方程组的求特解 方程Ax=b（A为非奇异） x=A-1b 两种方法: x=inv(A)b 采用求逆运算解方程 x=Ab 采用左除运算解方程 若A为奇异矩阵，则Ab给出出错信息,4.2.2 非齐次线性方程组的解法,34,恰定方程组的求特解,例: x1+2x2=8 2x1+3x2=13,=,A * x = b,x=in

17、v(A)*b x=Ab x = x = 2.00 2.00 3.00 3.00,4.2.2 非齐次线性方程组的解法,35,2、超定方程组的求特解一般求最小二乘解 方程 Ax=b ,mn时。 方程解 (A A)x=A b x=(A A)-1 A b 求逆法 x=Ab matlab用最小二乘法找一 个准确地基本解。,4.2.2 非齐次线性方程组的解法,36,超定方程组的求特解,例: x1+2x2=1 2x1+3x2=2 3x1+4x2=3 解1 x=ab 解2 x=inv(aa) a b x = x = 1.00 1.00 0 0.00,a *x = b,=,4.2.2 非齐次线性方程组的解法,3

18、7,3、欠定方程组的求特解 当方程数少于未知量个数时（mn）, 有无穷多个解存在。 matlab可求出两个解： 用除法求的解x是具有最多零元素的解 基于伪逆pinv求得的是具有最小长度或范数的解。,4.2.2 非齐次线性方程组的解法,38,欠定方程组的求特解,x1+2x2+3x3=1 2x1+3x2+4x3=2 x=ab x=pinv(a)b x = x = 1.00 0.83 0 0.33 0 -0.17,4.2.2 非齐次线性方程组的解法,39,例4-12 求方程组的解。,在Matlab中建立M文件如下 clear all A=1 1 1 1 -3 -1 1;1 0 0 0 1 1 0;-

19、2 0 0 -1 0 -1 -2; b=1,0,1; %输入矩阵A,b m,n=size(A); R=rank(A); B=A b; Rr=rank(B); %format rat if R=Rr elseif R=Rr diff(A) ans = 1 1 1 B(:,:,1)=1,2;3,4; B(:,:,2)=5,6;7,8; diff(B,1,2),ans(:,:,1) = 1 1 ans(:,:,2) = 1 1 diff(B,1,3) ans = 4 4 4 4,42,4.3.2 数值梯度,在matlab中求数值梯度的函数是gradient，该函数的调用格式为： Fx,Fy,Fz=g

20、radient(F) 如果F是一维的，那么只返回F的一维数值梯度Fx即；如果F是二维的，则返回Fx和Fy即，分别对应矩阵的列方向和行方向，依次类推。 =gradient(F,h1,h2,) 指定n个方向上相邻点之间的间距，如果不指定，缺省值为1。,43,4.3.2 数值梯度,如：A=round(10*rand(3,4) % round为四舍五入函数 A = 2 9 5 6 7 9 9 8 4 6 8 7 Fx,Fy=gradient(A) Fx = 7.0000 1.5000 -1.5000 1.0000 2.0000 1.0000 -0.5000 -1.0000 2.0000 2.0000

21、0.5000 -1.0000 Fy = 5.0000 0 4.0000 2.0000 1.0000 -1.5000 1.5000 0.5000 -3.0000 -3.0000 -1.0000 -1.0000,44,4.4 插值和拟合,4.4.1 插值 1、一维插值 一维插值就是对一维函数y=f(x)的数据进行插值，是最常用的插值运算，一维插值函数是interp1。 yi=interp1(x,y,xi,method) 输入参数x为原始数据点的横坐标向量，y为纵坐标向量或矩阵，method为插值方法选项。如果y是矩阵，那么插值按照y的列向量进行，返回值yi和矩阵y的列数相等，xi为插值点的横坐标，

22、yi时在xi指定位置计算出的插值结果。,45,4.4.1 插值,一维插值有四种方法，分别是： （1）邻近点插值(method=nearest) 将插值结果的值设置为最近数据点的值 （2）线性插值(method=linear) 在两个数据点之间连接直线，根据给定的插值点计算出它们在直线上的值，作为插值结果。缺省形式。 （3）三次样条插值(method=spline) 通过数据点拟合出三次样条曲线，根据给定的插值点计算出它们在曲线上的值，作为插值结果。 （4）立方插值(method=pchip/cubic) 通过三次多项式计算插值结果。,46,4.4.1 插值,由于在很多情况下，三次样条插值效果最

23、好，matlab 还专门提供了三次样条插值函数yi=spline(x,y,xi)，其中输入、输出参数含义同上。,47,4.4.1 插值,例4-13 一维插值函数插值方法的比较 clear x=0:2*pi;y=cos(x);xi=0:0.1:2*pi; %将插值方法定义成 单元数组 method=nearest,linear,spline,cubic lable=(a) method=nearest,(b) method=linear, (c)method=spline, (b) method=cubic; for i=1:4 yi=interp1(x,y,xi,methodi); %在一个图

24、形窗口绘制4幅图形 subplot(2,2,i),plot(x,y,ro,xi,yi,b),xlabel(lablei) end,48,4.4.1 插值,2、二维插值 二维插值与一维插值的基本思想相同，它是对两个自变量的函数z=f(x,y)进行插值。 matlab中二维插值函数为interp2，该函数的调用格式为： zi=interp1(x,y,z,xi,yi,method) 输入参数x，y为两个向量，z是矩阵，是由x和y确定的点的值z(i,:)=f(x,y(i) 和z(:,j)=f(x(j),y)，method为插值方法选项。,49,4.4.1 插值,二维插值有四种插值方法： （1）邻近点插

25、值(method=nearest) （2）双线性插值(method=linear) 该方法是interp2的缺省插值形式。 （3）三次样条插值(method=spline) （4）二重立方插值(method=cubic),50,4.4.1 插值,例4-14 二维插值方法比较 %二维插值方法比较 clear x,y,z=peaks(6); %生成双峰函数值 surf(x,y,z) % 画出表面图 xi,yi=meshgrid(-3:0.2:3,-3:0.2:3); % 生成供插值的数据 z1=interp2(x,y,z,xi,yi,nearest); z2=interp2(x,y,z,xi,yi,linear); z3=interp2(x,y,z,xi,yi,spline); z4=interp2(x,y,z,xi,yi,cubic); figure,surf(xi,yi,z1),title(nearest) % 绘制邻近