Question on implementing Jacobi Method

5 次查看(过去 30 天)
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.

回答(1 个)

Alan Stevens
Alan Stevens 2021-3-16
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.

类别

Help CenterFile Exchange 中查找有关 Programming 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by