This post is about comparing 2 codes to solve the 2D Laplace equation using finite difference method. A sample code mentioned under "Code 01" uses nested loops. We must, wherever possible avoid nested loops. The solution to the code is shown in Fig. 1. The second code uses vectorization instead of the nested loops. The vectorized version is mentioned under "Code 02".
Code 01
clear
clc
close all
%% Parameters
Lx = 1; % Length of the domain in the x-direction
Ly = 1; % Length of the domain in the y-direction
Nx = 201; % Number of grid points in the x-direction
Ny = 201; % Number of grid points in the y-direction
dx = Lx / (Nx - 1); % Grid spacing in the x-direction
dy = Ly / (Ny - 1); % Grid spacing in the y-direction
%% Initialize temperature matrix
T = zeros(Nx, Ny);
T(:, 1) = 100;
T(:, Nx) = 0;
T(1, :) = 25;
T(Ny, :) = 50;
%% Gauss-Seidel iteration
max_iter = 50000; % Maximum number of iterations
tolerance = 1e-10; % Convergence tolerance
error = inf; % Initialize error
iter = 0; % Iteration counter
while error > tolerance && iter < max_iter
T_old = T;
% Solve Laplace equation using Gauss-Seidel iterations
for i = 2:Nx-1
for j = 2:Ny-1
T(i, j) = ((T(i+1, j) + T(i-1, j))*dy^2 + (T(i, j+1) + T(i, j-1))*dx^2) / (2*(dx^2 + dy^2));
end
end
% Compute error
error = max(abs(T(:) - T_old(:)));
% Increment iteration counter
iter = iter + 1;
end
%% plotting
[X, Y] = meshgrid(0:dx:Lx, 0:dy:Ly);
contourf(X, Y, T')
axis equal
colormap jet
colorbar
clim([0 100])
title('Temperature Distribution Non Vec')
xlabel('x')
ylabel('y')
zlabel('Temperature (T)')
colorbar
Code 02
clear
clc
close all
%% Parameters
Lx = 1; % Length of the domain in the x-direction
Ly = 1; % Length of the domain in the y-direction
Nx = 201; % Number of grid points in the x-direction
Ny = 201; % Number of grid points in the y-direction
dx = Lx / (Nx - 1); % Grid spacing in the x-direction
dy = Ly / (Ny - 1); % Grid spacing in the y-direction
i = 2:Nx-1;
j = 2:Ny-1;
%% Initialize temperature matrix
T = zeros(Nx, Ny);
T(:, 1) = 100;
T(:, Nx) = 0;
T(1, :) = 25;
T(Ny, :) = 50;
%% Gauss-Seidel iteration
max_iter = 50000; % Maximum number of iterations
tolerance = 1e-10; % Convergence tolerance
error = inf; % Initialize error
iter = 0; % Iteration counter
while error > tolerance && iter < max_iter
T_old = T;
% Solve Laplace equation using Gauss-Seidel iterations
T(i, j) = ((T(i+1, j) + T(i-1, j))*dy^2 + (T(i, j+1) + T(i, j-1))*dx^2) / (2*(dx^2 + dy^2));
% Compute error
error = max(abs(T(:) - T_old(:)));
% Increment iteration counter
iter = iter + 1;
end
%% plotting
[X, Y] = meshgrid(0:dx:Lx, 0:dy:Ly);
contourf(X, Y, T')
axis equal
colormap jet
colorbar
clim([0 100])
title('Temperature Distribution Vec')
xlabel('x')
ylabel('y')
zlabel('Temperature (T)')
colorbar
Result
Results say that vectorized code is ~1.5x faster than nested looped code for 200x200 matrix. Simulation results are now presented.
 |
| Fig. 1, Vectorized VS Nested Loops |
Thank you for reading! I hope you learned something new! If you like this blog and want to hire me as your PhD student, please get in touch!
