Optimize evaluation of triple infinite series using FOR loops

I am programming the analytical solution of a concentration in a rectangular block with constant surface concentrations.

A scan of the analytical solution is attached as a picture.

I need concentrations along three vectors of discreticized x, y, and z and along a vector of time. I am currently evaluating the triple infinite series at a finite number of terms and using FOR loops. A higher number of terms yield more accurate results but the runtime starts to be impractical after 80 terms. For the results I want, I need 100 terms at least.

Below is my current code:

I made a small modification to the computation of alpha by multiplying the diffusion coefficient D by the tortuorsity tau_l and dividing them by the retardation factor R_l in my code. The pencil marks show how I split the programming of the terms inside the triple infinite series into variables p1vec, p2vec, and p3vec in my code.

%INPUTS
a = 5; b = 2; c = 4; %half fracture length
inf_end = 100; %at 80: 155 sec to complete; at 90: 272; at 100: over 4 mins 
nx = 21; ny = 9; nz = 17; 
%COMPUTE: 
l = 0:inf_end; m = 0:inf_end; n = 0:inf_end; 
x = linspace(-a,a,nx); y = linspace(-b,b,ny); z = linspace(-c,c,nz);
t = [.75 1 3]; dx = 2*a/(nx - 1); dy = 2*b/(ny - 1); dz = 2*c/(nz - 1); 
%Units: m2/yr, g/l or kg/m3, unitless, unitless
D = 0.0315576; C0 = 0.1; tau_l = 0.1; R_l = 1; kappa = D*tau_l/R_l;
%% Vectorized approach - compute concentrations
[L,M,N] = meshgrid(l,m,n); 
%non-surface only
[X,Y,Z] = meshgrid(x(2:end-1),y(2:end-1),z(2:end-1));
p1vec = (-1).^(L + M + N)./( (2*L + 1).*(2*M + 1).*(2*N + 1) );
alpha = (pi^2*D*tau_l/(R_l*4))*( ((2*L + 1)/a).^2 + ...
      ((2*M + 1)/b).^2 + ((2*N + 1)/c).^2 );
%Generate terms to sum, non-surface only (minus the two end-points in xyz)
for i = 1:length(x)-2        
    for j = 1:length(y)-2
        for k = 1:length(z)-2
            for el = 1:length(t)
                %Compute terms for conc calc
                p2vec{i,j,k} = cos( (2*L + 1)*pi*x(i+1) / (2*a) ).*...
                                      cos( (2*M + 1)*pi*y(j+1) / (2*b) ).*...
                                      cos ( (2*N + 1)*pi*z(k+1) / (2*c) );                
                p3vec{el} = exp(-t(el)*alpha);    
                terms{i,j,k,el} = p1vec.*p2vec{i,j,k}.*p3vec{el};                               
            end
        end
    end    
end
%Summation: Cn
for i = 1:(nx-2)*(ny-2)*(nz-2)*length(t)
    sumTS(i) = sum(terms{i},'all');    
end
Cn = (C0 - 64*C0*sumTS/pi^3)*10^3; %mg/l

I think I saved some time by vectorizing the l,m,n summation but the FOR loops at 100 terms (inf_end = 100) still took over 4 minutes to complete. Is there a way for me to compute p2vec faster?

Answer 1

[Stijn Haenen]

Maybe you can compute p2vec{i,j} outside the k-for loop and multiply it inside the for loop with the k-factor.