How can i vectorize the for loop in this code? is there any solution to calculate a matrix having both the variable indice?

clear all;          % clear all variables in memory
close all;

xl = 0;
xr = 1;                                                                
T =0.1;                                                               

Nt = 5000;                                      % Number of time steps
dt = T/Nt;
xn = 50;                                        % Number of space steps
dx = (xr-xl)/xn;
mu = dt/(dx)^2;                                 % Stability parameter (mu=<1)
if mu > 0.5                                                          % make sure dt satisy stability condition
error('mu should < 0.5!')
end
tic;

%  i = 1:xn+1
x(2:xn+1) =(0:xn-1)*dx;
u(2:xn+1,1) =sin(pi*x(2:xn+1)) + sin(2*pi*x(2:xn+1));
   j=2:Nt+1;
  tm(j)= (j-1)*dt;
  t = tm(j);
  u(1,j) = sin(pi*xl)*exp(-pi*pi*t)+sin(2*pi*xl)*exp(-4*pi*pi*t);
  u(xn+1,j) = sin(pi*xr)*exp(-pi*pi*t)+sin(2*pi*xr)*exp(-4*pi*pi*t);
for j=2:Nt+1                                                          % Time Loop
%  i=2:xn;                                                          % Space Loop
u(2:xn,j)=u(2:xn,j-1)+mu*(u(3:xn+1,j-1)-2*u(2:xn,j-1)+u(1:xn-1,j-1));
% end
end

xj = xl + (gpuArray(0:xn).*dx);
B= (sin(pi*xj))';
B_2D = repmat(B,1,Nt+1);
C= exp(-pi*pi*gpuArray(1:Nt+1)*dt);
C_2D = repmat(C,xn+1,1);
u_ex_a = B_2D.*C_2D;
D= (sin(2*pi*xj))';
D_2D = repmat(D,1,Nt+1);
E= exp(-4*pi*pi*gpuArray(1:Nt+1)*dt);
E_2D = repmat(E,xn+1,1);
u_ex_b = D_2D.*E_2D;
u_ex=u_ex_a+u_ex_b;

er = u_ex - u;
D = abs(u-u_ex).^2;
MSE = sum(D(:))/numel(u);

[Jan]

Omit the darn "clear all". It does not only "clear all variables in memory" but removes all loaded functions also, deletes all persistent variables and and classes. Reloading the M-files from disk is a massive waste of time, because each file needs to be parsed again in the first run.

Better use a function to keep your workspace clean.

Answer 1

[Jan]

for j = 2:Nt+1
  u(2:xn,j) = u(2:xn,j-1) + mu * (u(3:xn+1,j-1) - 2*u(2:xn,j-1) + u(1:xn-1,j-1));
end

The loop takes 0.04 seconds on my computer. Are you sure that this is a problem?

[Priya Khot]

Yes, it is correct, but I have to run this code using Parallel Computing toolbox of MATLAB for GPU computing i.e., gpuArray.That code is taking too much time and that is why I have to vectorize it

Answer 2

[Matt J]

Looks like a convolution to me. Use conv2?