Question on implementing Jacobi Method

function [x,info,relres] = myJacobi(A,b,MaxIter,TOL)

%Jacobi Method

A=zeros(n,n);

%x0=zeros(n,1);

x=zeros(n,1);

b=zeros(n,1);

MaxIter=10*n;

TOL=1.e-4;

k=1;

while k<=MaxIter

for i=1:n

x(i)=(b(i)-a(i,[1:i-1,i+1:n])*x([1:i-1,i+1:n]))/a(i,i);

if norm(A*x-b)/norm(b)<=TOL %abs(x-x0)<TOL

print(x);

end

end

k=k+1;

for i=1:n

%x0(i)=x(i);

print(x(i));

end

end

end

%This is demo codes I am trying to call the Jacobi function that I defined at the very beginning to solve for another system.

N=50;

A=gallery('poisson',N);

n=size(A,1);

xs=ones(n,1);

b=A*xs(:);

MaxIter=10*n;

TOL=1.e-4;

t=cputime;

[x,info,relres]=myJacobi(A,b,MaxIter,TOL); % Jacobi example

mycput = cputime - t;

relerr = norm(xs-x)/norm(xs);

fprintf('size(A,1)=%d,info=%d, relres=%8.5e, relerr=%8.5e\n',...

n, info, relres, relerr);

fprintf('time =%8.5e\n', mycput);

The above is my code. What I am doing is divided into two parts. In the first part, I am inplementing Jacobi method into Matlab. In the second part which is outside of the function that I defined, I am trying to calculate another system by calling the function myJacobi. However, I do not know why matlab keeps saying the line with N=50 is not right. Can somebody help me figure it out? What are things I might need to change? Thanks a lot.

Answer 1

[Alan Stevens]

More like this perhaps

N=10;  % 50;  
A=gallery('poisson',N);
n=size(A,1);
xs=ones(n,1);
b=A*xs(:);
MaxIter=10*n;
TOL=1.e-4;
t=cputime;
[x]=myJacobi(A,b,MaxIter,TOL,n); % Jacobi example
mycput = cputime - t;
relerr = norm(xs-x)/norm(xs);
fprintf('size(A,1) = %d, relerr = %8.5e\n',n, relerr);
fprintf('time = %8.5e\n', mycput);

function [x] = myJacobi(A,b,MaxIter,TOL,n)
%Jacobi Method
x=zeros(n,1);
k=1;
err = 1;
while k<=MaxIter && err>TOL
    for i=1:n
        x(i)=(b(i)-A(i,[1:i-1,i+1:n])*x([1:i-1,i+1:n]))/A(i,i);
    end
    err = norm(A*x-b)/norm(b); 
    k=k+1;
end
disp(k)
end

Note: you hadn't implemented any calculation for info and relres, so they've been removed in the version here.