1、第六章 MATLAB的数值计算, matlab 具有出色的数值计算能力，占据世界上数值计算软件的主导地位,一、代数方程组求解,matlab中有两种除运算左除和右除。 对于方程ax=b，a 为amn矩阵，有三种情 况： 当n=m时，此方程成为“恰定”方程 当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=abx = x = 2.00 2.
2、003.00 3.00,=,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=23x1+4x2=3a=1 2;2 3;3 4;b=1;2;3;解1 x=ab 解2 x=inv(aa) a b x = x =1.00 1.000 0.00,=,a x = b,3.欠定方程组的解,当方程数少于未知量个数时,即不定 情况,有无穷多个解存在。 matlab可求出两个解
3、： 一个是用除法求的解，是具有最多零元素的解 一个是具有最小长度或范数的解，这个解是基于伪逆pinv求得的。,x1+2x2+3x3=12x1+3x2+4x3=2a=1 2 3;2 3 4;b=1;2;x=ab x=pinv(a)b x = x =1.00 0.830 0.330 -0.17,a x = b,借助matlab提供的函数，处理多项式是一件非常简单的事情，很容易对多项式进行积分、微分以及求根的操作。,一元多项式在代数中占有非常重要的地位。在实际应用中如对实验数据的插值、微商和曲线拟合等，都要大量用到多项式；在矩阵分析时，也要用到一元多项式的概念。多项式函数是形式最简单的函数，也是最容
4、易计算的函数，从理论上讲，它可以表示绝大多数复杂函数。在许多计算机的计算和编程中，很多函数值如sin(x),cos(x)等的计算都是先将函数进行Tailor展开为多项式进行逼近计算的，并且都能达到很高的精度。,二、 多项式运算及其求根,鉴于MATLAB无零下标，故把多项式的一般形式表达为：,在MATLAB中可以用长度n+1的行向量表示为：,1. 多项式求根,命令格式：x=roots(A)。这里A为多项式的系数A(1),A(2),A(N),A(N+1)；解得的根赋值给数组X，即X(1),X(2), ,X(N)。 【例6】试用ROOTS函数求多项式x4+8x3-10的根 这是一个4次多项式，它的五
5、个系数依次为：1,8,0,0,-10。下面先产生多项式系数的向量A，然后求根： A=1 8 0 0 -10 A =1 8 0 0 -10 x=roots(A),x =-8.0194 -0.5075 + 0.9736i-0.5075 - 0.9736i1.0344,2. 多项式的建立,若已知多项式的全部根，则可以用POLY函数建立起该多项式；也可以用POLY函数求矩阵的特征多项式。POLY函数是一个MATLAB程序，调用它的命令格式是： A=poly(x) 若x为具有N个元素的向量，则poly(x)建立以x为其根的多项式，且将该多项式的系数赋值给向量A。在此种情况下，POLY与ROOTS互为逆函
6、数；若x为NN的矩阵x，则poly(x)返回一个向量赋值给A，该向量的元素为矩阵x的特征多项式之系数：A(1),A(2),A(N),A(N+1)。,【例7】试用POLY函数对例7.8所求得的根，建立相应的多项式。 x=-8.0194 -0.5075 + 0.9736i -0.5075 - 0.9736i 1.0344; z=poly(x) z =1.0000 8.0000 0.0000 0.0000 -9.9996,3. 求多项式的值,POLYVAL函数用来求代数多项式的值，调用的命令格式为： Y=polyval(A,x)本命令将POLYVAL函数返回的多项式的值赋值给Y。若x为一数值，则Y也
7、为一数值；若x为向量或矩阵，则对向量或矩阵中的每个元素求其多项式的值。,【例8】以4次多项式、分别取x=1.2和下面的矩阵的23个元素为自变量计算该多项式的值。 A=1 8 0 0 -10; % 例7.8的4次多项式系数 x=1.2; % 取自变量为一数值 y1=polyval(A,x) y1 = -5.8976 x=-1 1.2 -1.4;2 -1.8 1.6 % 给出一个矩阵x y1=polyval(A,x) y1 = -17.0000 5.8976 -28.110470.0000 -46.1584 29.3216,4. 多项式的四则运算,(1)多项式加、减 对于次数相同的若干个多项式，可
8、直接对多项式系数向量进行加、减的运算。如果多项式的次数不同，则应该把低次的多项式系数不足的高次项用零补足，使同式中的各多项式具有相同的次数。,(2)多项式乘法 若A、B是由多项式系数组成的向量，则CONV函数将返回这两个多项式的乘积。调用它的命令格式为： C=conv(A,B) 命令的结果C为一个向量，由它构成一个多项式。,求4次多项式与多项式2x2-x+3的乘积。 A=1 8 0 0 -10; B=2 -1 3 B =2 -1 3 C=conv(A,B) C =2 15 -5 24 -20 10 -30 本例的运行结果是求得一个6次多项式 2x6+15x5-5x4+24x3-20x2+10x
9、-30,(3)多项式除法 当A、B是由多项式系数组成的向量时，DECONV函数用来对两个多项式作除法运算。调用的命令格式为： Q,r=deconv(A,B) 本命令的结果：多项式A除以多项式B获商多项式赋予Q（也为多项式系数向量）；获余项多项式赋予r（其系数向量的长度与被除多项式相同，通常高次项的系数为0）。 DECONV是CONV的逆函数，即有A=conv(B,Q)+r。,【例10】试用4次多项式与多项式2x2-x+3相除。 A=1 8 0 0 -10; B=2 -1 3; P,r=deconv(A,B) P = 0.5000 4.2500 1.3750 r = 0 0 0 -11.3750
10、 -14.1250 商多项式P为0.5x2+4.25x+1.375， 余项多项式r为 -11.375x-14.125。,polyder( ) 多项式微分,k = polyder(p) 求多项式p的微分多项式,k = polyder(a,b) 求a x b的微分,q,d = polyder(b,a)：求有理分式b/a的微分q：分子多项式d：分母多项式,(4).多项式微分,例：a=1 2 3 4 5; poly2str(a,x) ans = x4 + 2 x3 + 3 x2 + 4 x + 5 b=polyder(a) b = 4 6 6 4 poly2str(b,x) ans =4 x3 + 6
11、 x2 + 6 x + 4,10、多项式的积分,g = 1,6,20,48,69,72,44h = polyder(g)q = polyint(h,44),q = polyint(p,k) 对多项式p积分，k为积分常数 q = polyint(p) 对多项式p积分，积分常数为0,也可以用p./(length(p):-1:1),k来完成积分其中k为常数项例：a=1 2 3 4 5; b=polyder(a) b = 4 6 6 4 c= b./(length(b):-1:1),0 c = 1 2 3 4 0,三、 拟合与插值,1. 多项式拟合：用一个多项式来逼近一组给定的数据。拟合的准则是最小二
12、乘法 x0=0:0.1:1; y0=-.447 1.978 3.11 5.25 5.02 4.66 4.01 4.58 3.45 5.35 9.22; p=polyfit(x0,y0,3) p = 56.6915 -87.1174 40.0070 -0.9043 xx=0:0.01:1;yy=polyval(p,xx); plot(xx,yy,-b,x0,y0,or),polyfit函数来求得最小二乘拟合多项式的系数， polyval函数按所得的多项式计算所给出的点上的函数近似值。polyfit函数的调用格式为： P,S=polyfit(X,Y,m)函数根据采样点X和采样点函数值Y，产生一个m
13、次多项式P及其在采样点的误差向量S。其中X,Y是两个等长的向量，P是一个长度为m+1的向量，P的元素为多项式系数polyval(P,x)按多项式的系数计算x点多项式的值.,已知数据表t,y， t=0,1.7,2.9,3.5,4.8,6.9,7.3,8.6,9.4,10.2 y=5.6,7.8,14,18,19.6,18.6,17.6,15.4,12.5,14 试求2次拟合多项式p(t)，然后求ti=1,1.5,2,2.5,9.5,10各点的函数近似值。 解：源程序如下 T=0,1.7,2.9,3.5,4.8,6.9,7.3,8.6,9.4,10.2 Y=5.6,7.8,14,18,19.6,1
14、8.6,17.6,15.4,12.5,14 P=polyfit(T,Y,2) ti=1:0.5:10 yi=polyval(P,ti) plot(T,Y,r,ti,yi,b),多项式拟合时的注意事项,在进行曲线拟合时对多项式阶次的选择是任意的。,虽然高阶的多项式可以更准确地拟合数据（仅指该曲线与给定的数据之间的均方误差最小），但在进行曲线拟合时，并不需要采用太高阶的多项式，这主要基于以下原因：,1、越是高阶的多项式其数值特性越差，计算起来也越耗时；,2、随着多项式阶次的升高，拟合的曲线变的越来越不平滑，通常会出现用户不愿意看到的局部波形；,3、由于数据本身的近似性，因此在进行数据拟合时没有必要
15、仅仅考虑使拟合的曲线无限接近数据点，而要在曲线的阶次合均方误差之间综合考虑，因为越是高阶的多项式在物理实现时越困难。,从数学原理上，n+1个数据点可以惟一定义一个n次曲线（或n阶多项式）。,2.插值 插值是根据数据点的规律，找到一个多项式可以连接两个点，插值得出两个数据点之间的数值。 当不能很快地求出所需中间点的函数时，插值是一个非常有价值的工具。 Matlab提供了一维、二维、 三次样条等许多插值选择,y1=interp1(x,y,x1,方法); 例： x = 0:10; y = sin(x); xx = 0:.25:10; yy = spline(x,y,xx); yy=interp1(x
16、,y,xx,spline); plot(x,y,o,xx,yy),Y1=interp1(X,Y,X1,method)函数根据X,Y的值，计算函数在X1处的值。 X,Y是两个等长的已知向量，分别描述采样点和样本值， X1是一个向量或标量，描述欲插值的点， Y1是一个与X1等长的插值结果。 method是插值方法，允许的取值有 linear、nearest、cubic、spline。注意：X1的取值范围不能超出X的给定范围，否则，会给出“NaN”错误。,用不同的插值方法计算在/2点的值。 x=0,0.5,0.8,1.5,1.9,2.3 y=4,5,8,9,6,3 y1=spline(x,y,pi/
17、2) y1 =8.5931e+000 y2=interp1(x,y,pi/2,linear) y2 =8.4690e+000,某观测站测得某日6:00时至18:00时之间每隔2小时的室内外温度()，用3次样条插值分别求得该日室内外6:30至17:30时之间每隔2小时各点的近似温度()。解：设时间变量h为一行向量，温度变量t为一个两列矩阵，其中第一列存放室内温度，第二列储存室外温度。命令如下： h =6:2:18; t=18,20,22,25,30,28,24;15,19,24,28,34,32,30; XI =6.5:2:17.5 YI=interp1(h,t,XI,spline) %用3次样
18、条插值计算,二维数据插值在MATLAB中，提供了解决二维插值问题的函数interp2，其调用格式为： Z1=interp2(X,Y,Z,X1,Y1,method)X,Y是两个向量，分别描述两个参数的采样点， Z是与参数采样点对应的函数值， X1,Y1是两个向量或标量，描述欲插值的点。 Z1是根据相应的插值方法得到的插值结果。 method的取值 linear、nearest、cubic、spline X,Y,Z也可以是矩阵形式。 同样，X1,Y1的取值范围不能超出X,Y的给定范围，否则，会给出“NaN”错误。,某实验对一根长10米的钢轨进行热源的温度传播测试。用x表示测量点0:2.5:10(米
19、)，用h表示测量时间0:30:60(秒)，用T表示测试所得各点的温度()。试用线性插值求出在一分钟内每隔20秒、钢轨每隔1米处的温度TI。x=0:2.5:10; h=0:30:60; T=95,14,0,0,0;88,48,32,12,6;67,64,54,48,41; xi=0:10; hi=0:20:60; TI=interp2(x,h,T,xi,hi),四、泛函命令,函数句柄方便地实现函数间互相调用 兼容函数加载的所有方式 拓宽子函数包括局部函数的使用范围 提高函数调用的可靠性 减少程序设计中的冗余 提高重复执行的效率 数组、结构数组、细胞型数组结合定义数据,函数句柄的创建和显示定义函数
20、句柄的定义只需在提示符后添加相应函数的函数名f_h=plot 函数句柄的内容通过functions显示functions(f_h) ans = function: plottype: overloadedfile: MATLAB built-in functionmethods: 1x1 struct,matlab的函数类型simple 未加载的内部函数、m文件overloaded 加载的内部函数、m文件subfunction 子函数 constructor matlab类的创建函数private 局部函数,函数句柄与函数名字符串转换func2str 函数名字符串转换函数句柄funh=str2
21、func(det)funh = detstr2func 函数句柄转换函数名字符串func2str(funh)ans =det,(1)内联函数 inline(string,arg1,arg2) 例：f=inline(sin(x),x);f(1:10)(2)函数句柄 h=fun 或 h=str2func(fun) 例：f=aaa;f(0:pi/100;2*pi),1、寻优函数求单变量函数最小值x,y=fminbnd(h_fun,x1,x2)x,y=fminbnd(fun,x1,x2) 其中：y是计算出的函数最小值x是函数最小值的自变量值h_fun是函数句柄fun为函数名字符串x1、x2是自变量取值
22、范围,例：x = fminbnd(fun_name,2,5)x = fminbnd(sin(x)+3,2,5)x = fminbnd(x) sin(x)+3,2,5)f=inline(sin(x)+3,x); x = fminbnd(f,2,5),2、寻优函数求多变量函数无约束最小值（单纯形法）x,y=fminsearch(h_fun,x0)x,y=fminsearch(fun,x0) 其中：y是计算出的函数最小值向量x是函数最小值的自变量值向量x0是最小值点的初始猜测值,求: min f(x)=100(x2-x12)2+ (1-x1)2banana = (x)100*(x(2)-x(1)2)
23、2+(1-x(1)2; x0=-1.2, 1; x,fval = fminsearch(banana, x0)x = 1.0000 1.0000 fval = 8.1777e-010,47,3、求过零点： x=fzero(h_fun,x0) x=fzero(fun,x1,x2) 求函数在x0点附近或x1、x2之间的过零点（或等于某个常数点）的自变量值,4、数值积分： x=quad(h_fun, x1,x2) x=quad8(fun,x1,x2)x1、x2为积分上下限,5、执行函数： y1,y2,=feval(h_fun, x1,x2,) y1,y2,=feval(fun, x1,x2,)直接执
24、行函数,习题,1已知某班的5名学生的三门课成绩列表如下：学生序号 1 2 3 4 5 高等数学 78 89 64 73 68 外语 83 77 80 78 70 MATLAB语言 82 91 78 82 68 试写出有关命令，先分别找出三门课的最高分及其学生序号；然后找出三门课总分的最高分及其学生序号。2针对上小题的成绩表，求出其三门课总分存入数组ZF，再利用SORT命令对之按降序排序，同时把相应的学生序号存入数组XH。,习题,3今有多项式P1(x)=x4-2x+1,P2(x)=x2+4x-0.5，要求先求得P(x)=P1(x)+P2(x)，然后计算xi=0.2*i各点上的P(xi)（i=0,1,2,5）值。4. 练习多项式拟合与插值,
