`mmf.math`¶

`lerch_`
`fermi_`
`fermi`
`lerch`
`differentiate`	Differentiation.
`earray`	Extensible numpy arrays.
`fft`	Wrappers for the FFT. Uses the :pkg:`anfft` package if available (which
`integrate`	Integration Utilities.
`interp`	Interpolation library for use with scipy.
`linalg`	Linear Algebra Routines
`multigrid`	Multigrid Solvers.
`ode`
`special`	Special Functions.
`fermi_functions`	Compute various properties of free fermions.
`ind2sub`(shape, ind)	Return the tuple of subscripts indexing into an array of specified shape from a linear index.
`sub2ind`(shape, sub)	Return the linear index corresponding to the tuple of subscripts sub into an array of specified shape.
`almost_equal`(a, b[, absTol, relTol])	Checks if a almost equals b to within the tolerances specified.
`almost_zero`(a[, absTol])	Checks if a is almost zero to within the tolerances specified.
`EllipticF`(theta, m)	Incomplete elliptic integrals of the first kind.
`EllipticE`(theta, m)	Incomplete elliptic integrals of the second kind.
`ellipsoid_area`(a[, b, c])	The surface area of the ellipsoid: x^2/a^2+y^2/b^2+z^2/c^2 = 1.
`ellipse_perimeter`(a[, b])	The perimeter of the ellipse: x^2/a^2+y^2/b^2 = 1.
`cbrt`(c3, c2, c1, c0)	Return the three roots (r1,r2,r3) to the cubic equation c3x^3+c2x^2+c1*x+c0=0.
`step`(x[, x0, dx, d])	Return the `d`th derivative of a smoothed step function.
`LambertW`(x[, k, abs_tol, rel_tol])	Return $w=W_k(x)$ where $x = w e^w$ to the desired precision.
`S`(d)	Return the surface-area of a d-dimensional sphere (without the r**(d-1) factor.
`LerchPhi`(z, s, a)	Return the Lerch Trancident function $\Phi(z,s,a)$ on a restricted
`Li`(s, z)	Return the polylogarithm $\Li_s(z)$ .
`Fi`(j, z)	Return the complete Fermi-Dirac integral $\operatorname{Fi}_{j}(z)$ .
`fermi_integral`(mu, m, T[, pmax, singlet, ...])	Return density of a single species of free fermions of mass m

class mmf.math.fermi_functions¶

Bases: object

Compute various properties of free fermions.

Examples

>>> fermi_functions.P(1, 1, 1)               
array(0.23482812162868...)
>>> fermi_functions.N(1, 1, 1)               
array(0.160318814464720...)
>>> fermi_functions.N_s(1, 1, 1)              
array(0.075970412559876...)
>>> fermi_functions.dN_dmu(1, 1, 1)          
array(0.193444769294...)
>>> fermi_functions.dN_dm(1, 1, 1)           
array(-0.0425632327406...)

__init__()¶: x.__init__(...) initializes x; see help(type(x)) for signature

classmethod N(mu, m, T)¶: Return the density.

classmethod N_s(mu, m, T)¶: Return the scalar density.

classmethod P(mu, m, T)¶: Return the pressure.

classmethod dN_dm(mu, m, T)¶: Return the derivative of the density wrt. m.

classmethod dN_dmu(mu, m, T)¶: Return the derivative of the density wrt. mu.

classmethod get_imt()¶: Return the class IMT instance. This is expensive to compute and so is “memoized”.

mmf.math.ind2sub(shape, ind)¶

Return the tuple of subscripts indexing into an array of specified shape from a linear index.

>>> shape = (2,3,4)
>>> ind = 3
>>> ind - sub2ind(shape,ind2sub(shape,ind))
0
>>> a = np.random.random_sample(shape)
>>> b = a.ravel()
>>> c = np.zeros(len(b))
>>> for n in range(len(b)):
...     c[n] = a[ind2sub(a.shape,n)]
>>> max(abs(c-b))
0.0

mmf.math.sub2ind(shape, sub)¶

Return the linear index corresponding to the tuple of subscripts sub into an array of specified shape.

>>> shape = (2,3,4)
>>> sub = (1,2,1)
>>> a = np.random.random_sample(shape)
>>> b = a.ravel()
>>> a[sub] - b[sub2ind(a.shape,sub)]
0.0
>>> c = np.zeros(len(b))
>>> for n0 in range(shape[0]):
...  for n1 in range(shape[1]):
...   for n2 in range(shape[2]):
...    c[sub2ind(shape,(n0,n1,n2))] = a[n0,n1,n2]
>>> max(abs(c-b))
0.0

mmf.math.almost_equal(a, b, absTol=2.2204460492503131e-16, relTol=2.2204460492503131e-16)¶

Checks if a almost equals b to within the tolerances specified.

>>> almost_equal(0,0)
True
>>> almost_equal(0.1,0)
False
>>> almost_equal(1.0,1.02,absTol=0.1,relTol=0.1)
True
>>> almost_equal(1.0,1.2,absTol=0.1,relTol=0.1)
False

mmf.math.almost_zero(a, absTol=2.2204460492503131e-16)¶

Checks if a is almost zero to within the tolerances specified.

>>> almost_zero(0.02,absTol=0.1)
True
>>> almost_zero(0.02,absTol=0.01)
False

mmf.math.EllipticF(theta, m)¶

Incomplete elliptic integrals of the first kind.

EllipticF(theta,m) = int_0^theta 1/sqrt{1-msin^2(t)} d{t}

mmf.math.EllipticE(theta, m)¶

Incomplete elliptic integrals of the second kind.

EllipticE(theta,m) = int_0^theta sqrt{1-msin^2(t)} d{t}

mmf.math.ellipsoid_area(a, b=None, c=None)¶

The surface area of the ellipsoid:

x^2/a^2+y^2/b^2+z^2/c^2 = 1.

ellipsoid_area(a) == ellipsoid_area(a,a,a) ellipsoid_area(a,b) == ellipsoid_area(a,b,b)

mmf.math.ellipse_perimeter(a, b=None)¶

The perimeter of the ellipse:

x^2/a^2+y^2/b^2 = 1.

ellipse_perimeter(a) == ellipse_perimeter(a,a)

mmf.math.cbrt(c3, c2, c1, c0)¶: Return the three roots (r1,r2,r3) to the cubic equation c3*x^3+c2*x^2+c1*x+c0=0. If the coefficients are real, r1 is guarenteed to be real.

mmf.math.step(x, x0=0.0, dx=0.0, d=0)¶

Return the `d`th derivative of a smoothed step function.

This is based on the following function

$\frac{1 + \tanh\left(\frac{\pi}{2}\tan\frac{\pi x}{2}\right)}{2}$

Examples

>>> import numpy as np
>>> step(-1.0)
0.0
>>> step(0.0)
0.5
>>> step(1.0, 1e-6, 1e-6/2)
1.0
>>> step(np.array(1.0), 1e-6, 1e-6/2)
1.0

The step function is implemented so that it is C(inf) smooth if dx > 0, strictly 1.0 for x > x0 + dx and strictly 0.0 for x < x0 - dx.

>>> step(np.array([0.9 - 1e-15,1.0,1.1 + 1e-15]),1.0,0.1)
array([ 0. ,  0.5,  1. ])

mmf.math.LambertW(x, k=0, abs_tol=2.2204460492503131e-16, rel_tol=2.2204460492503131e-16)¶

Return $w=W_k(x)$ where $x = w e^w$ to the desired precision.

Parameters :

x : float

k : -1, 1 or 0

Used to select the branch for x < 0.

Notes

Proceeds using the iteration:

$w \mapsto w - \frac{w e^w - x} {(w+1)e^w - \frac{(w+2)(w e^w - x)}{2w + 2}}$

which is based on Halley’s method. Close to the solution, the absolute error is approximately:

$w - w^* \approx \frac{x - w e^w}{(1+w)e^w}.$

The branch $k$ is chosen via the initial guess. This is done with the first two terms in the asymptotic approximation for the non-principle branches

$W_{k}(x) \approx \ln_k{x} - \ln(\ln_k{x}$

where $\ln_k{x} = \ln{x} + 2\pi i k$ and $\ln{x}$ is the principle valued function. This does not work for the branches $k=\pm 1$ and the principle branch $k=0$ near $0$ and $-e^{-1}$ . Near the branch point $-e^{-1}$ we use the series approximation

$W(x) \approx -1 - p - \tfrac{1}{3} p^2 + \tfrac{11}{72} p^3 + \cdots$

where $p = \pm\sqrt{2ex + 1}$ . If $\Im{x} >=0$ , then the positive $p$ is a good approximation of the $k=-1$ branch, otherwise the positive $p$ gives a good approximation for the $k=1$ branch.

For the principle branch, we use the approximation due to Winitzki??:

$W(x) \approx \frac{2\ln(1+Bp) - \ln[1 + C\ln(1+Dp)] +E} {1+(2\ln(1+Bp) + A)^-1}$

where $A=4.6948$ , $B=0.8842$ , $C=0.9294$ , $D=0.5106$ , and $E=-1.213$ .

Examples

>>> LambertW(-0.2, k=0)             
(-0.2591711018...
>>> LambertW(-0.2, k=1)             
(-3.7223204849...+7.38723021...j)
>>> LambertW(-0.2, k=-1)            
(-2.54264135777...
>>> LambertW(0.0001, k=0)           
(9.9990001...
>>> LambertW(1.0,0)                 
(0.5671432904...
>>> LambertW(1.0,1)                 
(-1.533913319...+4.375185153...j)
>>> LambertW(1.0,-1)                 
(-1.533913319...-4.375185153...j)
>>> LambertW(1.0,2)                 
(-2.4015851048...+10.776299516...j)
>>> LambertW(-10.0,0)               
(1.3699809685...+2.140194527...j)

Todo

Fix initial guess near 0. This does not work for k!=0 branches:

>>> W = LambertW(0.001,k=-1)        
>>> np.allclose(W*np.exp(W), 0.001) 
True

mmf.math.S(d)¶

Return the surface-area of a d-dimensional sphere (without the r**(d-1) factor.

Examples

>>> S(3)
12.566370614359...
>>> S(2)
6.2831853071795...
>>> S(1)
2.0...

mmf.math.LerchPhi(z, s, a)¶

Return the Lerch Trancident function $\Phi(z,s,a)$ on a restricted domain.

$\Phi(z,s,a) = \sum_{n=0}^\infty \frac{z^n}{(n+a)^s}$

Computes the result using the integral representation

$\Phi(z,s,a) = \frac{1}{\Gamma(s)}\int_0^\infty \d{t}\; \frac{t^{s-1}e^{-at}}{1-ze^{-t}}$

Examples

>>> abs(LerchPhi(1,2,3) - (np.pi**2/6.0 - 5.0/4.0)) < 1e-12
True
>>> abs(LerchPhi(0.5,2,3) - 0.157924211720100047221250) < 1e-12
True

mmf.math.Li(s, z)¶: Return the polylogarithm $\Li_s(z)$ .

$\Li_s(z) &= \sum_{k=1}^{\infty} \frac{z^k}{k^s} = z\Phi(z, s, 1),$

mmf.math.Fi(j, z)¶: Return the complete Fermi-Dirac integral $\operatorname{Fi}_{j}(z)$ .

$\operatorname{Fi}_{j}(z) &= \frac{1}{\Gamma(j + 1)}\int_{0}^{\infty}\d{t}\; \frac{t^{j}}{e^{t - z} + 1} = e^{z}\Phi(-e^{z},j+1,1)$

mmf.math.fermi_integral(mu, m, T, pmax=None, singlet=False, abs_tol=None, rel_tol=None)¶

Return density of a single species of free fermions of mass m at temperature T and chemical potential mu in natural units ( $\hbar = c = 1$ ). Does not include any degeneracy factors:

$\int_{0}^{p_{\text{max}}}\d{p}\; \frac{p^{2}}{1+e^{(\sqrt{p^2 + m^2} - \mu)/T}} - \frac{p^{2}}{1+e^{(\sqrt{p^2 + m^2} + \mu)/T}}$

Parameters :

mu : float

Chemical potential

m : float

Mass

T : float

Temperature

pmax : float

Upper limit on integration.

singlet : bool

If True, then the singlet density will be computed:

$\int_{0}^{p_{\text{max}}}\d{p}\; \frac{m}{\sqrt{p^2 + m^2}} \left( \frac{p^{2}}{1+e^{(\sqrt{p^2 + m^2} - \mu)/T}} - \frac{p^{2}}{1+e^{(\sqrt{p^2 + m^2} + \mu)/T}} \right)$

Examples

>>> import math, numpy as np
>>> mu = 3.0
>>> m = 2.0
>>> T = 0.0
>>> pF = math.sqrt(mu*mu-m*m)
>>> N = fermi_integral(mu, m, T)
>>> np.allclose(N, pF**3/6/pi**2)
True

>>> import math, numpy as np
>>> mu = 30.0
>>> m = 22.0
>>> T = 0.00001
>>> pF = math.sqrt(mu*mu-m*m)
>>> N = fermi_integral(mu, m, T)
>>> np.allclose(N, pF**3/6/pi**2)
True

Next topic

This Page

`mmf.math`¶

Navigation

Next topic

This Page

Quick search

mmf.math¶

Navigation

`mmf.math`¶