Trying to plot number of iterations vs gridsize for steady state

4 views (last 30 days)
I need to compute the number of iterations taken to reach the steady-state within a given tolerance. I then need to plot the final number of iterations against the grid size but I cannot for the life of me work out how I am supposed to do that. How am I suppose to set up and then plot those two things?
Code:
n=2;
nx=2^n;
nz=nx;
a=25;
b=25;
x = linspace(0, a, nx);
z = linspace(0, b, nz);
[X, Z] = meshgrid(x,z);
Tnp1 = zeros(nx, nz);
Tnp1(:,1) = 20;
Tnp1(:,end) = 20;
Tnp1(1,:) = 20+380*sin((X.*pi)/25)+205*sin((X.*5*pi)/25);
Tnp1(end,:) = 20;
err = 1;
tol = 1e-8;
k=0;
while err > tol
Tn = Tnp1;
k=k+1;
for i = 2:nx-1
for j = 2:nz-1
Tnp1(i,j) = (1/4)*(Tn(i+1,j)+Tn(i-1,j)+Tn(i,j+1)+Tn(i,j-1));
end
end
err = max(abs(Tnp1(:) - Tn(:)));
end
T = Tnp1;
plot(X,k)
end
end

Answers (1)

SALAH ALRABEEI
SALAH ALRABEEI on 5 Jun 2021
Edited: SALAH ALRABEEI on 5 Jun 2021
Hi, this looks like a diffusion equation solved by FDM. you had to mistakes, I already corrected them, n is the step size or grid size, and the u should have used x not X in the boundary condition conditions. See the code below
clear N = [2,4]; for kk = 1: length(N) n=N(kk); nx=2^n; nz=nx; a=25; b=25; x = linspace(0, a, nx); z = linspace(0, b, nz); [X, Z] = meshgrid(x,z); Tnp1 = zeros(nx, nz); Tnp1(:,1) = 20; Tnp1(:,end) = 20; Tnp1(1,:) = 20+380*sin((x.*pi)/25)+205*sin((x.*5*pi)/25); Tnp1(end,:) = 20; err = 1; tol = 1e-8; k=0; while err > tol Tn = Tnp1; k=k+1; for i = 2:nx-1 for j = 2:nz-1 Tnp1(i,j) = (1/4)*(Tn(i+1,j)+Tn(i-1,j)+Tn(i,j+1)+Tn(i,j-1)); end end err = max(abs(Tnp1(:) - Tn(:))); end T = Tnp1; %subplot(2,2,kk) %surf(X,Z,T),grid %title([' Profile at grid = ' num2str(kk)])
Num_iter(kk) = k; end
figure(2) plot(N,Num_iter),grid
  2 Comments
Jayden Cavanagh
Jayden Cavanagh on 5 Jun 2021
Edited: Jayden Cavanagh on 5 Jun 2021
It is a steady-state temperature field in a 2D rectangle solved using G-S method. So it isn't exactly what I am looking for. For my problem I have to choose the number of divisions in the domain to be ndiv = 2^n for n = 2 : 9 in both the x and z directions. n=2 was just a way to reduce processing power when testing. While your code works for N=[2,4] I get error code
Requested array exceeds the maximum possible variable size.
Error in linspace (line 44)
y = d1 + (0:n1).*(d2 - d1)./n1;
Error in hope (line 9)
x = linspace(0, a, nx);
when I try and run it for the maximum value of n=8 and nx=2^8. Any way the code could be changed to be able to composate for that?
SALAH ALRABEEI
SALAH ALRABEEI on 5 Jun 2021
Edited: SALAH ALRABEEI on 5 Jun 2021
%
clear
N = [2,4];
for kk = 1: length(N)
n=N(kk);
nx=2^n;
nz=nx;
a=25;
b=25;
x = linspace(0, a, nx);
z = linspace(0, b, nz);
[X, Z] = meshgrid(x,z);
Tnp1 = zeros(nx, nz);
Tnp1(:,1) = 20;
Tnp1(:,end) = 20;
Tnp1(1,:) = 20+380*sin((x.*pi)/25)+205*sin((x.*5*pi)/25);
Tnp1(end,:) = 20;
err = 1;
tol = 1e-8;
k=0;
while err > tol
Tn = Tnp1;
k=k+1;
for i = 2:nx-1
for j = 2:nz-1
Tnp1(i,j) = (1/4)*(Tn(i+1,j)+Tn(i-1,j)+Tn(i,j+1)+Tn(i,j-1));
end
end
err = max(abs(Tnp1(:) - Tn(:)));
end
T = Tnp1;
%subplot(2,2,kk)
%surf(X,Z,T),grid
%title([' Profile at grid = ' num2str(kk)])
Num_iter(kk) = k;
end
figure(2)
plot(N,Num_iter),grid

Sign in to comment.

Categories

Find more on Performance and Memory in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!