第7章 MATLAB解方程与函数极值

MATLAB R2006a 基础教程

第7章 MATLAB解方程与函数极值

第7章

MATLAB解方程与函数极值

7.1 线性方程组求解7.2 非线性方程数值求解 7.3 常微分方程初值问题的数值解法 7.4 函数极值

中国石油大学(北京)电子信息工程系

MATLAB R2006a 基础教程

第7章 MATLAB解方程与函数极值

7.1 线性方程组求解 7.1.1 直接解法1.利用左除运算符的直接解法

对于线性方程组Ax=b,可以利用左除运算符“\” 求解： x=A\b

中国石油大学(北京)电子信息工程系

MATLAB R2006a 基础教程

第7章 MATLAB解方程与函数极值

例7-1 用直接解法求解下列线性方程组。命令如下：

A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4]; b=[13,-9,6,0]'; x=A\b

中国石油大学(北京)电子信息工程系

MATLAB R2006a 基础教程

第7章 MATLAB解方程与函数极值

2.利用矩阵的分解求解线性方程组

矩阵分解是指根据一定的原理用某种算法将一个 矩阵分解成若干个矩阵的乘积。常见的矩阵分解 有LU分解、QR分解、Cholesky分解，以及Schur分 解、Hessenberg分解、奇异分解等。

中国石油大学(北京)电子信息工程系

MATLAB R2006a 基础教程

第7章 MATLAB解方程与函数极值

(1) LU分解 矩阵的LU分解就是将一个矩阵表示为一个交换下 三角矩阵和一个上三角矩阵的乘积形式。线性代 数中已经证明，只要方阵A是非奇异的，LU分解总 是可以进行的。 MATLAB提供的lu函数用于对矩阵进行LU分解，其 调用格式为： [L,U]=lu(X):产生一个上三角阵U和一个变换形 式的下三角阵L(行交换),使之满足X=LU。注意， 这里的矩阵X必须是方阵。 [L,U,P]=lu(X):产生一个上三角阵U和一个下三 角阵L以及一个置换矩阵P,使之满足PX=LU。当然 矩阵X同样必须是方阵。 实现LU分解后，线性方程组Ax=b的解x=U\(L\b)或 x=U\(L\Pb),这样可以大大提高运算速度。中国石油大学(北京)电子信息工程系

MATLAB R2006a 基础教程

第7章 MATLAB解方程与函数极值

例7-2 用LU分解求解例7-1中的线性方程组。命令如下：

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\(L\b)或采用LU分解的第2种格式，命令如下：

[L,U ,P]=lu(A); x=U\(L\P*b)

中国石油大学(北京)电子信息工程系

MATLAB R2006a 基础教程

第7章 MATLAB解方程与函数极值

(2) QR分解 对矩阵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))。中国石油大学(北京)电子信息工程系

MATLAB R2006a 基础教程

第7章 MATLAB解方程与函数极值

例7-3 用QR分解求解例7-1中的线性方程组。命令如下：

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\(Q\b)或采用QR分解的第2种格式，命令如下：

[Q,R,E]=qr(A); x=E*(R\(Q\b))

中国石油大学(北京)电子信息工程系

MATLAB R2006a 基础教程

第7章 MATLAB解方程与函数极值

(3) Cholesky分解 如果矩阵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)。中国石油大学(北京)电子信息工程系

MATLAB R2006a 基础教程

第7章 MATLAB解方程与函数极值

例7-4 用Cholesky分解求解例7-1中的线性方程组。命令如下：

A=[2,1,-5,1;1,-5,0,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为非正定矩阵 。中国石油大学(北京)电子信息工程系

MATLAB R2006a 基础教程

第7章 MATLAB解方程与函数极值

7.1.2 迭代解法迭代解法非常适合求解大型系数矩阵的方程组。在数值分 析中，迭代解法主要包括 Jacobi迭代法、Gauss-Serdel 迭代法、超松弛迭代法和两步迭代法。 1.Jacobi迭代法 对于线性方程组Ax=b,如果A为非奇异方阵，即 aii≠0(i=1,2,…,n),则可将A分解为A=D-L-U,其中D为 对角阵，其元素为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 这就是Jacobi迭代公式。如果序列{x(k+1)}收敛于x,则x 必是方程Ax=b的解。中国石油大学(北京)电子信息工程系

MATLAB R2006a 基础教程

第7章 MATLAB解方程与函数极值 Jacobi迭代法的MATLAB函数文件Jacobi.m如下：function [y,n]=jacobi(A,b,x0,eps) if nargin==3 eps=1.0e-6; elseif nargin<3 error return end D=diag(diag(A)); %求A的对角矩阵 L=-tril(A,-1); %求A的下三角阵 U=-triu(A,1); %求A的上三角阵 B=D\(L+U); f=D\b; y=B*x0+f; n=1; %迭代次数 while norm(y-x0)>=eps x0=y; y=B*x0+f; n=n+1; end中国石油大学(北京)电子信息工程系

MATLAB R2006a 基础教程

第7章 MATLAB解方程与函数极值

例7-5

用Jacobi迭代法求解下列线性方程组。设迭代初值为0,迭代精度为10-6。在命令中调用函数文件Jacobi.m,命令如下：

A=[10,-1,0;-1,10,-2;0,-2,10]; b=[9,7,6]'; [x,n]=jacobi(A,b,[0,0,0]',1.0e-6)

中国石油大学(北京)电子信息工程系

MATLAB R2006a 基础教程

第7章 MATLAB解方程与函数极值

2.Gauss-Serdel迭代法 在Jacobi迭代过程中，计算时，已经得到，不必再用，即 原来的迭代公式Dx(k+1)=(L+U)x(k)+b可以改进为 Dx(k+1)=Lx(k+1)+Ux(k)+b,于是得到： x(k+1)=(D-L)-1Ux(k)+(D-L)-1b 该式即为Gauss-Serdel迭代公式。和Jacobi迭代相比， Gauss-Serdel迭代用新分量代替旧分量，精度会高些。

中国石油大学(北京)电子信息工程系

MATLAB R2006a 基础教程

第7章 MATLAB解方程与函数极值 Gauss-Serdel迭代法的MATLAB函数文件gauseidel.m如下：function [y,n]=gauseidel(A,b,x0,eps) if nargin==3 eps=1.0e-6; elseif nargin<3 error return end D=diag(diag(A)); %求A的对角矩阵 L=-tril(A,-1); %求A的下三角阵 U=-triu(A,1); %求A的上三角阵 G=(D-L)\U; f=(D-L)\b; y=G*x0+f; n=1; %迭代次数 while norm(y-x0)>=eps x0=y; y=G*x0+f; n=n+1; end中国石油大学(北京)电子信息工程系

MATLAB R2006a 基础教程

第7章 MATLAB解方程与函数极值

例7-6 用Gauss-Serdel迭代法求解下列线性方程组。设迭代初值为0,迭代精度为10-6。在命令中调用函数文件gauseide …… 此处隐藏：3228字，全部文档内容请下载后查看。喜欢就下载吧 ……