《刘卫国全套配套课件matlab程序设计与应用第3版第8章matlab方程求解》由会员分享，可在线阅读，更多相关《刘卫国全套配套课件matlab程序设计与应用第3版第8章matlab方程求解（56页珍藏版）》请在金锄头文库上搜索。
1、第8章 MATLAB方程求解 8.1 线性方程组求解 8.2 非线性方程数值求解 8.3 最优化问题求解 8.4 常微分方程初值问题的数值求解,第8章 MATLAB方程求解,8.1 线性方程组求解,在MATLAB中，关于线性方程组的解法一般分为两类： 一类是直接法，就是在没有舍入误差的情况下，通过有限步的矩阵初等运算来求得方程组的解； 另一类是迭代法，就是先给定一个解的初始值，然后按照一定的迭代算法进行逐步逼近，求出更精确的近似解。,第8章 MATLAB方程求解,8.1.1 线性方程组的直接解法,在MATLAB中，这些算法已经被编成了现成的库函数或运算符，因此，只需调用相应的函数或运算符即可完
2、成线性方程组的求解。 1利用左除运算符的直接解法 线性方程组求解最简单的方法就是使用左除运算符“”，系统会自动根据输入的系数矩阵判断选用哪种方法进行求解。对于线性方程组Ax=b，可以利用左除运算符“”求解：x=Ab。,第8章 MATLAB方程求解,例8-1 用直接解法求解下列线性方程组。 A=2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4; b=13,-9,6,0; x=Ab 程序运行结果为： x = -66.5556 25.6667 -18.7778 26.5556,第8章 MATLAB方程求解,2利用矩阵的分解求解线性方程组 （1）LU分解 MATLAB提供的lu
3、函数用于对矩阵进行LU分解，其调用格式为： L,U=lu(X)：产生一个上三角阵U和一个变换形式的下三角阵L（行交换），使之满足X=LU。注意，这里的矩阵X必须是方阵。 L,U,P=lu(X)：产生一个上三角阵U和一个下三角阵L以及一个置换矩阵P，使之满足PX=LU。当然矩阵X同样必须是方阵。 当使用第一种格式时，矩阵L往往不是一个下三角阵，但可以通过行交换成为一个下三角阵。,第8章 MATLAB方程求解,设： 则对矩阵A进行LU分解的命令如下： A=1,-1,1;5,-4,3;2,1,1; L,U=lu(A) L = 0.2000 -0.0769 1.0000 1.0000 0 0 0.40
4、00 1.0000 0 U = 5.0000 -4.0000 3.0000 0 2.6000 -0.2000 0 0 0.3846 为检验结果是否正确，输入命令： LU=L*U LU = 1 -1 1 5 -4 3 2 1 1,第8章 MATLAB方程求解,利用第二种格式对矩阵A进行LU分解。 L,U,P=lu(A) L = 1.0000 0 0 0.4000 1.0000 0 0.2000 -0.0769 1.0000 U = 5.0000 -4.0000 3.0000 0 2.6000 -0.2000 0 0 0.3846 P = 0 1 0 0 0 1 1 0 0 LU=L*U %这种分
5、解其乘积不为A LU = 5 -4 3 2 1 1 1 -1 1 inv(P)*L*U %考虑矩阵P后的结果 ans = 1 -1 1 5 -4 3 2 1 1 实现LU分解后，线性方程组Ax=b的解x=U(Lb)或x=U(LPb)，这样可以大大提高运算速度。,第8章 MATLAB方程求解,例8-2 用LU分解求解例8-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(Lb) 程序运行结果为： x = -66.5556 25.6667 -18.7778 26.5556 或采用LU分解的第二
6、种格式，命令如下： L,U,P=lu(A); x=U(LP*b),第8章 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。,第8章 MATLAB方程求解,设： 则对矩阵A进行QR分解的命令如下： A=1,-1,1;5,-4,3;2,7,10; Q,R=qr(A) Q =
7、 -0.1826 -0.0956 -0.9785 -0.9129 -0.3532 0.2048 -0.3651 0.9307 -0.0228 R = -5.4772 1.2780 -6.5727 0 8.0229 8.1517 0 0 -0.5917 QR=Q*R QR = 1.0000 -1.0000 1.0000 5.0000 -4.0000 3.0000 2.0000 7.0000 10.0000,第8章 MATLAB方程求解,利用第二种格式对矩阵A进行QR分解： Q,R,E=qr(A) Q = -0.0953 -0.2514 -0.9632 -0.2860 -0.9199 0.2684
8、 -0.9535 0.3011 0.0158 R = -10.4881 -5.4347 -3.4325 0 6.0385 -4.2485 0 0 0.4105 E = 0 0 1 0 1 0 1 0 0 Q*R/E %验证A=Q*R*inv(E) ans = 1.0000 -1.0000 1.0000 5.0000 -4.0000 3.0000 2.0000 7.0000 10.0000 实现QR分解后，线性方程组Ax=b的解x=R(Qb)或x=E(R(Qb)。,第8章 MATLAB方程求解,例8-3 用QR分解求解例8-1中的线性方程组。 程序如下： A=2,1,-5,1;1,-5,0,7;
9、0,2,1,-1;1,6,-1,-4; b=13,-9,6,0; Q,R=qr(A); x=R(Qb) 程序运行结果为： x = -66.5556 25.6667 -18.7778 26.5556 或采用QR分解的第二种格式，命令如下： Q,R,E=qr(A); x=E*(R(Qb) 将得到与上面同样的结果。,第8章 MATLAB方程求解,（3）Cholesky分解 如果矩阵X是对称正定的，则Cholesky分解将矩阵X分解成一个下三角阵和一个上三角阵的乘积。设上三角阵为R，则下三角阵为其转置，即X=RR。MATLAB函数chol(X)用于对矩阵X进行Cholesky分解，其调用格式为： R=
10、chol(X)：产生一个上三角阵R，使RR=X。若X为非对称正定，则输出一个出错信息。 R,p=chol(X)：这个命令格式将不输出出错信息。若X为对称正定的，则p=0，R与上述格式得到的结果相同，否则p为一个正整数。如果X为满秩矩阵，则R为一个阶数为q=p-1的上三角阵，且满足RR=X(1:q,1:q)。,第8章 MATLAB方程求解,设： 则对矩阵A进行Cholesky分解的命令如下： A=2,1,1;1,2,-1;1,-1,3; R=chol(A) R = 1.4142 0.7071 0.7071 0 1.2247 -1.2247 0 0 1.0000 可以验证RR=A，命令如下： R*
11、R ans = 2.0000 1.0000 1.0000 1.0000 2.0000 -1.0000 1.0000 -1.0000 3.0000,第8章 MATLAB方程求解,利用第二种格式对矩阵A进行Cholesky分解： R,p=chol(A) R = 1.4142 0.7071 0.7071 0 1.2247 -1.2247 0 0 1.0000 p = 0 结果中p=0，这表示矩阵A是一个正定矩阵。如果试图对一个非正定矩阵进行Cholesky分解，则将得出错误信息，所以，chol函数还可以用来判定矩阵是否为正定矩阵。 实现Cholesky分解后，线性方程组Ax=b变成RRx=b，所以x
12、=R(Rb)。,第8章 MATLAB方程求解,例8-4 用Cholesky分解求解例8-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) 错误使用 chol 矩阵必须为正定矩阵。 命令执行时，出现错误信息，说明A为非正定矩阵。,第8章 MATLAB方程求解,8.1.2 线性方程组的迭代解法,迭代法是一种不断用变量的旧值递推新值的过程，是用计算机解决问题的一种基本方法。它利用计算机运算速度快、适合做重复性操作的特点，让一组指令重复执行，在每次执行这组指令时，都从变量的原值推出它的新值。 迭
13、代解法非常适合求解大型稀疏矩阵的方程组。在数值分析中，迭代解法主要包括 Jacobi迭代法、Gauss-Serdel迭代法、超松弛迭代法和两步迭代法。,第8章 MATLAB方程求解,迭代法的思想 :,将方程改写为：,这种形式的好处是将一组x代入右端，可以立即得到另一组x。如果两组x相等，那么它就是方程组的解，不等时可以继续迭代。可以构造方程的迭代公式为：,第8章 MATLAB方程求解,1Jacobi迭代法 对于线性方程组Ax=b，如果A为非奇异方阵，且aii0（i=1，2，n），则可将A分解为A=D-L-U，其中D为对角阵，其元素为A的对角元素，L与U为A的下三角阵取反和上三角阵取反：,于是A
14、x=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的解。,第8章 MATLAB方程求解,Jacobi迭代法的MATLAB函数文件jacobi.m如下： function y,n=jacobi(A,b,x0,ep) if nargin=3 ep=1.0e-6; elseif nargin=ep x0=y; y=B*x0+f; n=n+1; end,第8章 MATLAB方程求解,例8-5 用Jacobi迭代法求解下列线性方程组。设迭代初值为0，迭
15、代精度为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) 程序运行结果为： x = 0.9958 0.9579 0.7916 n = 11,第8章 MATLAB方程求解,2Gauss-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迭代用新分量代替旧分量，精度会高些。,第8章 MATLAB方程求解,Gauss-Serdel迭代法的MATLAB函数文件gauseidel.m如下： function y,n=gauseidel(A,b,x0,ep) if nargin=3 e
