14
1 Gamma Functions, Beta Functions, and Related
Ci(z, a) ± i Si(z, a) = exp
±
1
2
iπa
γ (∓i z, a).
(1.35)
From Eq. (1.35) we obtain easily
Ci(z, a) =
exp
1
2
iπa
γ (−i z, a) + exp
−
1
2
iπa
γ (i z, a)
/2 and
(1.36)
Si(z, a) =
exp
1
2
iπa
γ (−i z, a) − exp
−
1
2
iπa
γ (i z, a)
/2. (1.37)
The MATLAB-function q = integral(fun,xmin,xmax) (for details
see the MATLAB-documentation) allows directly to integrate Eqs. (1.32), and (1.33),
with “fun” the function handle representing the integrand and “xmin, xmax” the
integral limits. integral is based on an adaptive quadrature.
%% Example: Computing the Boehmer integrals
%
via incomplete gamma function
a0 = exp(i * rand) % a0 = 0.6861 + 0.7275i
z = exp(i * rand) % z = 0.6171 + 0.7869i
%% For sufficiently large a --> no problem
format long
a = 2 * a0;
% direct integration
Cih = @(t) t.^(a-1) . * cos(t);
Ci = integral(Cih,0,z)
% 0.153560594714140 + 0.036271725299627i
%
Sih = @(t) t.^(a-1) . * sin(t);
Si = integral(Sih,0,z)
% 0.006099267216393 + 0.100534030211523i
% incomplete gamma function
[objp,ergp] = incgammaC(i * z,a);
[objm,ergm] = incgammaC(-i * z,a);
Cig = (exp(-1/2 * pi * a * i). * ergp +
exp(1/2 * pi * a * i). * ergm)/2
% 0.153560588891605 + 0.036271722524393i
Sig = (exp(1/2 * pi * a * i). * ergm -
exp(-1/2 * pi * a * i). * ergp)/(2 * i)
% 0.006099267002845 + 0.100534030518022i
The results based on the MATLAB-function integral and the results computed
with the incomplete gamma function are in very good agreement.
%% If a is too small --> insufficient accuracy
a= a0/10
%
Cih = @(t) t.^(a-1) . * cos(t);
Ci = integral(Cih,0,z)
%Warning: Reached the limit on the maximum number of