第四章 数值计算 全文 文档投稿网

第四章

数值运算基础•

数值微积分•

线性方程组求解•

多多项项式式•

线性代数•

数数据据分分析析•

插插值值•

数字信号处理•近似数值极限及导数差分和梯度:dx=diff(X)//dx=X(2:n)-X(1:n-1)FX=grandient(F)//FX(1)=F(2)-F(1)FX(end)=F(end)-F(end-1),FX(2:end-1)=(F(3:end)-F(1:end-2))/2[FX,FY]=grandient(F)数值微积分•例1:已知x=sin(t),求该函数在区间[0,2Π]的近似导函数解1:利用导数定义（注意:自变量的增量不能太小）解2:采用diff和grandient方法（注意:diff方法缺少最后一个点）近似数值极限及导数•数值求和与近似数值积分Sx=sum(X)Scs=cumsum(X)//列方向累计求和St=trapz(x,y)//采用梯形法沿列方向计算y对x的积分Sct=cumtrapz(x,y)//采用梯形法沿列方向计算y对x的累计积分，即:Sct(k)=数值积分•例3:求积分，其中y(t)=0.2+sint解:分别利用sum和trapz计算•S1=quad(fun,a,b,tol)•S1=quad1(fun,a,b,tol)•S2=dblquad(fun,xmin,xmax,ymin,ymax,tol)•S3=triplequad(fun,xmin,xmax,ymin,yax,zmin,zmax,tol)•例:fx="exp(-x.^2)";ic=quad(fx,0,1,1e-8)精度可控的数值求积矩阵的秩矩阵的迹矩阵的行列式将矩阵A初等变换成阶梯阵RAX=0span

(Z)=span

(A)A的特征值与特征向量AV=VD线性方程组求解rank(A)trace(A)det(A)[R,ci]=rref(A)X=null(A)Z=orth(A)[V,D]=eig

(A)•

矩阵计算•1直接解法对于线性方程组Ax=b，可以利用左除运算符“\”求解:x=A\b例1用直接解法求解下列线性方程组。命令如下:A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-b=[13,-9,6,0]";x=A\b线性方程组求解2利用矩阵的分解求解线性方程组矩阵分解是指根据一定的原理用某种算法将一个矩阵分解成若干个矩阵的乘积。常见的矩阵分解有LU分解、QR分解、Cholesky分解,以及Schur分解、Hessenberg分解、奇异分解等。矩阵的LU分解就是将一个矩阵表示为一个交换下三角矩阵和一个上三角矩阵的乘积形式。线性代数中已经证明，只要方阵A是非奇异的，LU分解总是可以进行的。MATLAB提供的lu函数用于对矩阵进行LU分解，其调用格式为:[L,U]=lu(X):产生一个上三角阵U和一个变换形式的下三角阵L(行交换)，使之满足X=LU。注意，这里的矩阵X必须是方阵。LU分解[L,U,P]=lu(X):产生一个上三角阵U和一个下三角阵L以及一个置换矩阵P，使之满足PX=LU。当然矩阵X同样必须是方阵。实现LU分解后，线性方程组Ax=b的解x=U\(L\b)或x=U\(L\Pb)这样可以大大提高运算速度。,例2用LU分解求解例7-1中的线性方程组。命令如下:A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-b=[13,-9,6,0]";[L,U]=lu(A);x=U\(L\b)或采用LU分解的第2种格式，命令如下:[L,U,P]=lu(A);x=U\(L\P*b)对矩阵X进行QR分解，就是把X分解为一个正交矩阵Q和一个上三角矩阵R的乘积形式。QR分解只能对方阵进行。MATLAB的函数qr可用于对矩阵进行QR分解，其调用格式为:[Q,R]=qr(X):产生一个正交矩阵Q和一个上三角矩阵R，使之满足X=QR。[Q,R,E]=qr(X):产生一个一个正交矩阵Q、一个上三角矩阵R以及一个置换矩阵E，使之满足XE=QR。实现QR分解后，线性方程组Ax=b的解x=R\(Q\b)或x=E(R\(Q\b))。QR分解例3用QR分解求解例7-1中的线性方程组。命令如下:A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-b=[13,-9,6,0]";[Q,R]=qr(A);x=R\(Q\b)或采用QR分解的第2种格式，命令如下:[Q,R,E]=qr(A);x=E*(R\(Q\b))如果矩阵X是对称正定的，则Cholesky分解将矩阵X分解成一个下三角矩阵和上三角矩阵的乘积。设上三角矩阵为R，则下三角矩阵为其转置，即X=R"R。MATLAB函数chol(X)用于对矩阵X进行Cholesky分解，其调用格式为:R=chol(X):产生一个上三角阵R，使R"R=X。若X为非对称正定，则输出一个出错信息。[R,p]=chol(X):这个命令格式将不输出出错信息。当X为对称正定的，则p=0，R与上述格式得到的结果相同；否则p为一个正整数。如果X为满秩矩阵，则R为一个阶数为q=p-1的上三角阵，且满足R"R=X(1:q,1:q)。实现Cholesky分解后，线性方程组Ax=b变成R‘Rx=b，所以x=R\(R’\b)。Cholesky分解命令如下:A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-b=[13,-9,6,0]";R=chol(A)???Errorusing==>cholMatrixmustbepositivedefinite命令执行时，出现错误信息，说明A为非正定矩阵。用Cholesky分解求解例1中的线性方例4

程组。•

多多项项式式表表达达方方式式的的约约定定•

创创建建多多项项式式的的方方法法•

多多项项式式运运算算函函数数多项式•MATLAB约定降幂多项式P(x)=a0xn+a1xn-1+…+an-1x+an用以下系数矢量（系数行向量）表示:p=[a0,a1,…,an-1,an],即把多项式的各项系数依降幂次序排放在行向量的元素位置上。•注意:假如多项式中缺某幂次项，则应认为该幂次项的系数为零。多项式表达方式的约定•系数矢量的直接输入法–在命令窗直接输入多项式的系数矢量，然后利用转换函数poly2sym将多项式由系数矢量形式转换为符号形式。•特征多项式输入法–由矩阵的特征多项式取得，由函数poly实现。•由根矢量创建多项式–由给定的根矢量创建多项式，由函数poly实现。–注意:•若希望创建实系数多项式，根矢量的复数根必须共轭成对。•含复数根的根矢量所创建的多项式系数矢量的系数中，有可能带有很小（在截断误差数量级）的虚部，此时可采用取实部的命令（real）把虚部滤掉。创建多项式的方法指

令含

义pA=po

lyval(p

,S)按数组运算规则计算多项式值

。P为多项式

，

S

为矩阵PM=po

lyvalm(p

,S)按矩阵运算规则计算多项式值

。P为多项式

，

S

为矩阵r=roots(p)求多项式p的根p=

conv(p1

,p2)p是多项式p1和p2的乘积多项式[q

,r]=

deconv(p1

,p2)多项式p1除以p2商多项式为q

,而余多项式是r[r

,p

,k]=residue

(b

,a)部分分式展开：

b

,a分别是分子和分母多项式系

数向量

;r

,p

,k分别是留数

、极点

、直项p=po

ly(AR)求方阵AR的特征多项式p；

或求向量AR指定根

所对应的多项式p=po

lyfit

(x

,y

,n)求x

,y

向量给定的数据的n阶拟合多项式p多项式运算函数的调用格式•p1=conv(conv([1,0,2],[1,4]),[1,p2=[1011][q,r]=deconv(p1,p2)验证:qp1=conv(p2,q)qp2=qp1+rqp2==p11•A=[111213;141516;171819]pA=poly(A)ppA=poly2str(pA,’s’)比较:s=eig(A)s=45.3965-0.39650.0000r=roots(pA)r=45.3965-0.39650.0000例2:求3阶方阵A的特征多项式。r=

[1,2]p=poly(r)poly2str

(p,

"x

")r

=1

2p=1

-3

2an

s=x

^2

-

3

x

+

2例3:由给定根向量求多项式系数向量。•例1:对于给定数据对（x0,y0），求拟合三阶多项式，并图示拟合情况。x0=0:0.1:1;y0=[-.447,1.978,3.11,5.25,5.02,4.66,4.01,4.58,3.45,n=3p=polyfit(x0,y0,n)xx=x0;yy=polyval(p,xx);plot(xx,yy,"-b",x0,y0,".r","MarkerSize",20xlabel("x")图解拟合和最小二乘法x0=0:0.1:1;y0=[-.447,1.978,3.11,5.25,5.02,4.66,n=3p=polyfit(x0,y0,n)xx=x0;yy=polyval(p,xx);plot(xx,yy,"-b",x0,y0,".r","MarkerSixlabel("x")xx1=0:0.01:1;p2=interp1(x0,y0,xx1,"cubic")holdonplot(xx1,p2,":k","LineWidth"