# Optimize evaluation of triple infinite series using FOR loops

2 views (last 30 days)
Kien Pham on 14 Nov 2019
Answered: Stijn Haenen on 14 Nov 2019
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?

Stijn Haenen on 14 Nov 2019
Maybe you can compute p2vec{i,j} outside the k-for loop and multiply it inside the for loop with the k-factor.

### Categories

Find more on Calculus in Help Center and File Exchange

R2019b

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!