第05讲 线性代数方程组求解的迭代法

1、17:11:17解线性代数方程组解线性代数方程组的迭代法的迭代法17:11:17迭代法的一般形式迭代法的一般形式17:11:17Jacobi迭代法迭代法算法算法17:11:18Jacobi迭代法迭代法程序程序1A=10 -1 -2; -1 10 -2; -1 -1 5;b=72; 83; 42;n=3;x0=0; 0; 0;x=0; 0; 0;epsilon = 1.0e-6;N=10;k=1;17:11:18Jacobi迭代法迭代法程序程序2 for i=1:n x(i)=b(i); for j=1:n if j=i x(i) = x(i) - A(i,j)*x0(j); end end x

2、(i)=x(i)/A(i,i); end17:11:18Jacobi迭代法迭代法程序程序3while k=N for i=1:n end if norm(x-x0)epsilon return end k=k+1; x0=x;end17:11:18Gauss-Seidel迭代法迭代法算法算法17:11:18Gauss-Seidel迭代法迭代法程序程序clcclearA=10 -1 -2; -1 10 -2; -1 -1 5;b=72; 83; 42;n=3;x0=0; 0; 0;x=0; 0; 0;epsilon = 1.0e-6;N=10;k=1;while k=N if norm(x-x0

3、)epsilon return end k=k+1; x0=x;endx(i)=b(i);for j=1:(i-1) x(i) = x(i) - A(i,j)*x(j);endfor j=(i+1):n x(i) = x(i) - A(i,j)*x0(j);endx(i)=x(i)/A(i,i);x(n)=b(n);for j=1:(n-1) x(n)=x(n)-A(n,j)*x(j);endx(n)=x(n)/A(n,n);x(1)=b(1);for j=2:n x(1)=x(1)-A(1,j)*x0(j);endx(1)=x(1)/A(1,1);(0)111112()/njjjxba xa

4、1(0)11()/(2,1)iniiijjijjiijj ixba xa xain 11()/nnnnjjnnjxba xa17:11:18Gauss-Seidel迭代法迭代法程序改进程序改进while k=N MaxError = 0; for i=1:n xi=b(i); for j=1:n if i=j xi = xi - A(i,j)*x(j); end end xi=xi/A(i,i); MaxError=max(MaxError abs(xi-x(i); x(i)=xi; end x if MaxErrorepsilon return end k=k+1; x0=x;end(0)j

5、jxjixjiixix17:11:18松弛法松弛法算法算法17:11:19松弛法松弛法程序程序while k=N MaxError = 0; for i=1:n xi=b(i); for j=1:n if i=j xi = xi - A(i,j)*x(j); end end xi=xi/A(i,i); xi = x(i) + omega*(xi-x(i); MaxError=max(MaxError abs(xi-x(i); x(i)=xi; end x if MaxErrorepsilon return end k=k+1; x0=x;endclcclearA=2 -1 0; -1 2 -1

6、; 0 -1 2;b=1; 0; 1.8;n=3;x=1; 1; 1;omega=1.4;epsilon = 1.0e-6;N=10;k=1;17:11:19线性方程组求解与极小化问题线性方程组求解与极小化问题17:11:19最速下降法最速下降法算法算法17:11:19最速下降法最速下降法程序程序clcclearA=4 3 0; 3 3 -1; 0 -1 4;b=24; 30; -24;n=3;x=0; 0; 0;epsilon = 1.0e-6;N=100;k=1;while k=N r=b-A*x; lambda=r*r/(r*A*r); x=x+lambda*r; x if lambda

7、*norm(r)epsilon return end k=k+1;end17:11:19共轭梯度法共轭梯度法原理原理17:11:20共轭梯度法共轭梯度法原理原理17:11:20向量的正交化向量的正交化17:11:20共轭向量组共轭向量组17:11:20共轭梯度法算例共轭梯度法算例 17:11:20共轭梯度法的共轭梯度法的Matlab程序程序 clear allclc% 定义矩阵和向量A=4 3 0; 3 4 -1; 0 -1 4;b=24; 30; -24;%设置初值x0 = 0; 0; 0;%计算在处的负梯度r0 = b - A*x0d0 = r0;%计算新的极小点lambda0 = r0*

8、d0/(d0*A*d0)x1 = x0 + lambda0*d0%计算处的负梯度r1 = b - A*x1% 将正交化beta0 = -(r1*A*d0)/(d0*A*d0)d1 = r1 + beta0*d0%计算新的极小点lambda1 = r1*d1/(d1*A*d1)x2 = x1 + lambda1*d1%计算处的负梯度r2 = b - A*x2% 将正交化beta1 = -(r2*A*d1)/(d1*A*d1)d2 = r2 + beta1*d1%计算新的极小点lambda2 = r2*d2/(d2*A*d2)x3 = x2 + lambda2*d217:11:21运行结果运行结果

9、x1 = 3.5258 4.4072 -3.5258x2 = 2.8580 4.1490 -4.9542x3 = 3.0000 4.0000 -5.000017:11:21共轭梯度法的特性共轭梯度法的特性l理论上可以证明，最多经过n次迭代就能得到准确解。l考虑实际计算时的舍入误差，经过n次迭代一般得不到准确解。l可以将经过n次迭代后得到的解作为初值，再进行一轮n次迭代。17:11:21改进后的程序改进后的程序clear allclcA=4 3 0; 3 4 -1; 0 -1 4;b=24; 30; -24;epsilon = 1.0e-6;n = length(b);x = 0; 0; 0;for loop=1:2*n r = b - A*x; if norm(r) epsilon; break; end if mod(l