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!