Jacobi반복. Gauss‐Seidel 반복

카테고리 없음

Jacobi반복. Gauss‐Seidel 반복

mathboy 2011. 12. 29. 16:27

function [x, iflag, itnum] = Jacobi(A,b,x0,delta,max_it)

% function [x, iflag, itnum] = Jacobi(A,b,x0,delta,max_it)

% A program implementing the Jacobi iteration method to solve

% the linear system Ax=b.

% Input

% A: square coefficient matrix

% b: right side vector

% x0: initial guess

% delta: error tolerance for the relative difference between

% two consecutive iterates

% max_it: maximum number of iterations to be allowed

% Output

% x: numerical solution vector

% iflag: 1 if a numerical solution satisfying the error

% tolerance is found within max_it iterations

% -1 if the program fails to produce a numerical

% solution in max_it iterations

% itnum: the number of iterations used to compute x

% initialization

n = length(b);

iflag = 1;

k = 0;

% create a vector with diagonal elements of A

diagA = diag(A);

% modify A to make its diagonal elements zero

A = A-diag(diag(A));

% iteration

while k < max_it

k = k+1;

x = (b-A*x0)./diagA

relerr = norm(x-x0,inf)/(norm(x,inf)+eps);

x0 = x;

if relerr < delta

break

end

itnum = k;

if (itnum == max_it)

iflag = -1

end

function [x, iflag, itnum] = GS(A,b,x0,delta,max_it)

% function [x, iflag, itnum] = GS(A,b,x0,delta,max_it)

% A program implementing the Gauss-Seidel iteration method to solve

% the linear system Ax=b.

% Input

% A: square coefficient matrix

% b: right side vector

% x0: initial guess

% delta: error tolerance for the relative difference between

% two consecutive iterates

% max_it: maximum number of iterations to be allowed

% Output

% x: numerical solution vector

% iflag: 1 if a numerical solution satisfying the error

% tolerance is found within max_it iterations

% -1 if the program fails to produce a numerical

% solution in max_it iterations

% itnum: the number of iterations used to compute x

% initialization

n = length(b);

iflag = 1;

k = 0;

x = x0;

% iteration

while k < max_it

k = k+1;

x(1) = (b(1)-A(1,2:n)*x0(2:n))/A(1,1);

for i = 2:n

if i < n

x(i) = (b(i)-A(i,1:i-1)*x(1:i-1) ...

-A(i,i+1:n)*x0(i+1:n))/A(i,i);

else

x(n) = (b(n)-A(n,1:n-1)*x(1:n-1))/A(n,n);

end

relerr = norm(x-x0,inf)/(norm(x,inf)+eps);

x0 = x

if relerr < delta

break

end

itnum = k;

if (itnum == max_it)

iflag = -1

end

>> A=[5 7 6 5;7 10 8 7;6 8 10 9;5 7 9 10];

>> b=[1; -1; -1; 1];

>> x0=[100; -100; -50; 50];

>> x0=[0;0;0;0];

>> Jacobi(A,b,x0,0.00005,10)

～～～～～～～～～～～～～～～～～

iflag =

-1

ans =

1.0e+194 *

-2.750072354811125

-1.913739097811525

-1.939784316974014

-1.801629708784825

>> GS(A,b,x0,0.00005,100000)

～～～～～～～～～～～～～～～～～

x0 =

1.0e+002 *

1.338540688530760

-0.807102131160679

-0.344551426577325

0.206797431466687