第四章 数值计算

第四章数值计算第一页，共五十页，2022年，8月28日1数值计算内容:数值计算包括:多项式运算,线性方程组求解,矩阵特征值问题的解,卷积,数据分析,泛函指令的使用,信号处理和系统分析.4.1多项式运算多项式运算是数学中最基本的运算之一.多项式一般可以表示为:f(x)=a0xn+a1xn-1+a2xn-2+…+an-1x+an对于这种表示形式，很容易用一个行向量表示，即：T=[a0,a1,a2,…an-1,an]在MATLAB中,多项式正是用这样一个行向量表示的,系数是按降序排列的．幂第二页，共五十页，2022年，8月28日24.1.1多项式构造多项式可以直接用向量表示,因此,构造多项式最简单的方法是直接输入向量.例4.1.1-1直接输入向量构造f(x)=2x5+5x4+4x2+x+4T=[2,5,0,4,1,4];fx=poly2sym(T)函数poly2sym是符号工具箱中的函数,在用此种方式构造多项式时,必须把多项式各项的系数写完整,而不管此项的系数是否为0fx=2*x^5+5*x^4+4*x^2+x+4第三页，共五十页，2022年，8月28日3r=[1,2,3,4];T1=poly(r);y=poly2sym(T1)y_class=class(y)例4.1.1-2用多项式的根构造多项式,根为r=[1,2,3,4]T1=1-1035-5024y=x^4-10*x^3+35*x^2-50*x+24y_class=sym第四页，共五十页，2022年，8月28日44.1.2多项式的运算多项式的运算主要包括多项式的四则运算,导数运算,估值运算,求根运算以及多项式的拟合等．1>多项式的四则运算多项式四则运算主要是多项式的加,减,乘,除.其中,多项式的加减运算要求两个相加,减的多项式的大小必须相等,此时加,减才有效.当两个相加,减的多项式阶次不同时必须通过补0使其相同.加/减---+/-乘/除---conv,deconvT=deconv(T1,T3)[T,r]=deconv(T1,T3)商多项式余式

T3为分母第五页，共五十页，2022年，8月28日5

T1=[2,5,0,4,1,4]; T2=[0,0,5,1,3,2]; T3=[5,1,3,2];

%除法运算中分母多项式第一个系数不能为0 T=T1+T2;

%

必须是同维的才能相加 T_add=poly2sym(T) T=T1-T2; T_sub=poly2sym(T) T=conv(T1,T2);

%乘法不要求同维 T_mul=poly2sym(T) [A_coe,A_r]=deconv(T1,T3); T_coe=poly2sym(A_coe) T_rem=poly2sym(A_r)

例4.1.1-3多项式的加减乘除运算f1(x)=2x5+5x4+4x2+x+4,f2(x)=5x3+x2+3x+2第六页，共五十页，2022年，8月28日6例4.1.1-4多项式求值,求上式f1(x)在x=0.5处的函数值

T1=[2,5,0,4,1,4]; x=0.5; y=polyval(T1,x)y=5.8750第七页，共五十页，2022年，8月28日73>多项式的导数运算---polyder例4.1.1-6求多项式f1(x)=2x5+5x4+4x2+x+4的导数T1=[2,5,0,4,1,4];h=polyder(T1);poly2sym(h)2>多项式求根---roots例4.1.1-5求多项式f1(x)=2x5+5x4+4x2+x+4的根T1=[2,5,0,4,1,4];root=roots(T1);root=-2.77090.5611+0.7840i0.5611-0.7840i-0.4257+0.7716i-0.4257-0.7716i10*x^4+20*x^3+8*x+1第八页，共五十页，2022年，8月28日84>拟合和插值---polyfit,interp1例4.1.1-7对下列数据对(x0,y0)求三次拟合多项式并绘图.x0=0:0.1:1;y0=[-0.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,n);[p,s]=polyfit(x0,y0,n);x0,y0---给定数据对n----拟合出的多项式次数p---多项式向量s---偏差信息yi=interp1(x0,y0,xi,'cubic');xi,yi---得到的新的插值点cubic---插值方式第九页，共五十页，2022年，8月28日9x0=0:0.1:1;y0=[-0.447,1.978,3.11,5.25,5.02,4.66,4.01,4.58,3.45,5.35,9.22];n=3;%设定拟合次数为3[p,s]=polyfit(x0,y0,n);%得到拟合多项式向量和相关偏差信息xx=0:0.01:1;yy=polyval(p,xx);%按拟合曲线计算采样值n1=6;%设定拟合次数为6[p1,s1]=polyfit(x0,y0,n1);yy1=polyval(p1,xx);plot(xx,yy,'-b',xx,yy1,'-m',x0,y0,'.r','MarkerSize',20);y3=polyval(p,0.5);y6=polyval(p1,0.5);text(0.5,y3-0.3,'n=3');text(0.5,y6+0.2,'n=6');第十页，共五十页，2022年，8月28日10y=51.48x3-77.74x2y=348.21x6-1060.23x5+1297.68x4-758.90x3+181.00x2+1.00x+0.48第十一页，共五十页，2022年，8月28日11拟合只能在给定数据所限定的区间内使用,不要任意向外拓展．例4.1.1-8按上例所给数据研究插值,并观察插值和拟合的区别.x0=0:0.1:1;y0=[-0.447,1.978,3.11,5.25,5.02,4.66,4.01,4.58,3.45,5.35,9.22];xi=0:0.02:1;yi=interp1(x0,y0,xi,'cubic');subplot(1,2,1);plot(xi,yi,'-b',x0,y0,'.r','MarkerSize',20);xlabel('x'),ylabel('y');[p,s]=polyfit(x0,y0,3);xx=0:0.01:1;yy=polyval(p,xx);subplot(1,2,2);plot(xx,yy,'-r',x0,y0,'.r','MarkerSize',20);xlabel('x');ylabel('y');第十二页，共五十页，2022年，8月28日12插值拟合第十三页，共五十页，2022年，8月28日13插值方式:linear---线性插值nearest---最近插值cubic---三次多项式插值spline---样条插值pchip---分段三次Hermite插值v5cubic---MATLAB5.0版的三次多项式插值.拟合和插值的区别:曲线拟合是研究如何寻找“平滑”曲线最好地表现带噪声的“测量数据”,但并不要求拟合曲线通过这些“测试数据”点.插值是在认定所给数据完全正确的情况下,研究如何平滑地估算出“基准数据”之间其他点的函数值,因此插值一定通过“基准数据”第十四页，共五十页，2022年，8月28日14通俗地讲,拟合就是由已知点得到一条曲线,使该曲线最能反映点所代表的规律.比如做欧姆定理的实验的时候,由于实验中存在误差,最后拟合得到的曲线是一条直线,而且肯定只有部分点落在拟合的直线上,但此时该直线和测试点的方差最小.由拟合直线的斜率就可以知道电阻的阻值.拟合是探测事物变化规律的办法.插值就是根据函数上某些已知点(或实验数据),按一定规律(插值方法)寻求未知的点,比如已知一个常用对数y=log(x)表,是按照x=0.1:0.1:10制表的,如果按已知数据求y=log(2.897)就可以用插值得到.表制得越密,插值越准确．第十五页，共五十页，2022年，8月28日15例4.1.1-9已知通过实验得到电压和电流的数值如下：i0=0:0.5:10;v0=[0,98,202,306,395,506,592,708,789,890,1009,1108,1186,1310,1391,1510,1601,1716,1782,1920,2014];按一次,二次,三次多项式对数值关系进行拟合．说明：为保证较好的拟合效果，多项式阶数要取得适当，过低，残差较大，过高，拟合模型将包含噪声影响，通常要保证拟合阶数小于数据对的数目。第十六页，共五十页，2022年，8月28日16i0=0:0.5:10;v0=[0,98,202,306,395,506,592,708,789,890,1009,1108,...1186,1310,1391,1510,1612,1716,1782,1920,2014];n1=1;[p1,s]=polyfit(i0,v0,n1);n2=2;[p2,s]=polyfit(i0,v0,n2);n3=3;[p3,s]=polyfit(i0,v0,n3);ii=0:0.1:10;v1=polyval(p1,ii);v2=polyval(p2,ii);v3=polyval(p3,ii);subplot(1,3,1),plot(ii,v1,'-b',i0,v0,'.r','MarkerSize',8),xlabel('一次拟合');subplot(1,3,2),plot(ii,v2,'-b',i0,v0,'.r','MarkerSize',8),xlabel('二次拟合');subplot(1,3,3),plot(ii,v3,'-b',i0,v0,'.r','MarkerSize',8),xlabel('三次拟合');p1=[201.00,-2.89]P2=[0.31,197.90,2.02]P3=[0.03,-0.10,199.49,0.85]第十七页，共五十页，2022年，8月28日17p1=[201.00,-2.89]P2=[0.31,197.90,2.02]P3=[0.03,-0.10,199.49,0.85]第十八页，共五十页，2022年，8月28日18对于含n个未知数的n个方程构成的方程组Ax=b,在线性代数教科书中,最常介绍的解法有:Cramer法;逆阵法,即x=A-1b;高斯消元法;LU法．

对于维数不高,条件数不大的矩阵,以上四种所得结果不会有明显的差别.但前三种解法的更多的意义在理论上,而不在实际的数值计算上.在MATLAB中,方程采用LU法求解,并且出于算法稳定性的考虑,行列式和逆的计算也都是在LU分解的基础上进行的．MATLAB中的四条指令:[L,U,P]=lu(A)---矩阵的LU分解,使LU=PA.多种格式.det(A)---求矩阵A的行列式,它由detA=±detU＝±∏uii算得inv(A)---求矩阵A的逆,由A-1=U-1L-1P算得x=A\b---除法指令求方程的解(对恰定方程,采用LU法执行)4.2.1方程组求解4.2线性方程组的解第十九页，共五十页，2022年，8月28日19对于方程组Ax=b,采用x=A\b计算，如果方程组为yC=d,要使用右除，即指令为y=d/CAx=bx'A'=b'yC=dx=A\bx'=b'/A'y=d/C例4.2.1-1解下列方程组2x1+2x2+3x3=34x1+7x2+7x3=1-2x1+4x2+5x3=-7A=[2,2,3;4,7,7;-2,4,5],b=[3;1;-7]x=A\b,y=b'/A'求解方程时,尽量不要使用指令inv(A)*b进行,它不仅计算速度没有除法快,而且计算精度也没除法高.此外,除法还适用于“超定”和“欠定”方程．x=9.00-36.0031.00y=9.00-36.0031.00第二十页，共五十页，2022年，8月28日20方程组的解的三种情况:对于方程Ax=b,A为Am×n矩阵,有三种情况:当m=n时,此方程成为"恰定"方程,求解精确解

当m>n时,此方程成为“超定”方程,寻求最小二乘解

(直线拟合)当m<n时,此方程成为"欠定"方程,寻求基本解

matlab定义的除运算可以很方便地解上述三种方程1)恰定方程组的解方程组Ax=b(A非奇异)，解为x=A\bA=[2,-1,-1;3,4,-2;3,-2,4];b=[4;11;11];det(A),rank(A),rank([A,b])x=A\b例4.2.1-2求下列方程组的解2x1-x2-x3=43x1+4x2-2x3=113x1-2x2+4x3=11x=3.001.001.00第二十一页，共五十页，2022年，8月28日212)超定方程的解方程Ax=b,m>n时此时无解,但实用中可以求最小二乘解即:方程解(A'A)x=A'bx=(A'A)-1A'b---求逆法x=A\b---matlab用最小二乘法找一个近似解.x1+2x2=12x1+3x2=23x1+4x2=3例4.2.1-3求下列超定方程的解<1>求逆法：x=inv(A'*A)*A'*b<2>使用A\b的方式A=[1,2;2,3;3,4];b=[1;2;3];x1=inv(A'*A)*A'*b;x2=A\bx1=1.000.00x2=1.000第二十二页，共五十页，2022年，8月28日223)欠定方程当方程数少于未知量个数时,即不定情况,有无穷多个解存在.matlab可求出两个解：a.用除法求的解x是具有最多零元素的解b.具有最小长度或范数的解,这个解是基于伪逆pinv求得的.例4.2.1-4求下列欠定方程组的解x1+2x2+3x3=12x1+3x2+4x3=2A=[1,2,3;2,3,4];b=[1;2];x1=A\b;x2=pinv(A)*bx1=100x2=0.83330.3333-0.1667第二十三页，共五十页，2022年，8月28日234.2.2奇异值分解和矩阵结构相应的MATLAB指令如下:[U,S,V]=svd(A)---矩阵的奇异值分解三对组阵，使A=USVTnorm(A,flag)---计算矩阵A的范数，范数类型由flag指定.cond(A)---计算矩阵A的条件数.pinv(A)---求矩阵的广义逆flag---1,2,inf,‘fro’第二十四页，共五十页，2022年，8月28日244.2.3线性二乘问题的解对于超定方程,求其最小二乘解有三种方法：①正则方程法得解x=(A'

A)-1A'

b,数值精度差②广义逆法得解x=A+b,所得解可靠.③用矩阵除法得解x=A\b.可靠性稍差，但速度快.例4.2.3-1对于超定方程y=Ax,进行三种解法比较.(1)生成矩阵A及y,并用三种方法求解．A=gallery(5),A(:,1)=[];y=[1.7,7.5,6.3,0.83,-0.082]';x=inv(A'*A)*A'*y,xx=pinv(A)*y,xxx=A\yx=3.48115.15950.9534-0.0466xx=3.47595.19480.7121-0.1101xxx=3.46055.29870-0.2974由于Ａ阵缺一个列秩,除法给出的最小二乘解就只有三个非零元素.它是只有三个非零元素的所有最小二乘解中范数最小的．Gallery--------MATLAB设置的特殊矩阵.

第二十五页，共五十页，2022年，8月28日25(2)计算三个解的范数nx=norm(x),nxx=norm(xx),nxxx=norm(xxx)nx=6.2968nxx=6.2918nxxx=6.3356用广义逆所得的最小二乘范数最小(3)比较三种解法的方程误差e=norm(y-A*x),ee=norm(y-A*xx),eee=norm(y-A*xxx)e=0.6986ee=0.0474eee=0.0474误差的平方根最大，精度最差第二十六页，共五十页，2022年，8月28日264.2.4特征值分解和矩阵函数涉及矩阵特征值分解的常用的MATLAB指令如下:d=eig(A)---计算矩阵A的特征值,以向量形式存放[V,D]=eig(A)---计算矩阵的特征向量V和特征值对角阵D,使A*V=V*D[VR,DR]=cdf2rdf(VC,DC)---把复数对角形转化成实数块对角形.[VC,DC]=rsf2csf(VR,DR)---把实数块对角形转换成复数对角形[V,J]=jordan(A)---jordan分解,使A*V=V*J,J是对角阵．c=condeig(A)---计算各特征值的条件数．[V,D,C]=condeig(A)---相当于[V,D]=eig(A)和c=condeig(A)两条指令第二十七页，共五十页，2022年，8月28日27例4.2.4-1验证det(A)=ПλiA=[2,-1,-1;3,4,-2;3,-2,4]det(A)d=eig(A)ans=60d=2.0000+2.4495i2.0000-2.4495i6.0000d(1)*d(2)*d(3)ans=60.0000行列式特征值%矩阵行列式等于其所有特征值的积第二十八页，共五十页，2022年，8月28日284.3数据分析函数(1)若输入宗量x是向量,那么不管是行向量还是列向量,运算是对整个向量进行的．(2)若输入宗量x是二维数组，那么指令运算是按列进行的．MATLAB在进行数据分析时的约定:第二十九页，共五十页，2022年，8月28日294.3.1随机数发生器和统计分析指令函数含义rand(m,n)产生mxn维的[0,1]区间均匀分布随机数组．randn(m,n)产生mxn维的均值为0标准差为1的正态分布随机数组min(x)对矩阵各列分别求最小值max(x)对矩阵各列分别求最大值median(x)对矩阵各列分别求中位数mean(x)对矩阵各列分别求均值std(x)对矩阵各列分别求标准差var(x)对矩阵各列分别求方差C=conv(x)给出矩阵各列的协方差阵P=corrcoef(x)给出矩阵各列间的相关系数第三十页，共五十页，2022年，8月28日30例4.3-1基本统计示例randn('state',0);A=randn(1000,4);AMAX=max(A),AMIN=min(A)AMED=median(A);AMEAN=mean(A);ASTD=std(A)已知A=[1,2,3;4,5,6;7,8,9],演示max,min,median,mean,std,var,conv,corrcoef各列最大,最小值：max(A),min(A)各行最大,最小值：max(A'),min(A')各列中间值：median(A)各列平均值：mean(A)各列方差,标准差：var(A),std(A)A=[1,2,3;4,5,6;7,8,9]789第三十一页，共五十页，2022年，8月28日314.3.2差分和累计指令函数含义del2(U,hx,hy)五点Laplaciandiff(X,m,n)沿第n维求m次差分(二维：n=1;列,n=2,行)[DZx,DZy]=gradient(z,dx,dy)对Z求x,y方向梯度prod(X,n)沿第n维求乘积sum(X,n)沿第n维求和trapz(x,Y,n)用梯形法沿第n维求函数关于自变量x的积分cumprod(X,n)沿第n维求累计乘积cumsum(X,n)沿第n维求累计和cumtrapz(x,Y,n)用梯形法沿第n维求函数Y关于x的累计积分[XS,KK]=sort(X,n)沿第n维对X元素按模增大排列第三十二页，共五十页，2022年，8月28日32diff,trapz,cumtrapz指令的演示a=[1,2,3,4,5];b=[a',a'+1,a'+2,a'+3]diff(a)=diff(a,1)=[1,1,1,1]diff(b)=[1,1,1,1;1,1,1,1;1,1,1,1;1,1,1,1]trapz(b)=[12162024]cumtrapz(b)

00001.502.503.504.504.006.008.0010.007.5010.5013.5016.5012.0016.0020.0024.0012342345345645675678第三十三页，共五十页，2022年，8月28日33例4.3.2-1求1+2+3+…+100以及50!

x=1:1:100;sum_x=sum(x);sum_x=5050a=1:1:50;prod_a=prod(a)prod_a=3.0414e+064cumsum,cumprod的用法上两个函数分别是求向量的累计和和累计乘积.如果a=1:1:n,则cumsum(a)=[1,1+2,1+2+3,...,1+2+3+…+n]cumprod(a)=[1,1*2,1*2*3,…,1*2*3*…*n]第三十四页，共五十页，2022年，8月28日34例4.3.2-1求f(x)=3x2在区间[0,2]的积分dt=0.001;t=(0:dt:2)';y=3*t.^2;s1=dt*sum(y)s2=dt*trapz(y)s=dt*cumsum(y);s3=s(end)s=dt*cumtrapz(y);s4=s(end)matlab还有更精良的积分指令:quad,quda8s1=8.0060s2=8.000001s3=8.0060s4=8.000001第三十五页，共五十页，2022年，8月28日354.4MATLAB泛函指令在MATLAB中,凡以函数为输入宗量的指令,都被统称为泛函指令.最常见的有:z=fzero(fun,x0)---求一元函数零点指令的最简单格式．x=fsolve(fun,x0)---解非线性方程组的最简单格式x=fminbnd(fun,x1,x2)---求函数在区间(x1,x2)中极小值的指令最简格式．x=fminsearch(fun,x0)---单纯形法求多元函数极值点指令的最简格式.x=fminunc(fun,x0)---拟牛顿法求多元函数极值点指令的最简格式a=lsqnonlin(fun,a0)---解非线性最小二乘问题指令的最简格式.第三十六页，共五十页，2022年，8月28日36q=quad(fun,a,b)---采用递推自适应simpson法计算积分．q=quadl(fun,a,b)---采用递推自适应Lobatto法求数值积分．SS=dblquad(fun,inmin,inmax,outmin,outmax)---二重闭环数值积分[t,YY]=ode45(fun,tspan,Y0)---采用4,5阶Runge-Kutta方程解算ODE初值问题指令中被处理的函数fun,可以取:字符串表达式,内联函数,M函数文件的函数句柄第三十七页，共五十页，2022年，8月28日374.4.1求函数零点对于任意函数f(x)=0,零点情况复杂,没有一个通用解法.一般来说,零点的数值计算过程是:先猜测一个初始零点或该零点所在的区间,然后通过计算,使猜测值不断精确化,或使猜测区间不断收缩,直到达到预定的精度,终止计算.求函数的零点的步骤为:(1)利用matlab作图指令获取初步近似值具体做法:先确定一个零点可能存在的自变量区间,然后利用plot指令画出f(x)在该区间中的图形,观察f(x)与横轴的交点坐标,也可以用zoom对交点处局部放大再读数.借助ginput指令获得更精确的交点坐标值．第三十八页，共五十页，2022年，8月28日38(2)利用fzero指令求精确解z=fzero(fun,x0)---求一元函数零点指令的最简格式[z,f_z,exitflag]=fzero(fun,x0,options,p1,p2,…)---最完整格式

说明：●由于fzero是根据函数是否越横轴来决定零点的,因此本程序无法确定函数曲线仅触及横轴和不穿越的零点．●第二个输入量x0表示零点的初始猜测.他可以是标量或二元向量.当x0取标量时,该指令将在它两侧寻找一个与之最靠近的零点;当取二元向量[a,b]时,该指令将在区间[a,b]内寻找一个零点.●输入宗量options是优化迭代采用的参数.P1,P2是向函数fun传递的参数.●z是所求零点的自变量值.f_z是函数值,exitflag表明程序终止计算的条件.若exitflag>0,表明找到零点后退出．第三十九页，共五十页，2022年，8月28日39例4.4.1-1求函数f(x)=2x-10x的零点(1)绘制函数图形x=-10:0.1:10;y=2.^x-10*x;holdon;plot(x,y,'r','LineWidth',5);plot(x,zeros(size(x)),'k');%便于查看零点xlabel('x'),ylabel('y');title('曲线方程式为：y=2^x-10*x');(2)获取初步近似解zoomon;[tt,yy]=ginput(2);

zoomoff;第四十页，共五十页，2022年，8月28日40func='2^x-10*x';[x01,err1,exitflag1]=fzero(func,tt(1),[]);[x02,err2,exitflag2]=fzero(func,tt(2),[]);(3)利用fzero指令求精确解str=['方程的根为:

x1=',num2str(x01),'和x2=',num2str(x02)];disp(str);第四十一页，共五十页，2022年，8月28日414.4.2求函数极值点许多科学研究和工程计算问题都可以归结为一个极值问题,如能量最小，时间最短,最佳拟合,最短路径等.又因为f(x)的极小值问题等价于-f(x)的极大值问题,所以matlab只有处理极小值的指令．确切地说,这里讨论的只是“局域极值”问题.“全域最小”问题要复杂得多.至今没有一个系统性的方法可求解一般的“全域最小”问题.对于一元二元函数,作图观察对确定全域最小值有很好的应用价值.但更多元的函数,就很难使用作图法.matlab求函数极值的三条指令[x,fval,exitflag,output]=fminbnd(fun,x1,x2,opt,p1,p2,…)---求一元函数在区间(x1,x2)中极小值的指令的最完整格式．[x,fval,exitflag,output]=fminsearch(fun,x0,opt,p1,p2,…)---单纯形法求多元函数极值点指令的最完整格式．[x,fval,exitflag,output,grad,hessian]=fminunc(fun,x0,opt,p1,p2,…)---拟牛顿法求多元函数极值点指令的最完整格式．第四十二页，共五十页，2022年，8月28日42例4.4.2-1求例4.4.1-1中的函数的最小值(1)函数采用字符串表达式func='2^x-10*x';x1=0,x2=8;%通过观察,发现最小值在此区间[Xmin,Ymin]=fminbnd(func,x1,x2)Xmin=3.8507Ymin=-24.0800(2)函数采用M文件函数的句柄%M文件内容，保存为func.mfunctiony=f(x)y=2^x-10*x;%M文件内容x1=0;x2=8;Hf=@func;[Xmin,Ymin]=fminbnd(Hf,x1,x2)第四十三页，共五十页，2022年，8月28日434.4.3数值积分数值积分有闭型和开型算法之分．两者的区别在于：是否需要计算积分区间端点处的函数值．matlab中仅提供闭型数值积分指令:q=quad(fun,a,b,tol,trace,p1,p2,…)---采用递推自适应Simpson法计算积分q=qudal(fun,a,b,tol,trace,p1,p2,…)---采用递推自适应Lobatlo法求数值积分SS=dblquad(fun,inmin,inmax,outmin,outmax,tol,method)---二重积分指令

fun---被积函数,a,b---积分上下限tol---绝对误差控制量，缺省时积分的绝对精度为1x10-6trace---取非０时，将随积分的进程逐点画出被积函数第四十四