常州大学数值分析作业—第二章

1、20. 分别用Jacobi迭代法、Gauss-seidel迭代法解方程组解：Jacobi迭代法收敛，Gauss-seidel迭代法不收敛。27. 编写LU分解法、改进平方根法、追赶法的Matlab程序，并进行相关数值实验。3.将矩阵进行Doolittle和Crout分解解：Doolittle分解：结果如下，程序见后面。 Crout分解：结果如下，程序见后面。 7.用改进平方根法解方程组解：结果如下，程序见后面。 8（2）.用追赶法求解方程组解：结果如下，程序见后面。function x,iternum,flag=jacobi(A,b,x0,delta,max1)%检验输入参数，初始化if na

2、rgin2,error(more augments are needed);endif nargin3,x0=zeros(size(b);endif nargin4,delta=1e-13;endif nargin5,error(incorrect number of input);endn=length(b);x=0*b;flag=0;iternum=0;%用Jacobi迭代法解方程组for k=1:max1 iternum=iternum+1; for i=1:n if abs(A(i,i)eps error(A(i,i) equal to zero,divided by zero); e

3、nd x(i)=(b(i)-A(i,1:i-1,i+1:n)*x0(1:i-1,i+1:n)/A(i,i); end err=norm(x-x0); relerr=err/(norm(x)+eps); x0=x; if (errdelta)|(relerrdelta) flag=1; break; endendif flag=1 disp(The Jacobi method converges.) x=x;else disp(The Jacobi method does not converge with . ,num2str(max1), iterations)endreturnA=1 2

4、-2;1 1 1;2 2 1b=1;1;1A = 1 2 -2 1 1 1 2 2 1b = 1 1 1jacobi(A,b)The Jacobi method converges.ans = -3 3 1function x,iternum,flag=gseid(A,b,x0,delta,max1)%检验输入参数，初始化if nargin2,error(more augments are needed);endif nargin3,x0=zeros(size(b);endif nargin4,delta=1e-13;endif nargin5,error(incorrect number o

5、f input);endn=length(b);flag=0;iternum=0;%用Gauss-Seidel迭代法解方程组for k=1:max1 iternum=iternum+1; x=x0; for i=1:n if abs(A(i,i)eps error(A(i,i) equal to zero,divided by zero); end x0(i)=(b(i)-A(i,1:i-1,i+1:n)*x0(1:i-1,i+1:n)/A(i,i); end err=norm(x-x0); relerr=err/(norm(x0)+eps); if (errdelta)|(relerrdel

6、ta) flag=1; break; endendif flag=1 disp(The Gauss-seidel method converges.) x=x0;else disp(The Gauss-seidel method does not converge with . ,num2str(max1), iterations)endreturngseid(A,b)The Gauss-seidel method does not converge with 100 iterationsans = 1.0e+31 * -9.1905 9.2222 -0.0634A=1 2 -2;1 1 1;

7、2 2 1b=1;1;1A = 1 2 -2 1 1 1 2 2 1b = 1 1 1function L,U=Crout(A)%检验输入参数，初始化b=size(A);n=b(1);if b(1)=b(2) error(Matrix A should be a Square matrix.);endif n=rank(A) error(Matrix A should be FULL RANK.);end L=zeros(n,n);U=zeros(n,n); for i=1:n U(i,i)=1;end%将矩阵A进行Crout分解for k=1:n for i=k:n temp_sum=0;

8、for t=1:k-1 temp_sum=temp_sum+L(i,t)*U(t,k); end L(i,k)=A(i,k)-temp_sum; end for j=k+1:n temp_sum=0; for t=1:k-1 temp_sum=temp_sum+L(k,t)*U(t,j); end U(k,j)=(A(k,j)-temp_sum)/L(k,k); end endreturnA=1 0 2 0;0 1 1 1;2 0 -1 1;0 0 1 1A = 1 0 2 0 0 1 1 1 2 0 -1 1 0 0 1 1L,U=Crout(A)L = 1.0000 0 0 0 0 1.0

9、000 0 0 2.0000 0 -5.0000 0 0 0 1.0000 1.2000U = 1.0000 0 2.0000 0 0 1.0000 1.0000 1.0000 0 0 1.0000 -0.2000 0 0 0 1.0000function L,U=Crout(A)%检验输入参数，初始化b=size(A);n=b(1);if b(1)=b(2) error(Matrix A should be a Square matrix.);endif n=rank(A) error(Matrix A should be FULL RANK.);end L=zeros(n,n);U=zer

10、os(n,n); for i=1:n U(i,i)=1;end%将矩阵A进行Doolittle分解for k=1:n for i=k:n temp_sum=0; for t=1:k-1 temp_sum=temp_sum+L(i,t)*U(t,k); end L(i,k)=A(i,k)-temp_sum; end for j=k+1:n temp_sum=0; for t=1:k-1 temp_sum=temp_sum+L(k,t)*U(t,j); end U(k,j)=(A(k,j)-temp_sum)/L(k,k); end endreturnA=1 0 2 0;0 1 1 1;2 0 -

11、1 1;0 0 1 1A = 1 0 2 0 0 1 1 1 2 0 -1 1 0 0 1 1L,U=Doolittle(A)L = 1.0000 0 0 0 0 1.0000 0 0 2.0000 0 1.0000 0 0 0 -0.2000 1.0000U = 1.0000 0 2.0000 0 0 1.0000 1.0000 1.0000 0 0 -5.0000 1.0000 0 0 0 1.2000function x=improvecholesky(A,b,n)%检验输入参数，初始化L=zeros(n,n);D=diag(n,0);S=L*D;for i=1:n L(i,i)=1;e

12、ndfor i=1:n for j=1:n if (eig(A)=0)|(A(i,j)=A(j,i) disp(Matrix A should be a Positive-definite matrix.); break; end endendD(1,1)=A(1,1);%用改进平方根法解方程组for i=2:n for j=1:i-1 S(i,j)=A(i,j)-sum(S(i,1:j-1)*L(j,1:j-1); L(i,1:i-1)=S(i,1:i-1)/D(1:i-1,1:i-1); end D(i,i)=A(i,i)-sum(S(i,1:i-1)*L(i,1:i-1);endy=ze

13、ros(n,1);x=zeros(n,1);for i=1:n y(i)=(b(i)-sum(L(i,1:i-1)*D(1:i-1,1:i-1)*y(1:i-1)/D(i,i);endfor i=n:-1:1 x(i)=y(i)-sum(L(i+1:n,i)*x(i+1:n);endreturnA=4 1 -1 0;1 3 -1 0;-1 -1 5 2;0 0 2 4b=1;0;0;0A = 4 1 -1 0 1 3 -1 0 -1 -1 5 2 0 0 2 4b = 1 0 0 0n=4x=improvecholesky(A,b,n)n = 4x = 0.2821 -0.0769 0.051

14、3 -0.0256function L,U,x=pursue(a,b,c,f)n=length(b);u(1)=b(1);for i=2:n if(u(i-1)=0) l(i-1)=a(i-1)/u(i-1); u(i)=b(i)-l(i-1)*c(i-1); else break; endendL=eye(n)+diag(l,-1);U=diag(u)+diag(c,1);x=zeros(n,1);y=x;y(1)=f(1);for i=2:n y(i)=f(i)-l(i-1)*y(i-1);endif(u(n)=0) x(n)=y(n)/u(n);endfor i=n-1:-1:1 x(i)=(y(i)-c(i)*x(i+1)/u(i);enda=1 -1 -1 -1;b=2 2 2 2 2;c=-1 -1 -1 -1;f=1 0 0 0 0 ;L,U,x=pursue(a,b,c,f)L = 1.0000