《第五章线性方程组求解的数值方法ppt精选文档.ppt》由会员分享,可在线阅读,更多相关《第五章线性方程组求解的数值方法ppt精选文档.ppt(42页珍藏版)》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、第五章线性方程组求解的数值方法ppt本讲稿第一页,共四十二页线性方程组求解的数值方法线性方程组求解的数值方法本讲稿第二页,共四十二页5.1 Gauss5.1 Gauss消去法与矩阵的消去法与矩阵的LULU分解分解基本思想:基本思想:用逐次消去未知数的方法把原方程组化为三角形方程组再求解。消元:消元:用初等变换将原方程组的系数矩阵化为三角形矩阵(简称三角阵)再求解的方法。回代:回代:解出三角形方程组的最后一个方程,将求得的值逐步往前一个方程代入的方法。本讲稿第三页,共四十二页本讲稿第四页,共四十二页本讲稿第五页,共四十二页消元本讲稿第六页,共四十二页本讲稿第七页,共四十二页为什么选主元为什么选主
2、元本讲稿第八页,共四十二页本讲稿第九页,共四十二页避免方法:高斯主元消元法本讲稿第十页,共四十二页function X=uptrbk(A,b)%A是一个n阶矩阵。%b是一个n维向量。%X是线性方程组AX=b的解。N N=size(A);X=zeros(1,N+1);Aug=A b;for p=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实现本讲稿第十一页,共四十二页if Aug(p,p)=0 A是奇异阵,方程无惟一解 break end for k=p+1:
3、N m=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);本讲稿第十二页,共四十二页 A=0.001 2 3;-1 3.712 4.623;-2 1.072 5.643;b=1 2 3;uptrbk(A,b)比较计算结果比较计算结果:主元消元法 -0.4904 -0.0510 0.3675 用x=Ab -0.4904 -0.0510 0.3675 一般消元法 -0.4 -0.09980
4、0.4function X=backsub(A,b)%A是一个n阶上三角非奇异阵。%b是一个n维向量。%X是线性方程组AX=b的解。n=length(b);X=zeros(n,1);X(n)=b(n)/A(n,n);for k=n-1:-1:1X(k)=(b(k)-A(k,k+1:n)*X(k+1:n)/A(k,k);End用回代法求解上三角线性方程组AX=B,其中A为非奇异。本讲稿第十三页,共四十二页矩阵的LU分解本讲稿第十四页,共四十二页本讲稿第十五页,共四十二页本讲稿第十六页,共四十二页本讲稿第十七页,共四十二页本讲稿第十八页,共四十二页function X=lufact(A,b)%A为
5、为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;for p=1:N-1max1,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;if A(p,p)=0 A是奇异阵,方程组无惟一解 本讲稿第十九页,共四十二页 breakendfor k=p+1:N mult=A(k,p)/A(p,p);A(k,p)=mult;A(
6、k,p+1:N)=A(k,p+1:N)-mult*A(p,p+1:N);endendY(1)=b(R(1);for k=2:N Y(k)=b(R(k)-A(k,1:k-1)*Y(1:k-1);endX(N)=Y(N)/A(N,N);for k=N-1:-1:1 X(k)=(Y(k)-A(k,k+1:N)*X(k+1:N)/A(k,k);end本讲稿第二十页,共四十二页5.2 Cholesky5.2 Cholesky分解(平方根法分解(平方根法)本讲稿第二十一页,共四十二页本讲稿第二十二页,共四十二页本讲稿第二十三页,共四十二页5.3 5.3 向量范数与矩阵范数向量范数与矩阵范数本讲稿第二十四页
7、,共四十二页本讲稿第二十五页,共四十二页本讲稿第二十六页,共四十二页向量的收敛本讲稿第二十七页,共四十二页本讲稿第二十八页,共四十二页矩阵范数本讲稿第二十九页,共四十二页本讲稿第三十页,共四十二页常用的矩阵范数本讲稿第三十一页,共四十二页矩阵的收敛本讲稿第三十二页,共四十二页谱半径本讲稿第三十三页,共四十二页5.4 5.4 三种迭代法三种迭代法本讲稿第三十四页,共四十二页雅可比迭代法的Matlab程序function X=jacobi(A,b,P,delta,max1)%A是n维非奇异阵。%b是n维向量。%P是初值。%delta是误差界。%max1是给定的迭代最高次数。%X为所求的方程组AX=
8、b的近似解。N=length(b);for k=1:max1for j=1:N X(j)=(b(j)-A(j,1:j-1,j+1:N)*P(1:j-1,j+1:N)/A(j,j);end err=abs(norm(X-P);P=X;if(errdelta)break endendX=X;k,err;本讲稿第三十五页,共四十二页本讲稿第三十六页,共四十二页function X=gseid(A,b,P,delta,max1)%A是n维非奇异阵。%b是n维向量。%P是初值。%delta是误差界。%max1是给定的迭代最高次数。%X为所求的方程组AX=b的近似解。N=length(b);for k=1
9、:max1for j=1:N if j=1 X(1)=(b(1)-A(1,2:N)*P(2:N)/A(1,1);elseif j=N X(N)=(b(N)-A(N,1:N-1)*(X(1:N-1)/A(N,N);else X(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(errdelta)breakendendX=X;err,k塞德尔迭代法的迭代程序:本讲稿第三十七页,共四十二页本讲稿第三十八页,共四十二页本讲稿第三十九页,共四十二页5.5 5.5 迭代法的收敛性迭代法的收敛性本讲稿第四十页,共四十二页本讲稿第四十一页,共四十二页本讲稿第四十二页,共四十二页