《MATLAB程序设计教程(第二版)》第6章__MATLAB解方程与最优化问

《MATLAB程序设计教程(第二版)》第6章__MATLAB解方程与最优化问题求解

MATLAB程序设计教程(第二版) 程序设计教程(第二版) 程序设计教程刘卫国 主编

中国水利水电出版社

《MATLAB程序设计教程(第二版)》第6章__MATLAB解方程与最优化问题求解

第6章 MATLAB解方程与最优化问题求解

MATLAB线性方程组求解 MATLAB非线性方程数值求解 MATLAB常微分方程初值问题的数值解法 MATLAB最优化问题求解

《MATLAB程序设计教程(第二版)》第6章__MATLAB解方程与最优化问题求解

6.1 线性方程组求解 6.1.1 直接解法 1.利用左除运算符的直接解法 对于线性方程组Ax=b,可以利用左除运算符“\”求解： x=A\b

《MATLAB程序设计教程(第二版)》第6章__MATLAB解方程与最优化问题求解

例6-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程序设计教程(第二版)》第6章__MATLAB解方程与最优化问题求解

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

《MATLAB程序设计教程(第二版)》第6章__MATLAB解方程与最优化问题求解

(1) LU分解 矩阵的LU分解就是将一个矩阵表示为一个交换下 三角矩阵和一个上三角矩阵的乘积形式。线性 代数中已经证明，只要方阵A是非奇异的，LU 分解总是可以进行的。

《MATLAB程序设计教程(第二版)》第6章__MATLAB解方程与最优化问题求解

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程序设计教程(第二版)》第6章__MATLAB解方程与最优化问题求解

例6-2 用LU分解求解例6-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程序设计教程(第二版)》第6章__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程序设计教程(第二版)》第6章__MATLAB解方程与最优化问题求解

例6-3 用QR分解求解例6-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程序设计教程(第二版)》第6章__MATLAB解方程与最优化问题求解

(3) Cholesky分解 如果矩阵X是对称正定的，则Cholesky分解将矩阵X分 解成一个下三角矩阵和上三角矩阵的乘积。设上三 角矩阵为R,则下三角矩阵为其转置，即X=R'R。

《MATLAB程序设计教程(第二版)》第6章__MATLAB解方程与最优化问题求解

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程序设计教程(第二版)》第6章__MATLAB解方程与最优化问题求解

例6-4 用Cholesky分解求解例6-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程序设计教程(第二版)》第6章__MATLAB解方程与最优化问题求解

6.1.2 迭代解法 迭代解法非常适合求解大型系数矩阵的方程组。在数 值分析中，迭代解法主要包括 Jacobi迭代法、 Gauss-Serdel迭代法、超松弛迭代法和两步迭代法。

《MATLAB程序设计教程(第二版)》第6章__MATLAB解方程与最优化问题求解

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程序设计教程(第二版)》第6章__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