第五章

第五章数据分析§5.1线性方程组5.1.1线性方程求解在矩阵的表示方法中，线性方程的求解可以表述为：给定两个矩阵A和B，求X的唯一解使得：AX=B或XA=BX=A\B:表示求矩阵方程AX=B的解X=B/A:表示求矩阵方程XA=B的解(B/A)’=(A’\B’)矩阵A并不要求是方阵，若A是m*n的矩阵，则存在以下三种情况：1)m=n:适定方程组，寻求精确解；2)mn:超定方程组，寻求最小二乘解;3)mn:不定方程组，寻求基本解.对于不同的情况，左除算子采用不同的算法求解。5.1.1线性方程求解5.1.1线性方程求解1、适定方程组：一般形式的线性方程组为：Ax=b或AX=B求解命令是：x=A\bX=A\B2、超定方程组：对于超定方程，在MATLAB中，利用左除命令寻求它的最小二乘解，其调用格式是：x=A\bA*x并不是精确的等于b，但它们的差值小于原始数据的测量误差。3、不定方程组：解不唯一，MATLAB将寻求一个基本解，其中至多只有m个非零元素。5.1.2矩阵分解MATLAB中几种常用的矩阵分解方法，相关的函数见下表：命令功能说明cholcholesky分解对称正定矩阵的分解，用于求解方程组lu矩阵LU分解将矩阵采用LU分解，直接求解线性方程组qr正交三角分解将矩阵分解为正交矩阵或酉矩阵和上三角矩阵5.1.2矩阵分解Cholesky分解的函数chol的调用格式为：R=chol(X):X是对称正定矩阵，R是上三角阵(对角线为正)，使得R’*R=X，若X是非正定的，将给出错误信息。[R,p]=chol(X):不会给出错误信息，若X是正定的，p等于0，R同上；若X不是正定的，则p为正整数，R是上三角矩阵，它的阶数为q=p-1，Cholesky分解允许对线性方程组A*x=b进行如下替换：R’*R*x=b由于左除算子可以处理三角矩阵，因此可以得出：x=R\(R’\b)一、Cholesky分解5.1.2矩阵分解二、LU分解函数lu的调用格式：[L,U]=lu(A):U是上三角矩阵，L是“心理上”的下三角阵，实际上，它是下三角阵和置换矩阵的乘积，结果使得：A=L*U[L,U,P]=lu(A):L是下三角矩阵，上三角矩阵U和置换矩阵P，使得：P*A=L*ULU分解允许对线性方程A*x=b进行如下计算：x=U\(L\b)LU分解或Gaussian消去法，可以将任何方阵表示为一个下三角矩阵L和上三角矩阵U的乘积,即：A=LU5.1.2矩阵分解三、正交分解正交分解或QR分解，将方阵或矩形矩阵A分解为一个正交矩阵Q和一个上三角矩阵R的乘积，即：A=QR或AP=QR函数qr的用法如下：[Q,R]=qr(X):上三角矩阵R和正交矩阵Q，使得：X=Q*R(R与X同维)[Q,R,E]=qr(X):置换矩阵E，使得：X*E=Q*R[Q,R]=qr(X,0)和[Q,R,E]=qr(X,0)：其中，E是一个置换向量，使得：Q*R=X(:,E)。§5.2非线性数值计算5.2.1非线性函数最小值点1、求单变量函数最小值点求单变量函数最小值点的函数是:fminbnd，其调用格式是：[x,fval,exitflag,output]=fminbnd(fun,x1,x2,options,p1,p2,…)fun是被计算最小值点的单变量函数（目标函数）名称字符串；x1、x2是目标函数自变量的取值范围；p1、p2,…是向目标函数传递的附加参数；options是一个结构类型的变量，用于指定算法的优化参数。该结构的内容用函数optimset进行定义，若没有优化参数要设置，可以在options的位置用“[]”作为占位符。5.2.1非线性函数最小值点1、求单变量函数最小值点Display:显示的层次，off不显示输出内容，iter显示每次迭代的输出，final只显示最后一次的输出；MaxFunEvals:所允许的函数的最大求值次数；MaxIter:所允许的最大迭代数；TolX:终止迭代的容差值函数fminbnd使用了结构options的四个域，它们的含义如下：[x,fval,exitflag,output]=fminbnd(fun,x1,x2,options,p1,p2,…)输出参数：x是fun的最小值点;fval是目标函数在x处的值;exitflag描述函数fminbnd退出的情况，大于0表示函数收敛到了解x，等于0表示达到了所允许的函数最大求值次数。output是一个结构变量，包含了计算的信息，它共有三个域：5.2.1非线性函数最小值点output.algorithm:使用的算法；output.funcCount:函数求值次数；output.iterations:迭代次数。[x,fval,exitflag,output]=fminbnd(fun,x1,x2,options,p1,p2,…)例：求内联函数f(x,a)=x4+2x2+2ax+5的最小值。5.2.1非线性函数最小值点函数fminsearch用于求多变量函数的最小值点，其调用格式是：[x,fval,exitflag,output]=fminsearch(fun,x0,options,p1,p2,…)其中，x0是最小值点的初始值，一般是猜测的；结构options的域比函数fminbnd中增加一个---TolFun,指终止计算的目标函数容差值;退出标志exitflag也比函数fminbnd中增加一种情况：小于0，表示目标函数不收敛。2、求多变量函数最小值点例：求函数f(x)=100(x2-x12)2+(a-x1)2的最小值。5.2.2单变量函数零点求单变量函数零点的函数是fzero，其调用格式是：[x,fval,exitflag,output]=fzero(fun,x0,options,p1,p2,…)其中，参数x0指定搜索的起始点，同时也预定搜索零点的大致位置；结构options中只有两个域：Display和TolX输出参数exitflag大于0表示找到了函数的零点，小于0表示没有找到零点或在搜索过程中遇到了无穷大的函数值。例：求函数f(x)=x4-ax-5a+7的零点5.2.3常微分方程（组）数值解刚性问题，对于一个常微分方程组，如果其Jacobian矩阵的特征值相差十分悬殊，那么这个方程组就称为刚性方程组。七种函数对应解法的特点：1)ode45:四阶/五阶龙格—库塔法，单步解法，对很多问题都是首先试用的最好方法，但不能解决刚性问题；2)ode23:二阶/三阶龙格—库塔法，单步解法，在误差容许范围较宽且存在轻微刚性时比ode45效果好；3)ode113:可变阶次Adams-Bashforth-MoultonPECE算法，多步解法，比ode45更适合于误差容许范围要求较严格的情况，不能解决刚性问题.MATLAB提供了七种解法，分别由七个不同的函数来完成，它们是：ode45ode23ode113ode15sode23sode23tode23tb5.2.3常微分方程（组）数值解4)ode15s:可变阶次的数值微分公式算法，多步解法，若用ode45效果很差或者失败，可以考虑试用ode15s，且可解决刚性问题；5)ode23s:基于修正的Rosenbrock公式，单步解法，比ode15s更适用于误差容许范围较宽的情况，可以解决一些采用ode15s效果不好的刚性问题；6)ode23t:采用自由内插法实现的梯形规则，适用于系统存在适度刚性并要求无数值衰减的情况；7)ode23tb:龙格—库塔公式的第一级采用梯形规则、第二级采用Gear法。对于误差容限较宽的情况，比ode15s效果更好，可解决刚性问题。5.2.3常微分方程（组）数值解函数ode45的调用格式为：[T,Y]=ode45(‘F’,tspan,y0,options,p1,p2,…)输入参数：‘F’是常微分方程组的文件名字符串;tspan是向量，一般形如：[t0tfinal],指定求解的起始值和终止值，若需要得到指定时刻的结果，也可以采用:[t0,t1,…,tfinal];y0是初始状态向量；options是一些可选的综合参数，由函数odeset进行设置。输出参数：T为时间列向量，Y是数值解矩阵，它的每一行对应着列向量T中的一个时间值。例1：解常微分方程组，已知常微分方程组及初始条件如下：｛y1’=y2y3y1(0)=0y2’=-y1y3y2(0)=1y3’=-2y1y2y3(0)=-15.2.3常微分方程(组)数值解例2：解常微分方程组，已知二阶微分方程及初始条件：y’’=1000(1-y2)y’-yy(0)=0,y’(0)=1解：令y1=y,y2=y1’,则该方程化为如下的一阶常微分方程组：y1’=y2y2’=1000(1-y12)y2-y15.2.4函数的数值积分求解函数的数值积分的函数为quadl和quad，其调用格式为：quadl(fun,a,b,Tol,p1,p2,…)例1：在区间[0.01,1]上求的数值积分，设置绝对误差容限为1e-8。x例2：在区间[0,2]上求的数值积分。5213xx§5.3多项式5.3.1多项式表示法MATLAB采用行向量表示多项式，将多项式的系数按降幂次序存放在行向量中，多项式：p(x)=a0xn+a1xn-1+…+an-1x+an的系数行向量为：p=[a0a1…an],注意顺序必须是从高次幂向低次幂。为了将整个多项式输入MATLAB，可以采用类似于向量创建的方法，直接输入多项式的系数行向量。例如：多项式x4+3x2-5x+1可以采用系数向量[103–51]表示，注意多项式中缺少的幂次要用“0”补齐。5.3.2求多项式的根，由根创建多项式函数roots用于求多项式的根，其调用形式为：R=roots(C)函数poly由给定的根创建多项式，其调用格式为：C=poly(R)按照输入参数R的不同形式，该函数的输出参数C有以下两种情况：1)R是向量，那么C将是根向量为R的多项式系数行向量；2)R为N*N矩阵，那么返回值C将是矩阵R特征多项式的系数行向量，它由N+1个元素组成。例：求多项式x4-x3+9x-2的根5.3.3多项式乘和除函数conv用于实现向量的卷积运算，函数deconv用于实现向量的解卷运算。函数conv的调用格式是：C=conv(a,b)函数deconv的调用格式如下：[q,r]=deconv(c,a)例：多项式乘除运算：以多项式a(x)=x3-2x2+1和b(x)=3x2-4x+5为例5.3.4多项式导数函数polyder可计算多项式的导数,它既可以计算单个多项式的导数，又可以计算两个多项式积和商的导数。函数polyder的调用格式为：polyder(p):计算系数向量为p的多项式的导数；polyder(a,b)或c=polyder(a,b):计算多项式a*b的导数；[q,d]=polyder(a,b):计算多项式的商a/b的导数，用q/d表示。例：多项式a(x)=3x2+2x-7和b(x)=4x3-3x2+x-1的导数5.3.5多项式的值函数polyval用于计算给定的自变量值时，多项式的值函数polyval的调用格式：y=polyval(p,x)其中，p为多项式的系数向量，x是指定的自变量的值。函数polyvalm计算矩阵多项式的值，即以矩阵整体作为多项式的自变量进行计算。该函数用法为：y=polyvalm(p,x)例：求多项式的值，p(x)=5x2+4x+35.3.6多项式曲线拟合函数polyfit用于对给定的数据进行多项式曲线拟合，最后给出拟合多项式系数，此函数的调用格式为：p=polyfit(x,y,n)其中，x和y是数据的横纵坐标向量，n是规定的拟合多项式的阶数，p是返回多项式的系数向量。5.3.7部分分式展开函数residue用于将两个多项式之比用部分分式展开，该函数的调用格式是：[r,p,k]=residue(a,b):求多项式之比a/b的部分分式展开，r是留数、p是极点、k是直接项;当没有重根时，存在：当存在重根时，有：[a,b]=residue(r,p,k):从部分分式得出多项式系数向量。