Matrix calculation gets slower after each iteration in matlab

Question

I have a 1024*1024*51 matrix. I'll do calculations to change some value of the matrix within for loops (change the value of matrix for each iteration). I find that the computing speed gets slower and slower and finally my computer gets into trouble.However the size of the matrix doesn't change. Anyone can shed some light on this problem?

function ActiveContours3D(method,grad,im,mu,nu,lambda1,lambda2,TimeSteps)
epsilon = 10e-10;
tic
fid=fopen('Chr18_z_25of25tiles-C=0_c0_n000.raw','rb','ieee-le');   
Xdim = 1024;
Ydim = 1024;
Zdim = 51;
A = fread(fid,[Xdim Ydim*Zdim],'int16');
A = double(A);
size_of_A = size(A)
for(i=1:Zdim)
    u0_color(:,:,i) = A(1 : Xdim , (i-1)*Ydim+1 : i*Ydim);
end
fclose(fid)

time = toc

[M,N,P,color] = size(u0_color);
size(u0_color );
u0_color = double(u0_color);    % Convert u0_color values to double;
u0 = u0_color(:,:,:,1);           % Define the Grayscale volumetric image.
u0_color = uint8(u0_color);     % Necessary for color visualization
x = 1:M;
y = 1:N;
z = 1:P;
dx = 1
dy = 1
dz = 1
dim_approx = 2*M*N*P / sqrt(M*N*P);
if(method == 'Explicit')
    dt = 0.9 / ((2*mu/(dx^2)) + (2*mu/(dy^2)) + (2*mu/(dz^2)))    % 90% CFL
elseif(method == 'Implicit')
    dt = (10e7) * 0.9 / ((2*mu/(dx^2)) + (2*mu/(dy^2)) + (2*mu/(dz^2)))
end
[X,Y,Z] = meshgrid(x,y,z);
x0 = (M+1)/2;
y0 = (N+1)/2;
z0 = (P+1)/2;
r0 = min(min(M,N),P)/3;
phi = sqrt((X-x0).^2 + (Y-y0).^2 + (Z-z0).^2) - r0; 
phi_visualize = phi;            % Use this for visualization in 3D
phi = permute(phi,[2,1,3]);     % Use this for computations in 3D
write_to_binary_file(phi_visualize,0);       % record initial conditions

tic

for(n=1:TimeSteps)
n

c1 = C1_3d(u0,phi);
c2 = C2_3d(u0,phi);

% x
phi_xp = [phi(2:M,:,:); phi(M,:,:)];      % vertical concatenation
phi_xm = [phi(1,:,:);   phi(1:M-1,:,:)];        % (since x values are rows)
                        % cat(1,A,B) is the same as [A;B]
Dx_m = (phi - phi_xm)/dx;           % first derivatives
Dx_p = (phi_xp - phi)/dx;

Dxx = (Dx_p - Dx_m)/dx;     % second derivative

% y
phi_yp = [phi(:,2:N,:) phi(:,N,:)];       % horizontal concatenation
phi_ym = [phi(:,1,:)   phi(:,1:N-1,:)];         % (since y values are columns)
                        % cat(2,A,B) is the same as [A,B]
Dy_m = (phi - phi_ym)/dy;
Dy_p = (phi_yp - phi)/dy;

Dyy = (Dy_p - Dy_m)/dy;

% z
phi_zp = cat(3,phi(:,:,2:P),phi(:,:,P));
phi_zm = cat(3,phi(:,:,1)  ,phi(:,:,1:P-1));

Dz_m = (phi - phi_zm)/dz;
Dz_p = (phi_zp - phi)/dz;

Dzz = (Dz_p - Dz_m)/dz;

% x,y,z
Dx_0 = (phi_xp - phi_xm) / (2*dx);
Dy_0 = (phi_yp - phi_ym) / (2*dy);
Dz_0 = (phi_zp - phi_zm) / (2*dz);

phi_xp_yp = [phi_xp(:,2:N,:) phi_xp(:,N,:)];
phi_xp_ym = [phi_xp(:,1,:)   phi_xp(:,1:N-1,:)];
phi_xm_yp = [phi_xm(:,2:N,:) phi_xm(:,N,:)];
phi_xm_ym = [phi_xm(:,1,:)   phi_xm(:,1:N-1,:)];

phi_xp_zp = cat(3,phi_xp(:,:,2:P),phi_xp(:,:,P));
phi_xp_zm = cat(3,phi_xp(:,:,1)  ,phi_xp(:,:,1:P-1));
phi_xm_zp = cat(3,phi_xm(:,:,2:P),phi_xm(:,:,P));
phi_xm_zm = cat(3,phi_xm(:,:,1)  ,phi_xm(:,:,1:P-1));

phi_yp_zp = cat(3,phi_yp(:,:,2:P),phi_yp(:,:,P));
phi_yp_zm = cat(3,phi_yp(:,:,1)  ,phi_yp(:,:,1:P-1));
phi_ym_zp = cat(3,phi_ym(:,:,2:P),phi_ym(:,:,P));
phi_ym_zm = cat(3,phi_ym(:,:,1)  ,phi_ym(:,:,1:P-1));    


if(grad == 'Dirac')
    Grad = DiracDelta(phi);                  % Dirac delta    
    %Grad = 1;
elseif(grad == 'Grad ')
    Grad = (((Dx_0.^2)+(Dy_0.^2)+(Dz_0.^2)).^(1/2));   % |grad phi|
end

if(method == 'Explicit')
    % CURVATURE:   *mu*k|grad phi|*   (central differences):
    K = zeros(M,N,P);

    Dxy = (phi_xp_yp - phi_xp_ym - phi_xm_yp + phi_xm_ym) / (4*dx*dy);
    Dxz = (phi_xp_zp - phi_xp_zm - phi_xm_zp + phi_xm_zm) / (4*dx*dz);
    Dyz = (phi_yp_zp - phi_yp_zm - phi_ym_zp + phi_ym_zm) / (4*dy*dz);

    K = ( (Dx_0.^2).*Dyy - 2*Dx_0.*Dy_0.*Dxy + (Dy_0.^2).*Dxx ...
        + (Dx_0.^2).*Dzz - 2*Dx_0.*Dz_0.*Dxz + (Dz_0.^2).*Dxx ... 
        + (Dy_0.^2).*Dzz - 2*Dy_0.*Dz_0.*Dyz + (Dz_0.^2).*Dyy) ./ ((Dx_0.^2 + Dy_0.^2 + Dz_0.^2).^(3/2) + epsilon);

    phi_temp = phi + dt * Grad .* ( mu.*K + lambda1*(u0 - c1).^2 - lambda2*(u0 - c2).^2 );

elseif(method == 'Implicit')
    C1x = 1 ./ sqrt(Dx_p.^2 + Dy_0.^2 + Dz_0.^2 + (10e-7)^2);
    C2x = 1 ./ sqrt(Dx_m.^2 + Dy_0.^2 + Dz_0.^2 + (10e-7)^2);
    C3y = 1 ./ sqrt(Dx_0.^2 + Dy_p.^2 + Dz_0.^2 + (10e-7)^2);
    C4y = 1 ./ sqrt(Dx_0.^2 + Dy_m.^2 + Dz_0.^2 + (10e-7)^2);
    C5z = 1 ./ sqrt(Dx_0.^2 + Dy_0.^2 + Dz_p.^2 + (10e-7)^2);
    C6z = 1 ./ sqrt(Dx_0.^2 + Dy_0.^2 + Dz_m.^2 + (10e-7)^2);

    % m = (dt/(dx*dy)) * Grad .* mu;        % 2D
    m = (dt/(dx*dy)) * Grad .* mu;
    C = 1 + m.*(C1x + C2x + C3y + C4y + C5z + C6z);

    C1x_2x = C1x.*phi_xp + C2x.*phi_xm;
    C3y_4y = C3y.*phi_yp + C4y.*phi_ym;
    C5z_6z = C5z.*phi_zp + C6z.*phi_zm;

    phi_temp = (1 ./ C) .* ( phi + m.*(C1x_2x+C3y_4y) + (dt*Grad).*(lambda1*(u0 - c1).^2) - (dt*Grad).*(lambda2*(u0 - c2).^2) );
end

phi = phi_temp;

phi_visualize =  permute(phi,[2,1,3]);
write_to_binary_file(phi_visualize,n);     % record

end

time = toc

n = n
T = dt*n

clear
clear all

Answer 1

In general Matlab keeps track of all the variables in the form of matrix. When you work with lot of variables with many dimensions the RAM memory will be allocated for storing this variable. Hence on working with lots of variables that is gonna run for multiple iterations it is better to clear the variable from the memory. To do so use the command

clear variable_name_1, variable_name_2,... variable_name_3;

Although keeping all the variables keeps the code to look organised, however when you face such issues try clearing the unwanted variables.

Check this link to use clear command in detail: http://www.mathworks.in/help/matlab/ref/clear.html