# Optimize evaluation of triple infinite series using FOR loops

2 views (last 30 days)

Show older comments

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?

##### 0 Comments

### Answers (1)

Stijn Haenen
on 14 Nov 2019

##### 0 Comments

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!