How Can I Vectorize My MATLAB Code to Increase Speed in Solving a Heat Equation?

[range.rover]

Hi guys,I am new for programming. I made my first attempt to solve an heat equation by finite difference method and wrote a code for it in Matlab.I got a solution but i need help from you guys to (vectorize) increase the speed of my program.

((∂^2 T)/(∂x^2 )+(∂^2 T)/(∂y^2 )+(∂^2 T)/(∂z^2 ))(k)=-Q(x,y,z)

clear;
close all;
clc;
n = 5;% n grids and has n - 2 interior points per dimension

x = linspace(0,1,n);% line in x direction of 0 to 1 and divided into n parts

dx = x(2)-x(1); % distance between grids

y = x;

dy = dx;

k=1; %thermal conductivity

Q=300 ; % heat source
m=(n-2)^2;

%Right hand side matrix

B=zeros(m,1);

if (mod(m,2)==0)
B(m/2,1)=Q;
else
B((m+1)/2,1)=Q;
end

%Main Matrix

Tsol = zeros(n);
%boundry conditions
Tsol(1,1:n) = 0; %TOP
Tsol(n,1:n) = 0; %BOTTOM
Tsol(1:n,1) = 0; %LEFT
Tsol(1:n,n) = 0; %RIGHT

tic
error = 1; z = 0;
while error > 0.000001
z = z+1;
c=1;
Torg = Tsol;
for i = 2:n-1
for j = 2:n-1
Tsol(i,j) =abs(((dx^2*B(c,1))+((Torg(i+1,j)+Torg(i-1,j)+Torg(i,j+1)+Torg(i,j-1))*k))/(4*k));
c=c+1;
end
end

error = max(max(abs(Torg-Tsol)));
end


toc

%% plotting

figure, contour(x,y,Tsol),
title('Temperature (Steady State)'),xlabel('x'),ylabel('y'),colorbar



figure,pcolor(x,y,Tsol),shading interp,
title('Temperature (Steady State)'),xlabel('x'),ylabel('y'),colorbar

 

[kreil]

You initialize Tsol as a matrix of zeros, yet the next 4 lines set different boundary values to 0? Seems unnecessary.

Can't you just use del2() or gradient() here?

 

[AlephZero]

You initialize Tsol as a matrix of zeros, yet the next 4 lines set different boundary values to 0? Seems unnecessary.

"Initializing the solution" and "setting the boundary conditions" are two different steps in the solution algorithm. For this particular equation they happen to have the same numerical effect (setting some terms to zero twice) but optimizing that away would have a negligible effect on the total execution time.

To optimize the code, you first have to find out what is taking the time. For a small problem (n = 5) you only have about (n-2)² = 9 calculations of the elements of Tsol, in each iteration. So if your program runs for a long time (e.g. more than 1 second) you must be doing a huge number of iterations. Try counting the number of times the "while error > 0.000001" loop runs, and think about how to reduce that to a more "reasonable" number (e.g. 100 or less).

Converting the "Tsol" calculation to an array operation will help, but if you want to optimize a program effectively, you first need to find where the real bottlenecks are, not just optimize bits of code "because you can".

 

[range.rover]

I have already checked my program and found that my ittirations are going way beond than 100000,it is the reason why it is taking long time. I am trying a lot to vectorize my while loop but i am not able to do because of my lack of knowledge and experience in programming skills.Thats why i needed a helping hand,i will be glad to take your help in my first experience.

error = 1; z = 0;
while error > 0.000001
z = z+1;
c=1;
Torg = Tsol;
for i = 2:n-1
for j = 2:n-1
Tsol(i,j) =abs(((dx^2*B(c,1))+((Torg(i+1,j)+Torg(i-1,j)+Torg(i,j+1)+Torg(i,j-1))*k))/(4*k));
c=c+1;
end
end

error = max(max(abs(Torg-Tsol)));
end

 

[AlephZero]

Print out value of "error" in the first terations (say the first 100). Is it converging monotonically, or oscillating, or jumping about "randomly"?

 

[range.rover]

the error is exponentially decreasing and some times it is raising up at certain values of n,like you can see in the images down below at n = 25.

 

[AlephZero]

Very strange.

I don't understand the first value of "error". On the first iteration everything in the calculation of Tsol should be 0, apart from one value in the B array which is 300. dx = 0.25. So the only non-zero value in tsol should be
300 * (0.25)^2 / 4 = about 4.7. But on your plots it is about 0.67, not 4.7.

Maybe I misunderstood what your code does, or maybe you are not running the same code that you posted here.

Incidentally, why is the abs(...) function in the calculation of tsol?

 

[kreil]

OK, I invested some time and looked into this.

While it's easy to vectorize the for loops in a general sense, you were using the c variable as a sort of delta function to pick out the value from B for the heat source. This means that your loops were doing nothing for the first m/2 iterations, where the heat source is finally picked up and starts to propagate. So you can get rid of those m/2 unnecessary iterations by just initializing the heat source and adding it in differently. I did this by reshaping B into the interior of a matrix, A. Since A is the same size as Tsol, you can use it to initialize Tsol and then add in the heat source after each iteration. This means there are no wasted iterations; the very first one is productive.

Using MATLAB R2014a with n=25, the time differences on my laptop were 0.101265 seconds for your code, and 0.038494 seconds for the code below, so the improvement is significant. I believe this could get even more dramatic for larger n.

Hopefully this helps!

Note: While I still think del2 should be useful here (since it calculates vectorized finite differences via bsxfun), it uses single sided differences on the edges of the matrix, which messes with your boundary conditions.

n = 25;% n grids and has n - 2 interior points per dimension
x = linspace(0,1,n);% line in x direction of 0 to 1 and divided into n parts
dx = x(2)-x(1); % distance between grids
y = x;
dy = dx; 
k=1; %thermal conductivity
Q=300 ; % heat source
m=(n-2)^2;

%Right hand side matrix
B=zeros(m,1);
if (mod(m,2)==0)
B(m/2,1)=Q;
else
B((m+1)/2,1)=Q;
end

%% Main Matrix

%Initialize a heat source matrix
A = zeros(n);
A(2:n-1,2:n-1) = reshape(B,n-2,n-2);
A = (dx^2*A)/(4*k);

Tsol = A; 
%boundary conditions
Tsol(1,1:n) = 0; %TOP
Tsol(n,1:n) = 0; %BOTTOM
Tsol(1:n,1) = 0; %LEFT
Tsol(1:n,n) = 0; %RIGHT

tic
error = 1; 
z = 0;
i = 2:n-1;
j = 2:n-1;
while error > 0.000001
    z = z+1;
    Torg = Tsol;
    Tsol(i,j) = abs(((Torg(i+1,j)+Torg(i-1,j)+Torg(i,j+1)+Torg(i,j-1))*k))/(4*k);
    Tsol = Tsol + A;
    error = max(max(abs(Torg-Tsol)));
end
toc

figure, contour(x,y,Tsol),
title('Temperature (Steady State)'),xlabel('x'),ylabel('y'),colorbar
figure,pcolor(x,y,Tsol),shading interp,
title('Temperature (Steady State)'),xlabel('x'),ylabel('y'),colorbar