科学工程计算与matlab编程4.pptVIP

function y=jacobi(a,b,x0) D=diag(diag(a)); U=-triu(a,1); L=-tril(a,-1); B=D\(L+U); f=D\b; y=B*x0+f; n=1; while norm(y-x0)=1.0e-6 x0=y; y=B*x0+f; n=n+1; end n 例：用Jacobi方法求解， 设x(0)=0,精度为10-6。 a=[10 -1 0; -1 10 -2; 0 -2 10]; b=[9; 7; 6]; jacobi(a,b,[0;0;0]) n = 11 ans = 0.9958 0.9579 0.7916 4.3.2 Gauss-Seidel迭代法 由原方程构造迭代方程 x(k+1)=G x(k)+f 其中：G=(D-L)-1 U, f=(D-L)-1 b D=diag[a11,a22,…,ann], -L、-U分别为A的严格下、上三角部分（不包括对角线元素）. function y=seidel(a,b,x0) D=diag(diag(a));U=-triu(a,1);L=-tril(a,-1); G=(D-L)\U ;f=(D-L)\b; y=G*x0+f; n=1; while norm(y-x0)=1.0e-6 x0=y; y=G*x0+f; n=n+1; end n 例：对上例用Gauss-Seidel迭代法求解 a=[10 -1 0; -1 10 -2; 0 -2 10]; b=[9; 7; 6]; seidel(a,b,[0;0;0]) n = 7 ans = 0.9958 0.9579 0.7916 例：分别用Jacobi和G-S 法迭代求解,看是否收敛。 a=[1 2 -2; 1 1 1; 2 2 1]; b=[9; 7; 6]; jacobi(a,b,[0;0;0]) n = 4 ans = -27 26 8 seidel(a,b,[0;0;0]) n = 1011 ans = 1.0e+305 * -Inf Inf -1.7556 4.3.3 SOR迭代法 在很多情况下，J法和G-S法收敛较慢，所以考虑对G-S法进行改进。于是引入一种新的迭代法－逐次超松弛迭代法（Succesise Over-Relaxation),记为SOR法。 迭代公式为： X(k+1)= (D-wL)-1((1-w)D+wU)x(k) + w(D-wL)-1 b 其中：w最佳值在[1, 2）之间，不易计算得到，因此 w通常有经验给出。 function y=sor(a,b,w,x0) D=diag(diag(a));U=-triu(a,1);L=-tril(a,-1); M=(D-w*L)\((1-w)*D+w*U); f=(D-w*L)\b*w; y=M*x0+f; n=1; while norm(y-x0)=1.0e-6 x0=y; y=M*x0+f; n=n+1; end n 例：上例中，当w=1.103时，用SOR法求解原方程。 a=[10 -1 0; -1 10 -2; 0 -2 10]; b=[9; 7; 6]; sor(a,b,1.103,[0;0;0]) n = 8 ans = 0.9958 0.9579 0.7916 4.3.4 两步迭代法 当线性方程系数矩阵为对称正定时，可用一种特殊的迭代法来解决，其迭代公式为： （D-L）x(k+1/2) =U x(k) +b （D-U）x(k+1)=Lx(k+1/2) +b = x(k+1/2) =(D-L)-1 U x(k) + (D-L)-1 b x(k+1)= (D-U)-1 Lx(k+1/2) + (D-U)-1 b function y=twostp(a,b,x0) D=diag(diag(a));U=-triu(a,1);L=-tril(a,-1); G1==(D-L)\U; f1=(D-L)\b; G2==(D-U)\L; f1=(D-U)\b; y=G1*x0+f1; y=G2*y+f2;