How can I make my code run faster?
1 view (last 30 days)
Show older comments
a = 1;
b = 1;
rho = 1;
D = 3;
h =1;
%set matrix phi
Nx=150;
Ny=100;
phi = rand(Nx,Ny);
%assign random -1 and 1 for phi matrix
for i=1:Nx
for j=1:Ny
if phi(i,j)<=0.5
phi(i,j)=-1;
else
phi(i,j)=1;
end
end
end
Problem = input ('problem number');
switch Problem
case 1
%Five Point Stencil
dt = 1e-3;
for k = 1:dt:20
% apply Laplacian to matrix phi
v = Laplacian_2D (phi,h,5);
% Find MAtrix A in equation 4
A = ((b^4*phi.^3 - a*b^2*phi)-rho*v);
% Calculate Laplacian of Matrix A for equation 4
%left hand of equation 4
E4 = D * Laplacian_2D (A,h,5);
%First Order Runge Kutta with Equation 4
c1 = dt*E4;
phi = phi + c1;
%draw plot
imagesc (phi);
colorbar;
caxis ([-1,1])
title ('Five Point Stencil: First Order Runge Kutta');
drawnow;
end
case 2
%Nine Point Stencil
dt = 1e-4;
for k = 1:dt:(20/dt)
counter = 0;
% apply Laplacian to matrix phi
v = Laplacian_2D (phi,h,9);
% Find MAtrix A in equation 4
A = ((b^4*phi.^3 - a*b^2*phi)-rho*v);
% Calculate Laplacian of Matrix A for equation 4
v2 = Laplacian_2D (A,h,9);
%left hand of equation 4
E4 = D * v2;
%First Order Runge Kutta with Equation 4
c1 = dt*E4;
phi = phi + c1;
%draw plot
imagesc (phi);
colorbar;
caxis ([-1,1])
title ('Nine Point Stencil: First Order Runge Kutta');
drawnow;
end
case 3
%Second Order Runge Kutta with Five Point Stencil
dt = 5e-3;
for k = 1:dt:(20/dt)
% apply Laplacian to matrix phi
v = Laplacian_2D (phi,h,5);
% Find MAtrix A in equation 4
A = ((b^4*phi.^3 - a*b^2*phi)-rho*v);
% Calculate Laplacian of Matrix A for equation 4
v2 = Laplacian_2D (A,h,5);
%left hand of equation 4
E4 = D * v2;
%run again with c1
c1 = dt*E4;
v3 = Laplacian_2D (phi+0.5*c1,h,5);
R = ((b^4*(phi+0.5*c1).^3 - a*b^2*(phi+0.5*c1))-rho*v3);
v4 = Laplacian_2D (R,h,5);
E4_2 = D * v4;
%Second Order Runge Kutta with Equation 4
c2 = dt*E4_2;
phi = phi + c2;
%draw plot
imagesc (phi);
colorbar;
caxis ([-1,1])
title ('Five Point Stencil: Second Order Runge Kutta');
drawnow;
end
end
function v = Laplacian_2D(u,h,stencil)
Nx = 150;
Ny = 100;
v = zeros (Nx,Ny);
if stencil == 5 %insert for the five point stencil
%periodic boundary conditions
for i = 1:Nx
N = i-1;
if N==0
N=Nx;
end
S = i + 1;
if S== Nx+1
S=1;
end
for j = 1:Ny
E = j+1;
if E == Ny+1
E = 1;
end
W = j-1;
if W == 0
W=Ny;
end
v (i,j) = (1/h^2)*(u(N,j) + u(S,j) + u(i,W) + u(i,E) - 4*u(i,j));
end
end
elseif stencil == 9 %insert for the nnine point stencil
%periodic boundary conditions
for i = 1:Nx
N = i-1;
if N==0
N =Nx;
end
S = i + 1;
if S== Nx+1
S = 1;
end
for j = 1:Ny
E = j+1;
if E == Ny+1
E = 1;
end
W = j-1;
if W == 0
W=Ny;
end
v (i,j) = (1/(6*h^2))*(u(N,W)+ 4*u(N,j) + u(N,E) + 4*u(i,W) + 4*u(i,E) + u(S,W) + 4*u(S,j) + u(S,E) - 20*u(i,j));
end
end
else
error ('stencil must be 5 or 9');
end
end
3 Comments
Adam
on 15 Aug 2019
First step in working out how to make code faster is to understand what is slow, as Rik says.
Answers (1)
darova
on 15 Aug 2019
Shorter version of Laplacian_2D function. Is there any connection to gradient and laplacian built-n functions?
un = circshift(u, [-1 0]);
us = circshift(u, [+1 0]);
uw = circshift(u, [0 -1]);
ue = circshift(u, [0 +1]);
unw = circshift(u, [-1 -1]);
une = circshift(u, [-1 +1]);
usw = circshift(u, [+1 -1]);
use = circshift(u, [+1 +1]);
if stencil == 5
v = 1/h^2*(un + us + uw + ue - 4*u);
elseif stencil == 9 %insert for the nnine point stencil
v = 1/(6*h^2)*(4*(un+uw+ue+us) + unw+une+usw+use - 20*u);
else
error ('stencil must be 5 or 9');
end
This part can be shorter
% phi = rand(Nx,Ny);
% %assign random -1 and 1 for phi matrix
% for i=1:Nx
% for j=1:Ny
% if phi(i,j)<=0.5
% phi(i,j)=-1;
% else
% phi(i,j)=1;
% end
% end
% end
phi = randi([0 1],Nx,Ny)*2-1;
0 Comments
See Also
Categories
Find more on Startup and Shutdown 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!