How to quickly calculate the following function?

4 views (last 30 days)
Hi all,I want to calculate the function shown below:
This is Lerch transcendent function( https://en.wikipedia.org/wiki/Lerch_zeta_function ).
I use the following formula to calculate it, but this is too slow (I need to calculate it many times):
There is also no official function here to calculate it.
So I would like to know how to calculate it quickly, thank you all in advance!
  5 Comments

Sign in to comment.

Accepted Answer

Jan
Jan on 4 Sep 2022
Calculating the sum seems to be more stable. For cases, in which the integral method is successful, the sum is 10 to 200 times faster for the given test data:
function test_Lerch
% See: https://people.math.sc.edu/Burkardt/py_src/polpak/lerch_values.py
% and: https://people.sc.fsu.edu/~jburkardt/m_src/test_values/test_values.html
a_vec = [0.0E+00, ...
0.0E+00, ...
0.0E+00, ...
1.0E+00, ...
1.0E+00, ...
1.0E+00, ...
2.0E+00, ...
2.0E+00, ...
2.0E+00, ...
3.0E+00, ...
3.0E+00, ...
3.0E+00];
s_vec = [2, 3, 10, ...
2, 3, 10, ...
2, 3, 10, ...
2, 3, 10];
z_vec = [0.1000000000000000E+01, ...
0.1000000000000000E+01, ...
0.1000000000000000E+01, ...
0.5000000000000000E+00, ...
0.5000000000000000E+00, ...
0.5000000000000000E+00, ...
0.3333333333333333E+00, ...
0.3333333333333333E+00, ...
0.3333333333333333E+00, ...
0.1000000000000000E+00, ...
0.1000000000000000E+00, ...
0.1000000000000000E+00];
f_vec = [ ...
0.1644934066848226E+01, ...
0.1202056903159594E+01, ...
0.1000994575127818E+01, ...
0.1164481052930025E+01, ...
0.1074426387216080E+01, ...
0.1000492641212014E+01, ...
0.2959190697935714E+00, ...
0.1394507503935608E+00, ...
0.9823175058446061E-03, ...
0.1177910993911311E+00, ...
0.3868447922298962E-01, ...
0.1703149614186634E-04];
for k = 1:numel(z_vec)
z = z_vec(k);
s = s_vec(k);
a = a_vec(k);
f = f_vec(k);
fprintf('\nk: %d\n', k)
tic
for k = 1:1
y1 = Lerch_fun_integral(z,s,a);
end
toc
tic
for k = 1:1
y2 = Lerch_fun_sum(z,s,a);
end
toc
fprintf('Ref: %.16g\n', f);
fprintf('Int: %.16g delta: %16g\n', y1, abs(y1 - f));
fprintf('Sum: %.16g delta: %16g\n', y2, abs(y2 - f));
end
end
function out = Lerch_fun_integral(z,s,a)
f=@(t)((t.^(s-1)).*exp(-a*t))./(1-z.*exp(-t));
out=(1/gamma(s))*(integral(f,0,Inf));
end
function out = Lerch_fun_sum(z, s, a)
out = 0;
if z <= 0
return
end
lim = 1e-16;
k = 0;
z_k = 1;
term = Inf;
while abs(term) > lim * (1.0 + abs(out))
if a + k ~= 0
term = z_k / (a + k)^s;
out = out + term;
end
k = k + 1;
z_k = z_k * z;
end
end

More Answers (0)

Products


Release

R2018b

Community Treasure Hunt

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

Start Hunting!