MATLAB计算方法3解线性方程组计算解法

MATLAB计算方法3解线性方程组计算解法第一页，共73页。解线性方程组第三章线性方程组数值解法§3.1直接法第二页，共73页。一、Gauss消去法设有消元：第三页，共73页。用Matlab实现顺序Gauss消去法在Matlab程序编辑器中输入：functionx=nagauss(a,b,flag)%解线形方程组ax=b，a为系数矩阵，b为右端列向量，flag若为0，则显示中间过程，否则不显示，默认为0，x为解向量ifnargin<3,flag=0;endn=length(b);a=[a,b];%消元fork=1:(n-1)a((k+1):n,(k+1):(n+1))=a((k+1):n,(k+1):(n+1))-a((k+1):n,k)/a(k,k)*a(k,(k+1):(n+1));a((k+1):n,k)=zeros(n-k,1);ifflag==0,a,endend第四页，共73页。%回代x=zeros(n,1);x(n)=a(n,n+1)/a(n,n);fork=n-1:-1:1x(k,:)=(a(k,n+1)-a(k,(k+1):n)*x((k+1):n))/a(k,k);end第五页，共73页。第六页，共73页。程序运行结果：第七页，共73页。二列主元素Gauss消去法---计算结果可靠第八页，共73页。到此原方程组化为

第九页，共73页。第十页，共73页。(上三角方程组)（3.2)直到(n-1)原方程组化为以上为消元过程。第十一页，共73页。

（3.3)(n)回代求解公式

第十二页，共73页。系数矩阵为对称正定阵或严格对角占优阵的方程组按高斯消去法计算是数值稳定的，因而不必选主元。严格对角占优阵：至少有一个主对角线元素的绝对值严格大于此行或此列其他元素的绝对值之和。第十三页，共73页。Gauss列主元消去法：优点------计算结果更可靠；缺点------挑主元花机时更多,次序有变动，程序复杂。说明：(1)也可采用无回代的列主元消去法(叫Gauss---Jordan消去法)，但比有回代的列主元消去法的乘除运算次数多。(2)有回代的列主元消去法所进行的乘除运算次数为，量很小。第十四页，共73页。用Matlab实现选列主元Gauss消去法解线性方程组在Matlab程序编辑器中输入：functionx=nagauss2(a,b,flag)%a为系数矩阵；b为右端列向量；flag若为0，则显示中间过程，否则不显示ifnargin<3,flag=0;endn=length(b);a=[a,b];%选主元fork=1:(n-1)[ap,p]=max(abs(a(k:n,k)));p=p+k-1;ifp>k,t=a(k,:);a(k,:)=a(p,:);a(p,:)=t;end%消元

a((k+1):n,(k+1):(n+1))=a((k+1):n,(k+1):(n+1))-a((k+1):n,k)/a(k,k)*a(k,(k+1):(n+1));a((k+1):n,k)=zeros(n-k,1);

第十五页，共73页。ifflag==0,a,endend%回代x=zeros(n,1);x(n)=a(n,n+1)/a(n,n);fork=n-1:-1:1x(k,:)=(a(k,n+1)-a(k,(k+1):n)*x((k+1):n))/a(k,k);end第十六页，共73页。第十七页，共73页。程序运行结果:第十八页，共73页。

记，则

（三角因子分解）

Gauss消元，初等行变换，化原方程组为上三角型。三．矩阵三角分解法第十九页，共73页。定义3.1

叫的三角（因子）分解，其中是是上三角。下三角,为单位下三角阵（对角元全为1），为上三角阵，则称

为Doolittle分解;若是下三角，

是单位上三角，则称定义3.2

若

为Crout分解。定义3.2若第二十页，共73页。

为什么要讨论三角分解？若在消元法进行前能实现三角分解，则从而容易回代求解。第二十一页，共73页。1．直接三角分解法(Doolittle分解为例）

第二十二页，共73页。由矩阵乘法………第二十三页，共73页。

…….(k)第二十四页，共73页。第二十五页，共73页。第二十六页，共73页。例3．1第二十七页，共73页。用Matlab实现LU分解在Matlab程序编辑器中输入：function[L,U]=nalu(a)%a为可逆方阵；L返回单位下三角矩阵；U返回上三角矩阵n=length(a);U=zeros(n,n);L=eye(n,n);U(1,:)=a(1,:);L(2:n,1)=a(2:n,1)/U(1,1);fork=2:nU(k,k:n)=a(k,k:n)-L(k,1:k-1)*U(1:k-1,k:n);L(k+1:n,k)=(a(k+1:n,k)-L(k+1:n,1:k-1)*U(1:k-1,k))/U(k,k);end第二十八页，共73页。第二十九页，共73页。程序运行结果:第三十页，共73页。2．平方根法定理3.1

设A对称正定，则有非奇异下三角阵L，使分解方法：设第三十一页，共73页。第三十二页，共73页。第三十三页，共73页。第三十四页，共73页。3.列主元的三角分解法第三十五页，共73页。例4对下面的增广矩阵用列主元三角分解法分解。第三十六页，共73页。四、解三对角方程组——追赶法第三十七页，共73页。利用Gauss消元法得到同解的三角方程为第三十八页，共73页。第三十九页，共73页。第四十页，共73页。

,故有

(3.1)

第四十一页，共73页。解解

(3.1)—(3.3)叫追赶法，工作量小，非常有效。Ex:第三章4(3.2)（3.3）第四十二页，共73页。用Matlab实现追赶法求解线性方程组在Matlab程序编辑器中输入：functionx=pursue(a,b,flag)%a为系数矩阵；b为右端列向量；flag为0则显示过程ifnargin<3,flag=0;end%nargin为计算输入变量个数n=length(b);a=[a,b];fork=1:(n-1)a(k+1,k+1)=a(k+1,k+1)-a(k+1,k)/a(k,k)*a(k,k+1);a(k+1,n+1)=a(k+1,n+1)-a(k+1,k)/a(k,k)*a(k,n+1);a(k+1,k)=zeros(1,1);ifflag==0,a,endend第四十三页，共73页。%回代x(n)=a(n,n+1)/a(n,n);fork=n-1:-1:1x(k)=(a(k,n+1)-a(k,k+1)*x(k+1))/a(k,k)end第四十四页，共73页。第四十五页，共73页。程序运行结果:第四十六页，共73页。第四十七页，共73页。第四十八页，共73页。用Matlab实现改进平方根法求解线性方程组在Matlab程序编辑器中输入：functionx=squareroot(a,b,flag)ifnargin<3,flag=0;endn=length(b);a=[a,b];fork=1:(n-1)a(k+1,1)=a(k+1,1)/a(1,1);endafork=1:(n-2)forj=k+1:n+1a(k+1,j)=a(k+1,j)-a(1:k,j)*a(k+1,1:k);endforl=k+2:na(l,k+1)=(a(l,k+1)-a(1:k,k+1)*a(l,1:k))/a(k+1,k+1);

第四十九页，共73页。endaendfori=1:(n-1)a(n,n)=a(n,n)-a(i,n)*a(n,i);a(n,n+1)=a(n,n+1)-a(i,n+1)*a(n,i);endax=zeros(n,1);x(n)=a(n,n+1)/a(n,n);fork=n-1:-1:1x(k)=(a(k,n+1)-a(k,(k+1):n)*x((k+1):n))/a(k,k);end第五十页，共73页。第五十一页，共73页。程序运行结果：第五十二页，共73页。基本要求1.熟悉Gauss消去法及列主元Gauss消去法。2.熟悉矩阵的直接三角分解法的过程；2.熟悉平方根法的计算过程及其优缺点。作业：第三章9，11第五十三页，共73页。§3.2向量范数和矩阵范数第五十四页，共73页。第五十五页，共73页。第五十六页，共73页。第五十七页，共73页。第五十八页，共73页。在Matlab中可以用命令norm(A,p)来求范数。p为(2,1,inf)分别对应：2-范数，1-范数和∞-范数。举例如右图：第五十九页，共73页。一.简单迭代法1.迭代法建立.考虑

(矩阵B不唯一)对应写出

§3.3解线性方程组的迭代法产生向量序列若收敛,记第六十页，共73页。则于(3.4)两端取极限有：上式说明：

是解向量,从而当k充分大时注意:迭代阵B不唯一,影响收敛性。

解向量(3.4)叫简单迭代法,B叫迭代矩阵。第六十一页，共73页。第六十二页，共73页。二.Jacobi迭代法(k=0,1,2,…)(3.7)第六十三页，共73页。矩阵形式为：第六十四页，共73页。(3.8)

第六十五页，共73页。

称如下迭代法

(3.9)

为与(3.8)对应的迭代法，其迭代矩阵可用“代入法”求得。第六十六页，共73页。迭代法的收敛性第六十七页，共73页。例3.4分别用Jacobi迭代法和Seidel迭代法求解方程组取初始向量,要求

时迭代终止。解：因为系数矩阵严格对角占优，故Jacobi迭代法和Seidel方法对任意收敛。第六十八页，共73页。

首先，Jacobi迭代格式为第六十九页，共73页。计算结果可列表如下

注意：Seidel方法未必一定比Jacobi方法好。可参看作业中的题目。第七十页，共73页。

—松弛因子,=1即Seidel方法(3.10),