matlab第六讲

MATLAB的方程(组)解法第六讲王文健[e-mail]《MATLAB数据处理与应用》2011-2012学年选修课西南交通大学摩擦学研究所TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY2主要内容线性方程(组)的解法非线性方程(组)解法MATLAB统计分析西南交通大学摩擦学研究所TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY3简介线性方程(组)直接法—在没有舍入的情况下，通过有限步四则运算求得方程组的准确解迭代法—先给出一个解的初始值，然后按一定的法则逐步求出解的各个更准确的近似值的方法非线性方程(组)迭代法—不动点迭代法、Newton迭代法、BroydenBroyden迭代法西南交通大学摩擦学研究所TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY4线性方程(组)的解法线性方程的解法对于线性方程可以直接调用roots函数来求解例如：求解x3+5x2-8x+6=0P=[15-86];roots(P)求解2x9+43x7+x6-8x5+14x4-5x3+x2-10x+12=0a=[20431-814-51-1012];roots(a)西南交通大学摩擦学研究所TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY5线性方程(组)的解法线性方程组的解法直接法—利用符号运算“/”或“\”来求解直接法一般基于高斯消去法、主元素消去法、平方根法和追赶法等，在MATLAB中这些算法已被编成现成的库函数或运算符，可以直接调用进行求解067452296385243214324214321xxxxxxxxxxxxxx西南交通大学摩擦学研究所TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY6线性方程(组)的解法线性方程组的解法建立系数矩阵和常数矩阵A=[21-51;1-30-6;02-12;14-76]B=[89-50]'X=A\B西南交通大学摩擦学研究所TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY7线性方程(组)的解法线性方程组的解法利用矩阵的LU、QR和Cholesky分解法求解—这三种方法对求解大型方程组非常有用优点是运算速度快、可以节省磁盘空间、节省内存LU分解法—Gauss消去分解，可以把任意矩阵分解为下三角矩阵的基本变换形式和上三角矩阵的乘积A=LUL为下三角矩阵，U为上三角矩阵A*X=b变成L*U*X=b—X=U\(L\b)西南交通大学摩擦学研究所TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY8线性方程(组)的解法LU分解法A=[21-51;1-30-6;02-12;14-76]B=[89-50]'[L,U]=lu(A)X=U\(L\B)西南交通大学摩擦学研究所TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY9线性方程(组)的解法Cholesky分解法如果A为对称正定矩阵，则Cholesky分解可将矩阵A分解为上三角矩阵和其转置的乘积，即A=R'*R，其中R上三角矩阵A*X=b变成R'*R*X=bX=R\(R'\b)举例：A=[1648;45-4;8-422]B=[28526]'R=chol(A)X=R\(R'\B)西南交通大学摩擦学研究所TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY10线性方程(组)的解法QR分解法对于任何长方矩阵A，都可以进行QR分解，其中Q为正交矩阵，R为上三角矩阵，即A=QRA*X=b变成Q*R*X=bX=R\(Q\b)举例：A=[1648;45-4;8-422]B=[28526]'[Q,R]=qr(A)X=R\(Q\B)西南交通大学摩擦学研究所TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY11线性方程(组)的解法迭代解法在MATLAB中迭代法非常适合求解大型系数矩阵的方程组Jacobi迭代法Gauss-Serdel迭代法超松弛迭代法两步迭代法西南交通大学摩擦学研究所TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY12线性方程(组)的解法Jacobi迭代法线性方程组Ax=b，如果系数矩阵A为非奇异矩阵，则A可写成A=D.L.U的分解形式，其中D为对角矩阵，元素为A的对角元素，L和U分别为严格的下三角矩阵和上三角矩阵，故迭代形式为：xi+1=D-1(L+U)xi+D-1b对应的Jacobi迭代公式的矩阵为：x(k+1)=Bx(k)+f式中：B=D-1(L+U)=I-D-1Af=D-1b西南交通大学摩擦学研究所TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY13线性方程(组)的解法Jacobi迭代法举例：A=[10,-1,0;-110-2;0-210];b=[9;7;6];jacobi(A,b,[0;0;1])610272109103232121xxxxxxx西南交通大学摩擦学研究所TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY14线性方程(组)的解法Gauss-Seidel迭代法线性方程组Ax=b，Gauss-Seidel迭代形式为：X(k+1)=GX(k)+f式中：G为Gauss-Seidel迭代矩阵，G=-(D+L)-1b，f==(D+L)-1b，D为对角矩阵，L和U为严格的下三角和上三角矩阵编制Gauss-Seidel迭代法的M文件西南交通大学摩擦学研究所TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY15线性方程(组)的解法Gauss-Seidel迭代法举例：A=[10,-1,0;-110-2;0-210];b=[9;7;6];GS(A,b,[0;0;0])610272109103232121xxxxxxx西南交通大学摩擦学研究所TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY16线性方程(组)的解法SOR迭代法线性方程组Ax=b，SOR迭代形式为：X(k+1)=B0X(k)+f式中：B0为SOR迭代矩阵，B0=(D-wL)-1[(1-w)D+wU],f=w(D-wL)-1bD为对角矩阵，L和U为严格的下三角和上三角矩阵编制SOR迭代法的M文件西南交通大学摩擦学研究所TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY17线性方程(组)的解法SOR迭代法举例：A=[43-1;34-1;1-14];b=[19;30;27];w=1.03;SOR(A,b,[1;1;1],w)27430431934321321321xxxxxxxxx西南交通大学摩擦学研究所TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY18线性方程(组)的解法双步迭代法线性方程组Ax=b，如果A为对称正定矩阵，双步迭代法形式为：(D-L)x(k+1/2)=Ux(k)+b(D-U)x(k+1)=Lx(k+1/2)+b式中：D为对角阵，L和U为严格的下三角和上三角矩阵编制SOR迭代法的M文件西南交通大学摩擦学研究所TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY19线性方程(组)的解法双步迭代法举例：A=[10-122;-111-13;2-1103;23-18];b=[8;25;-11;17];TS_D(A,b,[0;0;0;0])1783321110225311822104321432143214321xxxxxxxxxxxxxxxx西南交通大学摩擦学研究所TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY20非线性方程的解法对比法假定函数在[a,b]区间上是单调函数，且f(a)f(b)0，则在区间上方程f(x)=0至少有一个根思想：通过判断函数f(x)的符号，逐步将有限区间缩小，使在足够小的区间内方程有且仅有一个根，近似根就是最终区间的重点位置举例：求解非线性方程xex-2=01.编制函数文件exam.m2.命令窗口输入：DF('exam',0,10,0.05)西南交通大学摩擦学研究所TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY21非线性方程的解法迭代法不定点迭代法思想：非线性方程f(x)=0，如果在区间上连续，可将函数表示为x=g(x)，迭代格式为：xk+1=g(xk)，如果limxk=x*，则迭代过程将收敛于方程的根举例：求解非线性方程ex-3x2+3=01.编制函数文件fun1.m2.命令窗口输入：[x,y,n]=iterate(10)西南交通大学摩擦学研究所TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY22非线性方程的解法迭代法切线迭代法—牛顿迭代法思想：迭代格式xk+1=xk-f(xk)/fg(xk)-给出导数举例：求解非线性方程ex-3x2=01.编制函数文件ff.m2.编制导数文件dff.m3.命令窗口输入：[x,n]=NT(2)西南交通大学摩擦学研究所TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY23非线性方程的解法迭代法割线迭代法思想：当不存在解析的导函数时，割线法可近似来求解导函数，迭代格式xk+1=xk-f(xk)(xk-xk-1)/(f(xk)-f(xk-1))举例：求解非线性方程ex-3x2=01.编制函数文件GX.m2.编制函数文件ff.m3.命令窗口输入：y=GX(2,1)西南交通大学摩擦学研究所TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY24非线性方程的解法求解非线性方程的MATALB函数求解非线性方程f(x)=0的函数fsolve调用格式：x=fsolve(fun,x0)—变量fun为方程表达式或函数文件名，x0为估计的方程根的初始值x=fsolve(fun,x0,options)—根据优化参数options进行最小化，默认值为0，如为1显示运算的中间步骤[x,fval]=fsolve(fun,x0)—返回值为目标函数在求得的根x对应的函数值西南交通大学摩擦学研究所TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY25非线性方程的解法求解非线性方程的MATALB函数求解非线性方程f(x)=0的函数fsolve调用格式：[x,fval,exitflag]=fsolve(…)—返回值为函数调用的退出标记exitflag，用于描述当前退出状态[x,fval,exitflag,output]=fsolve(…)—函数返回值为一个结构输出，用于描述优化的信息[x,fval,exitflag,output,jacobian]=fsolve(…)—函数返回值为在点x处函数的Jacobian数值西南交通大学摩擦学研究所TribologyRes