Hi.. this program contains no errors, I do not know why the compilaton is long !!!

clc
clear   
format short
% Define the parameters
eps = 0.2; % Perturbation Parameter
% Define the number of grid points in x and y direction
% nx = number of grid points in x direction
% ny = number of grid points in y direction
nx = 101;
ny = 65;
t = (2.0 * cos(pi/(nx * ny)))^2;
% Calculate a optimum value of SOR parameter
omega1 = (16.0+sqrt((256.0-(64.0 * t))))/(2.0 * t);
omega2 = (16.0-sqrt((256.0-(64.0 * t))))/(2.0 * t);
oopt = min(omega1,omega2);
if ( (oopt <= 1.0) || (oopt >= 2.0 ))
oopt = 1.0;
end
omega = oopt;
% Define the dimension of the domain
Lx = 1.6; % Length in x of the computation region
Ly = 1.0; % Length in y of the computation region
% Calculate the mesh size
hx = Lx/(nx-1); % Grid spacing in x
hy = Ly/(ny-1); % Grid spacing in y
% Generate the mesh/grid
x=zeros(nx,1);
for i = 2:nx
x(i)=x(i-1)+hx;
end
y=zeros(ny,1);
for j = 2:ny
y(j)=y(j-1)+hy; 
end
% Solve the Temperature equation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize the temperature field as zero
T = zeros(nx,ny);
% Max-Norm on Error {L-inf Error) Initialized
Linf = 1.0;
iteration = 0;
while (Linf > 1.0e-5)
iteration = iteration+1;
% Store the old values of T in Told
Told = T;
% Apply the boundary conditions
for i = 1:nx
% BC for the bottom boundary
T(i,1)=T(i,2)+hy;
% BC for the top boundary
if ( x(i) <= 0 )
T(i,ny) = 1;
else
T(i,ny) = 0;
end
end
for j = 1:ny
% BC for the Left boundary
T(1,j)=T(2,j);
% BC for the Left boundary
T(nx,j)=T(nx-1,j);
end
% Now compute the interior domain using 2nd order finite difference
for i = 2:nx-1
for j = 2:ny-1
term1 = ((eps/hx)^2) * (T(i-1,j)+T(i+1,j));
term2 = ((1.0/hy)^2) / (T(i,j-1)+T(i,j+1));
num = term1+term2;
den = 2.0 * (((eps/hx)^2)+((1.0/hy)^2));
Tgs = num/den;
T(i,j) = (omega * Tgs)+((1.0-omega) * Told(i,j));
end
end
% Calculate the error
Terr = abs(T-Told);
Linf = norm(Terr,2);
% Print the convergence history every 50 iterations
ccheck = round(iteration/100)-(iteration/100);
if ( (iteration == 1) || ( ccheck == 0) )
fprintf('%d \t %e \n', iteration, Linf)
end
end
% Plot the solutions
figure 
[X,Y] = meshgrid(x,y);
clevel = [-1 -0.005 0 0.053 0.152 0.252 0.352 0.451 0.551 0.65 0.75 0.949 1.049 1.148 1.248 1.348 1.447 1.547 1.646 1.746 1.848 1.945 2.0 3.0];
contourf(X',Y',T,clevel)
colorbar 
axis([-1.2 1.2 -0.8 1.8])
xlabel('x')
ylabel('y')
zlabel('T(x,y)')
hold on
yi = (0:0.01:1);
xi = zeros(length(yi));
plot(xi,yi,'black')
print -djpeg figures\chap4\temperature4

[OCDER]

You mean running it? or compiling it into a binary file? Compiling could be slow slow as matlab has to check for dependencies, convert, etc.

Also, can you do us big favor and hit that {} Code button, and explain how long it takes to compile?

[yassine amadane]

I mean running ...

[OCDER]

Oh ok. One last thing - highlight your question text, and push the {} Code button located above the text box you're using, so that we can see your code

Like this
for k = 1:100
end

[Walter Roberson]

Are you asking why your code is taking a long time to converge? Or are you asking for hints on vectorizing your code?

Answer 1

[Chad Greene]

I formatted the code to make it easier for folks to read.

Using the profiler for a few iterations indicates most of the time is spent inside the nested loops. I stopped the code on its 443rd time through the while loop, but by then, term1, term2, etc had already been calculated over 2.7 million times. I suggest rewriting the code ( vectorize) to remove these nested loops:

for i = 2:nx-1
for j = 2:ny-1
  ... 
end
end

[OCDER]

Thanks Chad for formatting the code! Yassine, the profiler that Chad is suggesting is really helpful for optimizing code, so please do give the document a read. I use it a lot to find bottlenecks quickly, and then the tic/toc or timeit method for a quick check to see if a new code is faster/slower.

Answer 2

[OCDER]

I'm getting that norm is the slowest step, follow closely by that nested for loop that Chad is pointing out. There's a sqrt in norm that is slowing things down. You can try rewriting the while loop to check Linf BEFORE the sqrt is taken, for instance:

 while Linf2 > (1.0e-5)^2

     ....

     Linf2 = sum(sum(T-Told).^2)); %Calc the norm^2 yourself, but careful about precision as it's not 100% the same as norm(T-Told, 2).
 end

This sped things up 26 times for me, but then the nested for loop becomes the bottleneck. To be even faster, this article could be helpful: https://www.mathworks.com/company/newsletters/articles/programming-patterns-maximizing-code-performance-by-optimizing-memory-access.html