Sie sind hier: ICP » R. Hilfer » Publikationen

2 Method of Calculation

[page 639, §1]   

2.1 Recursion Relation

[639.1.1] In the rest of this paper α and β are real numbers with α>0. [639.1.2] Calculating Eα,β⁢z for the case α>1 can be reduced to the case 0<α≤1 by virtue of the recursion relation [2, 21]

Eα,β⁢z=12⁢m+1⁢∑h=-mmEα/2⁢m+1,β⁢z1/2⁢m+1⁢e2⁢π⁢i⁢h/2⁢m+1(6)

valid for all β∈R, z∈C where m=α-1/2+1. [639.1.3] Here x denotes the largest integer smaller than x. [639.1.4] Because of the recursion we consider only the case 0<α≤1 in the following.

2.2 Regions

[639.2.1] For 0<α≤1 the complex z-plane is partitioned into four regions G0, G1, G2, G3. [639.2.2] In each region a different method is used for the calculation of the Mittag-Leffler function. [639.2.3] The central region

G0=D⁢r0¯=z∈C:z≤r0(7)

is the closure of the open disk D⁢r0=z∈C:z<r0 of radius r0<1 centered at the origin. [639.2.4] In our numerical calculations we have chosen r0=0.9 throughout. [639.2.5] The regions G1,G2 are defined by

G1=W+⁢5⁢π⁢α/6∖G0(8)
G2=C∖G1∪G0(9)

where

W+⁢ϕ=z∈C:arg⁡z<ϕ(10)

is the open wedge with opening angle ϕ,(0<ϕ<π). [639.2.6] For future reference we define also the (left) wedge

W-⁢ϕ=z∈C:ϕ<arg⁡z≤π(11)

[639.2.7] Finally the region G3 is defined as

G3=C∖D⁢r1=z∈C:z>r1(12)

where r1>r0 will be defined below.

2.3 Taylor series

[639.3.1] For z∈G0, i.e. for the central disk region (with r0=0.9), the Mittag-Leffler function is computed from truncating its Taylor series (1). [639.3.2] For given z∈G0 and accuracy ϵ we choose [page 640, §0]    N such that

Eα,β⁢z-∑k=0NzkΓ⁢α⁢k+β≤ϵ(13)

holds. [640.0.1] The dependence of N on the accuracy ϵ and other parameters is found as

N≥max⁡2-β/α+1,ln⁡ϵ⁢1-zln⁡z+1.(14)

2.4 Integral Representations

[640.1.1] We start from the basic integral representation [2, p. 210]

Eα,β⁢z=12⁢π⁢i⁢∫Cyα-β⁢eyyα-z⁢d⁢y(15)

where the path of integration C in the complex plane starts and ends at -∞ and encircles the circular disc y≤z1/α in the positive sense. [640.1.2] When this integral is evaluated several cases arise. [640.1.3] For z∈G1 we distinguish the cases β≤1 and β>1 and compute the Mittag-Leffler function from the integral representations [22]

Eα,β⁢z=A⁢z;α,β,0+∫0∞B⁢r;α,β,z,π⁢α⁢d⁢r,β≤1(16)
Eα,β⁢z=A⁢z;α,β,0+∫1/2∞B⁢r;α,β,z,π⁢α⁢d⁢r+∫-π⁢απ⁢αC⁢φ;α,β,z,1/2⁢d⁢φ,β>1(17)

where

A⁢z;α,β,x=1α⁢z1-β/α⁢exp⁡z1/α⁢cos⁡x/α(18)
B⁢r;α,β,z,ϕ=1π⁢A⁢r;α,β,ϕ⁢r⁢sin⁡ω⁢r,ϕ,α,β-ϕ-z⁢sin⁡ω⁢r,ϕ,α,βr2-2⁢r⁢z⁢cos⁡ϕ+z2(19)
C⁢φ;α,β,z,ρ=ρ2⁢π⁢A⁢ρ;α,β,φ⁢cos⁡ω⁢ρ,φ,α,β+i⁢sin⁡ω⁢ρ,φ,α,βρ⁢cos⁡φ+i⁢sin⁡φ-z(20)
ω⁢x,y,α,β=x1/α⁢sin⁡y/α+y⁢1+1-β/α.(21)

[640.1.4] For z∈G2 the integral representations read

Eα,β⁢z=∫0∞B⁢r;α,β,z,2⁢π⁢α/3⁢d⁢r,β≤1(22)
Eα,β⁢z=∫1/2∞B⁢r;α,β,z,2⁢π⁢α/3⁢d⁢r+∫-2⁢π⁢α/32⁢π⁢α/3C⁢φ;α,β,z,1/2⁢d⁢φ,β>1(23)

where the integrands have been defined in eq. (18)–(21) above.

[page 641, §1]   
[641.1.1] The integrand C⁢φ;α,β,z,ρ is oscillatory but bounded over the integration interval. [641.1.2] Thus the integrals over C can be evaluated numerically using any appropriate quadrature formula. [641.1.3] We use a robust Gauss-Lobatto quadrature.

[641.2.1] The integrals over B⁢r;α,β,z,ϕ involve infinite intervals. [641.2.2] For given z∈G1 (resp. z∈G2) and accuracy ϵ we approximate the integrals by truncation. [641.2.3] The error

∫Rmax∞B⁢r;α,β,z,ϕ⁢d⁢r<ϵ(24)

depends on the truncation point Rmax. [641.2.4] For ϕ=π⁢α we find

Rmax≥max⁡1, 2⁢z,-ln⁡π⁢ϵ6α,β≥0max⁡β+1α, 2⁢z,-2⁢ln⁡π⁢ϵ6⁢β+2⁢2⁢ββα,β<0,(25)

while for ϕ=2⁢π⁢α/3 we have

Rmax≥max⁡2α, 2⁢z,-2⁢ln⁡π⁢2β⁢ϵ12α,β≥0max⁡2⁢β+1α, 2⁢z,-4⁢ln⁡π⁢2β⁢ϵ12⁢β+2⁢4⁢ββα,β<0.(26)

2.5 Exponential Asymptotics

[641.3.1] The asymptotic expansions given in eqs. (21) and (22) on page 210 in Ref. [2] indicate that the Mittag-Leffler function exhibits a Stokes phenomenon. [641.3.2] The Stokes lines are the rays arg⁡z=±α⁢π and it was recently shown in [23] that a Berry-type smoothing applies. [641.3.3] The exponentially improved asymptotic expansions will be used here for computing Eα,β⁢z for z∈G3. [641.3.4] More precisely for 0<α<1 we use the following exponentially improved uniform asymptotic expansions [23] :

Eα,β⁢z=-∑k=1Mz-kΓ⁢β-α⁢k+O⁢z1/2-β/α⁢e-z1/α(27)

for z∈G3∩W-⁢α⁢π+δ

Eα,β⁢z=1α⁢z1-β/α⁢ez1/α-∑k=1Mz-kΓ⁢β-α⁢k+O⁢z1/2-β/α⁢e-z1/α(28)

for z∈G3∩W+⁢α⁢π-δ,

Eα,β⁢z=12⁢α⁢z1-β/α⁢ez1/α⁢erfc⁡c3/36-i⁢c2/6-c⁢z1/α/2
-∑k=1Mz-kΓ⁢β-α⁢k+O⁢z1/2-β/α⁢e-z1/α(29)

[page 642, §0]    for z∈G3 with -α⁢π-δ≤arg⁡z≤-α⁢π+δ and c=arg⁡z1/α+π, and

Eα,β⁢z=12⁢α⁢z1-β/α⁢ez1/α⁢erfc⁡c+i⁢c2/6-c3/36⁢z1/α/2
-∑k=1Mz-kΓ⁢β-α⁢k+O⁢z1/2-β/α⁢e-z1/α(30)

for z∈G3 with α⁢π-δ≤arg⁡z≤α⁢π+δ and c=arg⁡z1/α-π. [642.0.1] Here δ>0 is a small number, in practice we use δ=0.01. [642.0.2] In eqs. (2.5) and (2.5)

erfc⁡z=2π⁢∫z∞ey2/2⁢d⁢y(31)

denotes the complex complementary error function. [642.0.3] Note the difference between eqs. (27)–(2.5) and the expansions in [2, p. 210]. [642.0.4] We take

M=min⁡50,z1/α/α(32)

to truncate the asymptotic series. [642.0.5] Choosing for ϵ the machine precision we fix the radius r1 in G3 at

r1=1ϵ⁢Γ⁢β-α⁢M1/M.(33)

[642.0.6] The remainder estimates in eqs. (27)–(2.5) are only valid for a value of M that is chosen in an optimal sense as in eqs. (32) and (33).