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

20.分别用Jacobi 迭代法、Gauss-seidel 迭代法解方程组???

??=++=++=-+.

122,1,122321

321

321x x x x x x x x x

解：Jacobi 迭代法收敛????

?

?????-=??????????133321x x x ，Gauss-seidel 迭代法不收敛。 27.编写LU 分解法、改进平方根法、追赶法的Matlab 程序，并进行相关数值实验。

3.将矩阵?????

????

???-=11

00

11021110

0201A 进行Doolittle 和Crout 分解 解：Doolittle 分解：结果如下，程序见后面。?

?

??????????-???????????

??-==2.100015001110020112.00001020010

0001

LU A Crout 分解：结果如下，程序见后面。?

?

???

???????-???????????

??-==1000

2.0100

1110

0201

2.110

05020010

0001LU A 7.用改进平方根法解方程组????

?

???????=??????????????????????

???----000142002511013101144321x x x x 解：结果如下，程序见后面。????

?????

??

?--=????????????0256.00513.00769.02821.04321x x x x 8（2）.用追赶法求解方程组?????

????

???????=?????????????????????????????????--------00001210001210001210001210001254321x x x x x

解：结果如下，程序见后面。????

?

????

???---=????????????1429.02413.02857.03571.04321x x x x

%检验输入参数，初始化

if nargin<2,error('more augments are needed');end

if nargin<3,x0=zeros(size(b));end

if nargin<4,delta=1e-13;end

if nargin<5,max1=100;end

if nargin>5,error('incorrect number of input');end

n=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))

error('A(i,i) equal to zero,divided by zero');

end

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 (err

flag=1;

break;

end

end

if flag==1

disp('The Jacobi method converges.')

x=x;

else

disp(['The Jacobi method does not converge with '...

,num2str(max1),' iterations'])

end

return

A=[1 2 -2;1 1 1;2 2 1] b=[1;1;1]

A =

1 2 -2

1 1 1

2 2 1

b =

1

1

1 jacobi(A,b)

The Jacobi method converges. ans =

-3

3

1

%检验输入参数，初始化

if nargin<2,error('more augments are needed');end if nargin<3,x0=zeros(size(b));end if nargin<4,delta=1e-13;end if nargin<5,max1=100;end

if nargin>5,error('incorrect number of input');end n=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))

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 (err

if flag==1

disp('The Gauss-seidel method converges.') x=x0; else

disp(['The Gauss-seidel method does not converge with '... ,num2str(max1),' iterations']) end return

gseid(A,b) The Gauss-seidel method does not converge with 100 iterations ans = 1.0e+31 * -9.1905 9.2222 -0.0634 A=[1 2 -2;1 1 1;2 2 1] b=[1;1;1] A = 1 2 -2 1 1 1 2 2 1 b = 1 1 1

%检验输入参数，初始化

b=size(A);n=b(1);

if b(1)~=b(2)

error('Matrix A should be a Square matrix.'); end

if 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;

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

end

return

A=[1 0 2 0;0 1 1 1;2 0 -1 1;0 0 1 1]

A =

1 0

2 0

0 1 1 1

2 0 -1 1

0 0 1 1 [L,U]=Crout(A)

L =

1.0000 0 0 0

0 1.0000 0 0

2.0000 0 -5.0000 0

0 0 1.0000 1.2000 U =

1.0000 0

2.0000 0

0 1.0000 1.0000 1.0000

0 0 1.0000 -0.2000

0 0 0 1.0000

%检验输入参数，初始化

b=size(A);n=b(1);

if b(1)~=b(2)

error('Matrix A should be a Square matrix.'); end

if 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进行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

end

return

A=[1 0 2 0;0 1 1 1;2 0 -1 1;0 0 1 1]

A =

1 0

2 0

0 1 1 1

2 0 -1 1

0 0 1 1 [L,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.0000 U =

1.0000 0

2.0000 0

0 1.0000 1.0000 1.0000

0 0 -5.0000 1.0000

0 0 0 1.2000

function [x]=improvecholesky(A,b,n)

%检验输入参数，初始化

L=zeros(n,n);D=diag(n,0);S=L*D;

for i=1:n

L(i,i)=1;

end

for 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

end

end

D(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)');

end

y=zeros(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); end

for i=n:-1:1

x(i)=y(i)-sum(L(i+1:n,i)'*x(i+1:n));

end

return

A=[4 1 -1 0;1 3 -1 0;-1 -1 5 2;0 0 2 4] b=[1;0;0;0]

A =

4 1 -1 0

1 3 -1 0

-1 -1 5 2

0 0 2 4

b =

1

0 n=4

[x]=improvecholesky(A,b,n)

n =

4

x =

0.2821

-0.0769

0.0513

-0.0256

function [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;

end

end

L=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);

end

if(u(n)~=0)

x(n)=y(n)/u(n);

end

for i=n-1:-1:1

x(i)=(y(i)-c(i)*x(i+1))/u(i);

end

a=[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 0 0 0 0

0.5000 1.0000 0 0 0

0 -0.4000 1.0000 0 0

0 0 -0.6250 1.0000 0

0 0 0 -0.7273 1.0000 U =

2.0000 -1.0000 0 0 0

0 2.5000 -1.0000 0 0

0 0 1.6000 -1.0000 0

0 0 0 1.3750 -1.0000

0 0 0 0 1.2727 x =

0.3571

-0.2857

-0.2143

-0.1429