三讲matlab解方程和函数极值

《三讲matlab解方程和函数极值》由会员分享，可在线阅读，更多相关《三讲matlab解方程和函数极值（52页珍藏版）》请在金锄头文库上搜索。

1、第三讲解方程与函数极值 方程问题和极值问题始终 是数学问题中的核心问题 ！ 学习内容 l多项式运算 l线性方程组求解 l非线性方程数值求解 l常微分方程初值问题的数值解法 l函数极值 l线性插值 （一）多项式的表示方法 n对于多项式的表达式约定如下 n对于多项式 n对于上述多项式一般用以下行向量表示 （二）多项式的创建 1.系数向量直接输入法 由于MATLAB中多项式是以向量形式存储的,因此 最简单的多项式输入即向量输入。 例:输入多项式 p=1-5 6 -33； poly2sym(p)%此函数将多项式转换换为符号 多项式 （二）多项式的创建 2.特征多项式输入法 多项式创建的另一途径是从矩阵

2、求其特征多项 式获得，由函数poly实现 例如： a=1 2 3;2 3 4;3 4 5 p1=poly(a)%矩阵a对应的特征多项式 poly2sym(p1)%将多项式p1转换为符号多 项式 （二）多项式的创建 3.由根创立多项式 给定的根也可产生其相应的多项式,此功能还是 由函数poly实现 例：root=-5 -3+4i -3-4i; p=poly(root) poly2sym(p) 注：若要生成实系数多项式,则根中的复数比为对 称共轭复数 （三）多项式运算 1.求多项式的值 求多项式的值可以由两种形式 (1)输入变量值代入多项式计算时是以数组为单 元的,此时函数为polyval (2)

3、以矩阵为计算单元,进行矩阵式运算以求的 多项式的值,此时的函数为polyvalm 注：这两种计算在数值上由很大差别。 当进行矩阵运算时,变量矩阵须为方阵。 实例演示 p=1 11 55 125;b=1 1;1 1;c=5 poly2sym(p) polyval(p,b) poly(p,c) polyvalm(p,b) （三）多项式运算 2.求多项式的根 求多项式的根可以由两种方法 (1)直接调用MATLAB的函数roots求解多 项式的所有根 (2)通过建立多项式的伴随矩阵再求其特 征值的方法得到多项式的所有根 实例演示 n用两种方法解方程 的所有根 p=2 -5 6 -1 9; roots(

4、p) a=compan(p)%求多项式的特征矩阵 eig(a)%求特征矩阵a的特征根 （三）多项式运算 3.多项式的乘除法 (1)多项式的乘法由函数conv来实现,此函数同用 于向量的卷积 (2)多项式的除法由函数deconv来实现,此函数 与向量的解卷函数相同 实例演示 n计算两多项式的乘除法 np=2 -5 6 -1 9; npoly2sym(p); nd=3 -90 -18; nploy2sym(d) npd=conv(p,d)%多项式p与d相乘 npoly2sym(pd) npl=deconv(pd,d)%多项式pd除以d （三）多项式运算 4.多项式微分 多项式的微分可以用函数pol

5、yer进行 例：p=2 -5 6 -1 9；poly2sym(p); Dp=polyer(p) （三）多项式运算 5.多项式的拟合 多项式拟合是多项式运算的一个重要组成部分，在工 程及科研工作中得到了广泛的应用，其一方面可以 由矩阵的除法解超定方程来进行;另一方面在 MATLAB中还提供了专门的拟合函数polyfit,其调用格 式如下： (1)polyfit(X,Y,n)%其中X,Y为拟合数据,n为你和多项式 的阶数 (2)p,s=polyfit(X,Y,n)%其中p为拟合多项式的系数向 量，s为拟合多项式系数向量的结构信息 实例演示 例:用五阶多项式对0,pi/2上的正弦函数进行最小二乘拟合

6、 x=0:pi/20:pi/2;y=sin(x); a=ployfit(x,y,5);%用五阶多项式拟合y=sin(x),a为你和多项 式的系数 x1=0:pi/30:pi/2;y1=sin(x1); y2=a(1)*x1.5+ a(2)*x1.4+ a(3)*x1.3+ a(4)*x1.2+ a(5)*x1. +a (6) plot(x1,y1,b-,x2,y2,r*) legend(原曲线，拟合曲线) axis(0,7,-1.2,4)%坐标轴的限制 学习内容 l多项式运算 l线性方程组求解 l非线性方程数值求解 l常微分方程初值问题的数值解法 l函数极值 l线性插值 线性方程组解的结构 n

7、一、齐次线性方程组解的结构 n齐次线性方程组的矩阵形式为AX=0 n其中A是mn阶矩阵；X为未知向量 二、线性方程组求解 1.直接解法 （1）利用左除运算符的直接解法 对于线性方程组Ax=b，可以利用左除运算符“”求 解: x=Ab 例:用直接解法求解下列线性方程组。 命令如下： A=2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4; b=13,-9,6,0; x=Ab 二、线性方程组求解 (2)利用矩阵的分解求解线性方程组 矩阵分解是指根据一定的原理用某种算法将 一个矩阵分解成若干个矩阵的乘积。 常见的矩阵分解有:特征值分解,奇异值分解 ,LU分解,QR分解,Chol

8、esky分解,Schur分解 Hessenberg分解 等 下面将对他们一一介绍 a)特征值分解 n矩阵的特征值分解也调用函数eig,其调用格式 为 n(1)V,D=eig(X)%得到矩阵X的特征值对角矩 阵D和其列为对应特征值的特征向量矩阵V,矩 阵的特征值分解为XV=VD n(2) V,D=eig(X,nobalance)%此形式为关闭平 衡算法的求解方法.平衡算法对于某些问题可以 得到更高的准确度。 n(3) )V,D=eig(A,B)%对矩阵A和B做广义特征 值分解,即AV=BVD 实例演示-矩阵的特征值分 解 n单位阵的特征值分解 na=-149 -50 -154; 537 180

9、546;-27 -9 -25 nv,d=eig(a) n双矩阵的特征值分解 nb=2 10 2;10 5 -8;2 -8 11; nv,d=eig(a,b) 附:复数特征值对角阵与实数特征 值对角阵的转换 n即使对于实阵,其特征值也可能出现复数.而在 实际使用中,常需要把这些共轭复数特征值转 换为实数块. nMATLAB提供两个函数进行转换 n(1)V,D=cdf2rdf(V,D)%将复数对角型转换 为实数块对角型 n(2)U,T=rdf2cdf(U,T)%将实数对角块转换为 复数对角型 实例演示-复数对角型矩阵 与实数对角型矩阵的转换 na=1 -3;2 2/3; nv,d=eig(a)%将

10、矩阵a进行特征值分解,其结 果复对角型 nvs,ds=cdf2rdf(v,d)%将复对角型转换为实对 角型 nb=vs*ds/vs%其结果与a相等 b)奇异值分解 n设矩阵A为一个mn阶的(实)矩阵(mn),则存在正交 矩阵V和U,使得: 此式称为A的奇异值分 解 n矩阵奇异值分解可由函数svd实现,其调用格式为 n(1)U,S,V=svd(X) n(2)U,S,V=svd(X,0) n作用:生成U,S,和V使得X=USV 实例演示 例：对矩阵a进行奇异值分解 a=1;1;U,S,V=svd(a) 注：进行奇异值分解时,矩阵的维数没有限制 c) LU分解 LU分解是高斯消去法的基础,它将一矩阵

11、表示为一个 交换下三角矩阵和一个上三角矩阵的乘积形式.线性代 数中已经证明,只要方阵A是非奇异的,LU分解总是可 以进行的。 lu函数用于对矩阵进行LU分解，其调用格式为 (1)L,U=lu(X)-产生一个上三角阵U和一个变换形式 的下三角阵L(行交换)，使之满足X=LU。注意,这里的 矩阵X必须是方阵。 (2)L,U,P=lu(X)-产生一个上三角阵U和一个下三角 阵L以及一个置换矩阵P，使之满足PX=LU。当然矩 阵X同样必须是方阵。 实现LU分解后，线性方程组Ax=b的解x=U(Lb)或 x=U(LPb)，这样可以大大提高运算速度。 实例演示 例 :用LU分解求解例1中的线性方程组。 命

12、令如下： A=2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4; b=13,-9,6,0; L,U=lu(A); X=U(Lb) 或采用LU分解的第2种格式，命令如下： L,U ,P=lu(A); X=U(LP*b) d)QR分解 n 在求解矩阵的特征值时,引入一种分解方法,即实 阵A可以写成A=QR,其中Q为一个正交矩阵和R 为一个上三角矩阵.QR分解只能对方阵进行. n QR分解由函数qr实现，其调用格式为: (1)Q,R=qr(X):产生一个一个正交矩阵Q和一个 上三角矩阵R,使之满足X=QR。 (2)Q,R,E=qr(X):产生一个一个正交矩阵Q,一个 上三角矩

13、阵R以及一个置换矩阵E,使之满足 XE=QR。 n 实现QR分解后，线性方程组Ax=b的解x=R(Qb) 或x=E(R(Qb)。 实例演示 例: 用QR分解求解前例中的线性方程组。 命令如下： A=2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4; b=13,-9,6,0; Q,R=qr(A); x=R(Qb) 或采用QR分解的第2种格式，命令如下： Q,R,E=qr(A); x=E*(R(Qb) e)Cholesky分解 n 如果矩阵X是对称正定的,则Cholesky分解将矩阵X分解成一 个对角元素为正的下三角矩阵R和上三角矩阵R的乘积. 使得 X=RR. n 函数c

14、hol(X)用于对矩阵X进行Cholesky分解,其调用格式为: (1)R=chol(X):产生一个上三角阵R,使RR=X.若X为非对称正定, 则输出一个出错信息。 (2)R,p=chol(X):这个命令格式将不输出出错信息.当X为对称正 定的,则p=0,R与上述格式得到的结果相同;否则p为一个正整 数.如果X为满秩矩阵,则R为一个阶数为q=p-1阶的上三角阵, 且满足RR=X(1:q,1:q). n 实现Cholesky分解后,线性方程组Ax=b变成RRx=b,所以 x=R(Rb)。 实例演示 例:用Cholesky分解求解前例中的线性方程组。 命令如下：A=2,1,-5,1;1,-5,0,

15、7;0,2,1,-1;1,6,-1,- 4; b=13,-9,6,0; R=chol(A) 结果为? Error using = chol Matrix must be positive definite%说明A为非正 定矩阵 ?可以执行命令R,p=chol(X) x=R(Rb) 2.迭代解法 n 迭代解法非常适合求解大型系数矩阵的方程组.在数 值分析中,迭代解法主要包括 :Jacobi迭代法、Gauss- Serdel迭代法、超松弛迭代法和两步迭代法。 1Jacobi迭代法:对于线性方程组Ax=b,如果A为非奇异 方阵,即aii0(i=1,2,n),则可将A分解为A=D-L-U， 其中D为对

16、角阵,其元素为A的对角元素,L与U为A的 下三角阵和上三角阵,于是AX=b化为:X=D- 1(L+U)x+D-1b 与之对应的迭代公式为:X(k+1)=D-1(L+U)x(k)+D-1b n 这就是Jacobi迭代公式.如果序列x(k+1)收敛于x， 则x必是方程Ax=b的解 附：Jacobi迭代法的 MATLAB函数文件 Jacobi.m如下： function y,n=jacobi(A,b,x0,eps) if nargin=3 eps=1.0e-6; elseif nargin=eps x0=y; y=B*x0+f; n=n+1; end 例: 用Jacobi迭代法求解下列线性方程组。设