第五章线性方程组求解的数值方法ppt

第五章线性方程组求解的数值方法ppt第1页，共42页，2023年，2月20日，星期三线性方程组求解的数值方法第2页，共42页，2023年，2月20日，星期三§5.1Gauss消去法与矩阵的LU分解基本思想：用逐次消去未知数的方法把原方程组化为三角形方程组再求解。消元：用初等变换将原方程组的系数矩阵化为三角形矩阵（简称三角阵）再求解的方法。回代：解出三角形方程组的最后一个方程，将求得的值逐步往前一个方程代入的方法。第3页，共42页，2023年，2月20日，星期三第4页，共42页，2023年，2月20日，星期三第5页，共42页，2023年，2月20日，星期三消元第6页，共42页，2023年，2月20日，星期三第7页，共42页，2023年，2月20日，星期三为什么选主元第8页，共42页，2023年，2月20日，星期三第9页，共42页，2023年，2月20日，星期三避免方法：高斯主元消元法第10页，共42页，2023年，2月20日，星期三functionX=uptrbk(A,b)%A是一个n阶矩阵。%b是一个n维向量。%X是线性方程组AX=b的解。[NN]=size(A);X=zeros(1,N+1);Aug=[Ab];forp=1:N-1[Y,j]=max(abs(Aug(p:N,p)));C=Aug(p,:);Aug(p,:)=Aug(j+p-1,:);Aug(j+p-1,:)=C;

高斯主元消去法的MATLAB实现

第11页，共42页，2023年，2月20日，星期三ifAug(p,p)==0'A是奇异阵，方程无惟一解'breakendfork=p+1:Nm=Aug(k,p)/Aug(p,p);Aug(k,p:N+1)=Aug(k,p:N+1)-m*Aug(p,p:N+1);endend%这里用到程序8-1中定义的函数backsub来进行回代。X=backsub(Aug(1:N,1:N)，Aug(1:N,N+1)；第12页，共42页，2023年，2月20日，星期三A=[0.00123;-13.7124.623;-21.0725.643];>>b=[123]';>>uptrbk(A,b)比较计算结果：主元消元法-0.4904-0.05100.3675

用x=A\b-0.4904-0.05100.3675

一般消元法-0.4-0.099800.4functionX=backsub(A,b)%A是一个n阶上三角非奇异阵。%b是一个n维向量。%X是线性方程组AX=b的解。n=length(b);X=zeros(n,1);X(n)=b(n)/A(n,n);fork=n-1:-1:1X(k)=(b(k)-A(k,k+1:n)*X(k+1:n))/A(k,k);End用回代法求解上三角线性方程组AX=B，其中A为非奇异。第13页，共42页，2023年，2月20日，星期三矩阵的LU分解第14页，共42页，2023年，2月20日，星期三第15页，共42页，2023年，2月20日，星期三第16页，共42页，2023年，2月20日，星期三第17页，共42页，2023年，2月20日，星期三第18页，共42页，2023年，2月20日，星期三functionX=lufact(A,b)%A为n阶矩阵。%b是n维向量。%X是所求的AX=b的解。[N,N]=size(A);X=zeros(N,1);Y=zeros(N,1);C=zeros(1,N);R=1:N;forp=1:N-1[max1,j]=max(abs(A(p:N,p)));C=A(p,:);A(p,:)=A(j+p-1,:);A(j+p-1,:)=C;d=R(p);R(p)=R(j+p-1);R(j+p-1)=d;ifA(p,p)==0'A是奇异阵，方程组无惟一解'

第19页，共42页，2023年，2月20日，星期三breakendfork=p+1:Nmult=A(k,p)/A(p,p);

A(k,p)=mult;A(k,p+1:N)=A(k,p+1:N)-mult*A(p,p+1:N);endendY(1)=b(R(1));fork=2:NY(k)=b(R(k))-A(k,1:k-1)*Y(1:k-1);endX(N)=Y(N)/A(N,N);fork=N-1:-1:1

X(k)=(Y(k)-A(k,k+1:N)*X(k+1:N))/A(k,k);end第20页，共42页，2023年，2月20日，星期三§5.2Cholesky分解（平方根法)第21页，共42页，2023年，2月20日，星期三第22页，共42页，2023年，2月20日，星期三第23页，共42页，2023年，2月20日，星期三§5.3向量范数与矩阵范数第24页，共42页，2023年，2月20日，星期三第25页，共42页，2023年，2月20日，星期三第26页，共42页，2023年，2月20日，星期三向量的收敛第27页，共42页，2023年，2月20日，星期三第28页，共42页，2023年，2月20日，星期三矩阵范数第29页，共42页，2023年，2月20日，星期三第30页，共42页，2023年，2月20日，星期三常用的矩阵范数第31页，共42页，2023年，2月20日，星期三矩阵的收敛第32页，共42页，2023年，2月20日，星期三谱半径第33页，共42页，2023年，2月20日，星期三§5.4三种迭代法第34页，共42页，2023年，2月20日，星期三雅可比迭代法的Matlab程序functionX=jacobi(A,b,P,delta,max1)%A是n维非奇异阵。%b是n维向量。%P是初值。%delta是误差界。%max1是给定的迭代最高次数。%X为所求的方程组AX=b的近似解。N=length(b);fork=1:max1forj=1:NX(j)=(b(j)-A(j,[1:j-1,j+1:N])*P([1:j-1,j+1:N]))/A(j,j);enderr=abs(norm(X'-P));P=X'；

if(err<delta)breakendendX=X';k,err;第35页，共42页，2023年，2月20日，星期三第36页，共42页，2023年，2月20日，星期三functionX=gseid(A,b,P,delta,max1)%A是n维非奇异阵。%b是n维向量。%P是初值。%delta是误差界。%max1是给定的迭代最高次数。%X为所求的方程组AX=b的近似解。N=length(b);fork=1:max1forj=1:Nifj==1X(1)=(b(1)-A(1,2:N)*P(2:N))/A(1,1);

elseifj==NX(N)=(b(N)-A(N,1:N-1)*(X(1:N-1))')/A(N,N);elseX(j)=(b(j)-A(j,1:j-1)*X(1:j-1)-A(j,j+1:N)*P(j+1:N))/A(j,j);endenderr=abs(norm(X'-P));P=X'；if(err<delta)breakendendX=X';err,k塞德