iteration method numerical analysis - matrix

hello every one, i am trying to form an iiteration loop to converge the best solution

my code is here, could you please tell me how to put this in iteration loop and can control the iteration number according to various conditions

%%% Point Gauss Seidel (PGS)  

close all, clc
format compact
%Defining the constants
L=1;                                % length
H=2;                                % Height
deltax=0.05
deltay=0.05
Beta=deltax/deltay
Beta2=Beta^2
jn=H/deltay                         % Maximum number of grid points along y
im=L/deltax                         % Maximum number of grid points along x
T1=100; T2=0; T3=0; T4=0;           % boundary conditions
Errormax=0.01;
y=2:-deltay:0;
x=0:deltax:L;

Told=zeros(jn+1,im+1);

% set boundary conditions
Told(1,1:im+1)=T1;
Told(jn+1,1:im+1)=T3;
Told(2:jn+1,1)=T2;
Told(2:jn+1,im+1)=T4;

% Iteration 1
for i=2:jn 
for j=2:im
Told(i,j)=(1/(2*(1+Beta2)))* (Told(i-1,j)+Told(i+1,j)+Beta2*(Told(i,j+1)+Told(i,j-1)) ) ;       
end
end
% Iteration 2

for i=2:jn 
for j=2:im
Told(i,j)=(1/(2*(1+Beta2)))* (Told(i-1,j)+Told(i+1,j)+Beta2*(Told(i,j+1)+Told(i,j-1)) ) ;       
end
end
% Iteration 3

for i=2:jn 
for j=2:im
Told(i,j)=(1/(2*(1+Beta2)))* (Told(i-1,j)+Told(i+1,j)+Beta2*(Told(i,j+1)+Told(i,j-1)) ) ;       
end
end

Told=flip(Told);
disp(' ');disp(' ');
disp([y' Told(:,find(abs(x-0.0) < 0.001)) Told(:,find(abs(x-0.2) < 0.001))...
    Told(:,find(abs(x-0.4) < 0.001)) Told(:,find(abs(x-0.6) < 0.001))...
    Told(:,find(abs(x-0.8) < 0.001)) Told(:,find(abs(x-1.0) < 0.001))])

[mehmet salihi]

i could do it like this, but i entered manually the iteration number to be 410, is there any way to know the iteraion number.

for k=1:410
for i=2:jn 
    Told(i,j)=(1/(2*(1+Beta2)))* (Told(i-1,j)+Told(i+1,j)+Beta2*(Told(i,j+1)+Told(i,j-1)));
for j=2:im
    Told(i,j)=(1/(2*(1+Beta2)))* (Told(i-1,j)+Told(i+1,j)+Beta2*(Told(i,j+1)+Told(i,j-1))) ;       
end
end
end

[Torsten]

% Initialize T_old to the initial guess values

T_old = ...;

T_new = T_old;

error = 1.0;

eps = 1e-4;

while error > eps

% calculate T_new in the double loop (write T_new everywhere in your code where you now use T_old)

error = max(abs(T_old-T_new),[],'all')

T_old = T_new;

end

[mehmet salihi]

the errormax is given in the question,

i did the code as you mention, could you please chech if it is correct , appretiate it.

in addition, i want to make the output variables of this code to be round from 4 decimal to 3. could you add the code for this

%%% Point Gauss Seidel (PGS)  

close all, clc
format short 
%Defining the constants
L=1;                                % length
H=2;                                % Height
deltax=0.05
deltay=0.05
Beta=deltax/deltay
Beta2=Beta^2
jn=H/deltay                         % Maximum number of grid points along y
im=L/deltax                         % Maximum number of grid points along x
T1=100; T2=0; T3=0; T4=0;           % boundary conditions
y=2:-deltay:0;
x=0:deltax:L;

T_old=zeros(jn+1,im+1);

% set boundary conditions
T_old(1,1:im+1)=T1;
T_old(jn+1,1:im+1)=T3;
T_old(2:jn+1,1)=T2;
T_old(2:jn+1,im+1)=T4;

% Initialize T_old to the initial guess values
for i=2:jn 
for j=2:im
T_old(i,j)=(1/(2*(1+Beta2)))* (T_old(i-1,j)+T_old(i+1,j)+Beta2*(T_old(i,j+1)+T_old(i,j-1)) ) ;       
end
end
T_new = T_old;

eps = 1e-4;
Errormax=0.01;
t=0;
while Errormax > eps
    t=t+1
    for i=2:jn 
        for j=2:im
            T_new(i,j)=(1/(2*(1+Beta2)))* (T_old(i-1,j)+T_old(i+1,j)+Beta2*(T_old(i,j+1)+T_old(i,j-1)) ) ;       
        end
    end

Errormax = max(abs(T_old-T_new),[],'all');
T_old = T_new;
end



disp(' ');disp(' ');
disp([y' T_old(:,find(abs(x-0.0) < 0.001)) T_old(:,find(abs(x-0.2) < 0.001))...
    T_old(:,find(abs(x-0.4) < 0.001)) T_old(:,find(abs(x-0.6) < 0.001))...
    T_old(:,find(abs(x-0.8) < 0.001)) T_old(:,find(abs(x-1.0) < 0.001))])

[Torsten]

1. Delete the double loop before the while loop. T_old is already initialized to 0 in the interior.
2. You now programmed Jacobi iteration in the while loop, not Gauss-Seidel (also the T_old's on the right-hand side of the equation have to be replaced by T_new.
3. a = round(a,3) rounds a to 3 decimal places.
4. I don't know how you are told to calculate Errormax. Maybe the relative error is better suited than the absolute error: Errormax = max(abs((T_old-T_new)./T_old,[],'all'))

[mehmet salihi]

dear, now it is ok for the matrix. but about round, not coming as i want, here i put the table to be printed.

round does not work as i wish

%%% Point Gauss Seidel (PGS)  

close all, clc
format short
%Defining the constants
L=1;                                % length
H=2;                                % Height
deltax=0.05
deltay=0.05
Beta=deltax/deltay
Beta2=Beta^2
jn=H/deltay                         % Maximum number of grid points along y
im=L/deltax                         % Maximum number of grid points along x
T1=100; T2=0; T3=0; T4=0;           % boundary conditions
y=2:-deltay:0;
x=0:deltax:L;

% initialize T_old to the initial guess values
T_old=zeros(jn+1,im+1);

% set boundary conditions
T_old(1,1:im+1)=T1;
T_old(jn+1,1:im+1)=T3;
T_old(2:jn+1,1)=T2;
T_old(2:jn+1,im+1)=T4;

T_new = T_old;

eps = 1e-5;
Errormax=0.01;
t=0;
while Errormax > eps
    t=t+1;
    for i=2:jn 
        for j=2:im
            T_new(i,j)=(1/(2*(1+Beta2)))* (T_new(i-1,j)+T_new(i+1,j)+Beta2*(T_new(i,j+1)+T_new(i,j-1)) ) ;       
        end
    end

Errormax = max(abs(T_old-T_new),[],'all');
T_old = T_new;
end

T_new=round(T_new,3);
T_new=flip(T_new); 
disp(' ');disp(' ');
disp([y' T_new(:,find(abs(x-0.0) < 0.001)) T_new(:,find(abs(x-0.2) < 0.001))...
    T_new(:,find(abs(x-0.4) < 0.001)) T_new(:,find(abs(x-0.6) < 0.001))...
    T_new(:,find(abs(x-0.8) < 0.001)) T_new(:,find(abs(x-1.0) < 0.001))])

[Torsten]

Why not ? Three decimal places everywhere ...

[mehmet salihi]

answer for your comment

1. I don't know how you are told to calculate Errormax. Maybe the relative error is better suited than the absolute error: Errormax = max(abs((T_old-T_new)./T_old,[],'all'))

could you please rechange the code error condition according to this. becasue when i use eps = 1e-5;

it gives 647 iterations which should be 574.

and about the round, when i use format short it gives 4 decimal with round command and when i use format short g it write 3 decimal but the integers also in 3 decimal. (100.000 to be 100 only)

thanks in advance for your help

[Torsten]

I think it should be no problem for you to do the remaining cosmetic changes on your own.

[mehmet salihi]

actually i tried by myself before asking you, but not worked, the iteration comes 0 or 1.

[Torsten]

Please show your actual code.

What do you mean by "the iteration comes 0 or 1" ?

[mehmet salihi]

now does not gove 0 or 1 but

please i write a note in error loop, read it.

%%% Point Gauss Seidel (PGS)  
close all, clc
format short g
%Defining the constants
L=1;                                % length
H=2;                                % Height
deltax=0.05
deltay=0.05
Beta=deltax/deltay
Beta2=Beta^2
jn=H/deltay                         % Maximum number of grid points along y
im=L/deltax                         % Maximum number of grid points along x
T1=100; T2=0; T3=0; T4=0;           % boundary conditions
y=2:-deltay:0;
x=0:deltax:L;
% initialize T_old to the initial guess values
T_old=zeros(jn+1,im+1);
% set boundary conditions
T_old(1,1:im+1)=T1;
T_old(jn+1,1:im+1)=T3;
T_old(2:jn+1,1)=T2;
T_old(2:jn+1,im+1)=T4;
T_new = T_old;
% eps = 1e-5;   %3.1 e-5
Errormax=0.01;
t=0;
while Error < Errormax
    t=t+1;
    for i=2:jn 
        for j=2:im
            T_new(i,j)=(1/(2*(1+Beta2)))* (T_new(i-1,j)+T_new(i+1,j)+Beta2*(T_new(i,j+1)+T_new(i,j-1)) ) ;       
        end
    end
    
        for i=2:jn 
            for j=2:im
                % when i make rhs errormax it works but does not give
                % correct iteration number
                Error = abs(T_old(i,j)-T_new(i,j));     
            end
        end

T_old = T_new;
end

T_new=round(T_new,3); T_new=flip(T_new); 
disp(' ');disp(' ');
disp([y' T_new(:,find(abs(x-0.0) < 0.001)) T_new(:,find(abs(x-0.2) < 0.001))...
    T_new(:,find(abs(x-0.4) < 0.001)) T_new(:,find(abs(x-0.6) < 0.001))...
    T_new(:,find(abs(x-0.8) < 0.001)) T_new(:,find(abs(x-1.0) < 0.001))])
disp(['number of iteration is ',int2str(t)])

[Torsten]

1. You did not give a value to Error before entering the while loop (should of course be larger than Errormax)
2. while Error > Errormax instead of while Error < Errormax
3. error = sum(sum(abs(T_old-T_new),2)) instead of your loop

[mehmet salihi]

i did the step what you told, and the itteration is 0.

could not you work on my code with your pc on matlab to check it. really will appretiate it.

thanks.

but as the bool told us Error < Errormax this have to be applied.

Error = 1;   %3.1 e-5
Errormax=0.01;
t=0;
 
while Error < Errormax
    t=t+1;
    for i=2:jn 
        for j=2:im
            T_new(i,j)=(1/(2*(1+Beta2)))* (T_new(i-1,j)+T_new(i+1,j)+Beta2*(T_new(i,j+1)+T_new(i,j-1)) ) ;       
        end
    end
    
        for i=2:jn 
            for j=2:im
                Error = sum(sum(abs(T_old-T_new)),2);     
            end
        end

T_old = T_new;
end

[Torsten]

Then the book is wrong.

And you will have to delete the double loop around the statement

Error = sum(sum(abs(T_old-T_new),2))

[mehmet salihi]

Error = 1;

ok now it is ok. it works good.

thanks very much.

yes the book is wrong, because ( for Error > Errormax case) when i set Error to 1 it gives correct iteration, and when i give Error 1000 it gives same itteration.

[Jan]

@mehmet salihi: Please use the button to format your code. This improves the readability. I've sone this tody for you.

        for i=2:jn 
            for j=2:im
                Error = sum(sum(abs(T_old-T_new)),2);     
            end
        end        

This piece of code repeats the calculation of Error jn*im times. This is a waste of time. Omit the loops.

[mehmet salihi]

yes dear @Jan, i removed the double loop,

ok i find the button. thanks