"""Differentiation."""
from __future__ import division, with_statement
from math import sin, cos
import itertools
import scipy.sparse
import scipy as sp
import numpy as np
#import pygsl.qrng
import mmf.objects
from mmf.math.integrate import Richardson
__all__ = ['differentiate', 'D', 'D1', 'D2', 'D_Hermitian']
[docs]def differentiate(f, x=0.0, d=1, h0=1.0,
l=1.4, nmax=10, dir=0,
p0=1, err=[0]):
r"""Evaluate the numerical dth derivative of f(x) using a Richardson
extrapolation of the finite difference formula.
Parameters
----------
f : function
The function to compute the derivative of.
x : {float, array}
The derivative is computed at this point (or at these points if
the function is vectorized.
d : int, optional
Order of derivative to compute. `d=0` is the function `f(x)`,
`d=1` is the first derivative etc.
h0 : float, optional
Initial stepsize. Should be on about a factor of 10 smaller
than the typical scale over which `f(x)` varies significantly.
l : float, optional
Richardson extrapolation factor. Stepsizes used are `h0/l**n`
nmax : int, optional
Maximum number of extrapolation steps to take.
dir : int, optional
If `dir < 0`, then the function is only evaluated to the
left, if positive, then only to the right, and if zero, then
centered form is used.
Returns
-------
df : {float, array}
Order `d` derivative of `f` at `x`.
Other Parameters
----------------
p0 : int, optional
This is the first non-zero term in the taylor expansion of either the
difference formula. If you know that the first term is zero (because of
the coefficient), then you should set `p0=2` to skip that term.
.. note:: This is not the power of the term, but the order. For centered
difference formulae, `p0=1` is the `h**2` term, which would vanish if
third derivative vanished at `x` while for the forward difference
formulae this is the `h` term which is absent if the second derivative
vanishes.
err[0] : float
This is mutated to provide an error estimate.
Notes
-----
This implementation uses the Richardson extrapolation to
extrapolate the answer. This is based on the following Taylor
series error formulae:
.. math::
f'(x) &= \frac{f(x+h) - f(x)}{h} - h \frac{f'')}{2} + \cdots\\
&= \frac{f(x+h) - f(x-h)}{2h} - h^2 \frac{f''}{3!} + \cdots\\
f''(x) &= \frac{f(x+2h) - 2f(x+h) + f(x)}{h^2} - hf^{(3)} + \cdots\\
&= \frac{f(x+h) -2f(x) + f(x-h)}{h^2}
- 2h^2 \frac{f^{(4)}}{4!} + \cdots\\
If we let $h = 1/N$ then these formula match the expected error
model for the Richardson extrapolation
.. math::
S(h) = S(0) + ah^{p} + \order(h^{p+1})
with $p=1$ for the one-sided formulae and $p=2$ for the centered
difference formula respectively.
See :class:`mmf.math.integrate.Richardson`
See Also
--------
:func:`mmf.math.integrate.Richardson` : Richardson extrapolation
:func:`Richardson` : Richardson extrapolation
Examples
--------
>>> x = 100.0
>>> assert(abs(differentiate(sin, x, d=0)-sin(x))<1e-15)
>>> assert(abs(differentiate(sin, x, d=1)-cos(x))<3e-15)
>>> assert(abs(differentiate(sin, x, d=2)+sin(x))<4e-14)
>>> assert(abs(differentiate(sin, x, d=3)+cos(x))<4e-12)
>>> assert(abs(differentiate(sin, x, d=4)-sin(x))<2e-10)
>>> differentiate(abs, 0.0, d=1, dir=1)
1.0
>>> differentiate(abs, 0.0, d=1, dir=-1)
-1.0
>>> differentiate(abs, 0.0, d=1, dir=0)
0.0
Note that the Richardson extrapolation assumes that `h0` is small
enough that the truncation errors are controlled by the taylor
series and that the taylor series properly describes the behaviour
of the error. For example, the following will not converge well,
even though the derivative is well defined:
>>> def f(x):
... return np.sign(x)*abs(x)**(1.5)
>>> df = differentiate(f, 0.0)
>>> abs(df) < 0.1
True
>>> abs(df) < 0.01
False
>>> abs(differentiate(f, 0.0, nmax=100)) < 3e-8
True
Sometimes, one may compensate by increasing nmax. (One could in
principle change the Richardson parameter p, but this is optimized
for analytic functions.)
The :func:`differentiate` also works over arrays if the function
`f` is vectorized:
>>> x = np.linspace(0, 100, 10)
>>> assert(max(abs(differentiate(np.sin, x, d=1) - np.cos(x))) < 3e-15)
"""
if 0 == d:
return f(x)
elif 1 == d:
if dir < 0:
def df(N, f=f, x=x, h0=h0):
h = float(h0)/N
h = (x + h) - x
return (f(x) - f(x-h))/h
elif dir > 0:
def df(N, f=f, x=x, h0=h0):
h = float(h0)/N
h = (x + h) - x
return (f(x + h) - f(x))/h
else:
def df(N, f=f, x=x, h0=h0):
h = float(h0)/N
h = (x + h) - x
return (f(x + h) - f(x - h))/(2*h)
elif 2 == d:
if dir < 0:
def df(N, f=f, x=x, h0=h0):
h = float(h0)/N
h = (x + h) - x
return (f(x - 2*h) - 2*f(x - h) + f(x))/(h*h)
elif dir > 0:
def df(N, f=f, x=x, h0=h0):
h = float(h0)/N
h = (x + h) - x
return (f(x + 2*h) - 2*f(x + h) + f(x))/(h*h)
else:
def df(N, f=f, x=x, h0=h0):
h = float(h0)/N
h = (x + h) - x
return (f(x + h) - 2*f(x) + f(x - h))/(h*h)
else:
df = lambda x:differentiate(f=f, x=x, d=d-2, h0=h0, dir=dir,
l=l, nmax=nmax)
return differentiate(df, x=x, d=2, h0=h0, dir=dir,
l=l, nmax=nmax, err=err)
if dir == 0:
p = 2
else:
p = 1
r = Richardson(df, ps=itertools.count(p*p0,p), l=l)
r.next()
d1 = r.next()
d0 = r.next()
err_old = abs(d1 - d0)
n = 2
for d in r:
n += 1
err[0] = abs(d - d0)
d0 = d
if (err[0] > err_old).any() or (n > nmax):
break
err_old = err[0]
return r.next()
[docs]def D1(x, sparse=False):
r"""Return D where D is the finite difference
operator over the lattice x.
Parameters
----------
type : 'left' | 'center' | 'right'
Which type of difference formula to use for interior
points.
Examples
--------
>>> x = np.linspace(0, 10, 100)
>>> d1 = D1(x)
>>> f = np.sin(x)
>>> abs(np.dot(d1, f)-np.cos(x)).max() < 0.002
True
Use sparse matrices if there is lots of data.
>>> d1s = D1(x,sparse=True)
>>> f = np.sin(x)
>>> abs(np.dot(d1, f)-np.cos(x)).max() < 0.002
True
Notes
-----
.. code-block:: none
restart;
with(CodeGeneration);
eq:={
f1=f+df*d1+ddf*d1^2/2,
f2=f+df*d2+ddf*d2^2/2};
DF:=simplify(subs(solve(eq, [df, ddf])[1], df));
simplify(coeff(DF, f));
simplify(coeff(DF, f1));
simplify(coeff(DF, f2));
eq:={
f1=f+df*d1+ddf*d1^2/2+dddf*d1^3/6,
f2=f+df*d2+ddf*d2^2/2+dddf*d2^3/6,
f3=f+df*d3+ddf*d3^2/2+dddf*d3^3/6};
DF:=simplify(subs(solve(eq, [df, ddf, dddf])[1], df));
Matlab(simplify(coeff(DF, f)));
Matlab(simplify(coeff(DF, f1)));
Matlab(simplify(coeff(DF, f2)));
Matlab(simplify(coeff(DF, f3)));
"""
inds = np.arange(len(x))
inds_1 = np.array([1] + range(0, len(x)-2) + [len(x)-2])
inds_2 = np.array([2] + range(2, len(x)) + [len(x)-3])
inds_3 = np.array([3] + range(2, len(x)) + [len(x)-4])
x = np.array(x)
d1 = x[inds_1] - x[inds]
d2 = x[inds_2] - x[inds]
if not sparse:
D = np.zeros((len(x),)*2, dtype=float)
D[inds, inds] = -((d1 + d2)/(d1*d2))
D[inds, inds_1] = (d2/d1/(d2-d1))
D[inds, inds_2] = (d1/d2/(d1-d2))
else:
data = np.hstack([
-((d1 + d2)/(d1*d2)),
(d2/d1/(d2-d1)),
(d1/d2/(d1-d2))])
rows = np.hstack([inds, inds, inds])
cols = np.hstack([inds, inds_1, inds_2])
D = sp.sparse.coo_matrix((data, (rows, cols)),
shape=(len(x),)*2, dtype=float)
D = D.tolil()
# Use higher order formula for endpoints
with np.errstate(divide='ignore'):
if 3 < len(x):
d3 = x[inds_3] - x[inds]
for i in [0, len(x)-1]:
D[i, inds[i]] = (-(d1*d2+d1*d3+d3*d2)/d2/d1/d3)[i]
D[i, inds_1[i]] = (d3*d2/(-d1*d2+d3*d2+d1*d1-d1*d3)/d1)[i]
D[i, inds_2[i]] = (-d1*d3/(-d1*d3+d1*d2-d2*d2+d3*d2)/d2)[i]
D[i, inds_3[i]] = (d2*d1/(-d1*d3+d1*d2+d3*d3-d3*d2)/d3)[i]
if sparse:
D = D.tocsr()
else:
D = np.mat(D)
return D
[docs]def D2(x):
r"""Return D2 where D2 is the finite difference
operator for the second derivative over the lattice x.
Examples
--------
>>> x = np.linspace(0, 10, 100)
>>> d1 = D2(x)
>>> f = np.sin(x)
>>> abs(np.dot(d1, f)+np.sin(x)).max() < 0.001
True
>>> x = np.cumsum(np.random.rand(100))
>>> x = x/x.max()*6.0
>>> d1 = D2(x)
>>> f = np.sin(x)
>>> abs(np.dot(d1, f)+np.sin(x)).max() < 0.05
True
Notes
-----
.. code-block:: none
restart;
with(CodeGeneration);
eq:={
f1=f+df*d1+ddf*d1^2/2,
f2=f+df*d2+ddf*d2^2/2};
DF:=simplify(subs(solve(eq, [df, ddf])[1], ddf));
Matlab(simplify(coeff(DF, f)));
Matlab(simplify(coeff(DF, f1)));
Matlab(simplify(coeff(DF, f2)));
eq:={
f1=f+df*d1+ddf*d1^2/2+dddf*d1^3/6,
f2=f+df*d2+ddf*d2^2/2+dddf*d2^3/6,
f3=f+df*d3+ddf*d3^2/2+dddf*d3^3/6};
DF:=simplify(subs(solve(eq, [df, ddf, dddf])[1], ddf));
Matlab(simplify(coeff(DF, f)));
Matlab(simplify(coeff(DF, f1)));
Matlab(simplify(coeff(DF, f2)));
Matlab(simplify(coeff(DF, f3)));
eq:={
f1=f+df*d1+ddf*d1^2/2+dddf*d1^3/6+ddddf*d1^4/24,
f2=f+df*d2+ddf*d2^2/2+dddf*d2^3/6+ddddf*d2^4/24,
f3=f+df*d3+ddf*d3^2/2+dddf*d3^3/6+ddddf*d3^4/24,
f4=f+df*d4+ddf*d4^2/2+dddf*d4^3/6+ddddf*d4^4/24};
DF:=simplify(subs(solve(eq, [df, ddf, dddf, ddddf])[1], ddf));
Matlab(simplify(coeff(DF, f)));
Matlab(simplify(coeff(DF, f1)));
Matlab(simplify(coeff(DF, f2)));
Matlab(simplify(coeff(DF, f3)));
Matlab(simplify(coeff(DF, f4)));
"""
inds = np.arange(len(x))
inds_1 = np.array([1] + range(0, len(x)-2) + [len(x)-2])
inds_2 = np.array([2] + range(2, len(x)) + [len(x)-3])
inds_3 = np.array([3] + range(2, len(x)) + [len(x)-4])
inds_4 = np.array([4] + range(2, len(x)) + [len(x)-5])
x = np.array(x)
d1 = x[inds_1] - x[inds]
d2 = x[inds_2] - x[inds]
D = np.zeros((len(x), len(x)), dtype=float)
for i in xrange(len(x)):
D[i, inds[i]] = (2.0/d1/d2)[i]
D[i, inds_1[i]] = (2.0/d1/(d1-d2))[i]
D[i, inds_2[i]] = (2.0/d2/(d2-d1))[i]
# Use higher order formula for endpoints
if 4 == len(x):
d3 = x[inds_3] - x[inds]
for i in [0, len(x)-1]:
D[i, inds[i]] = (2.0*(d1+d2+d3)/d1/d2/d3)[i]
D[i, inds_1[i]] = (-2.0*(d2+d3)/(d2*d3-d1*d2+d1*d1-d1*d3)/d1)[i]
D[i, inds_2[i]] = (2.0*(d1+d3)/(d1*d2-d1*d3-d2*d2+d2*d3)/d2)[i]
D[i, inds_3[i]] = (-2.0*(d1+d2)/(d1*d2-d1*d3-d2*d3+d3*d3)/d3)[i]
if 4 < len(x):
d3 = x[inds_3] - x[inds]
d4 = x[inds_4] - x[inds]
with np.errstate(divide='ignore'):
for i in [0, len(x)-1]:
D[i, inds[i]] = (2*(d1*d2+d1*d3+d1*d4+d2*d3+d2*d4+d3*d4)/\
d1/d3/d2/d4)[i]
D[i, inds_1[i]] = (2*(d2*d3+d2*d4+d3*d4)/\
(d2*d3*d1-d4*d2*d3-d2*d1*d1\
+d2*d4*d1-d3*d1*d1+d3*d4*d1\
+d1*d1*d1-d4*d1*d1)/d1)[i]
D[i, inds_2[i]] = (-2*(d3*d1+d4*d1+d3*d4)/\
(-d2*d3*d1+d3*d4*d1+d1*d2*d2\
-d2*d4*d1+d2*d2*d3-d4*d2*d3\
-d2*d2*d2+d4*d2*d2)/d2)[i]
D[i, inds_3[i]] = (2*(d1*d2+d4*d1+d2*d4)/\
(-d2*d4*d1+d2*d3*d1-d1*d3*d3+d3*d4*d1\
-d3*d3*d2+d4*d2*d3+d3*d3*d3\
-d4*d3*d3)/d3)[i]
D[i, inds_4[i]] = (-2*(d1*d2+d3*d1+d2*d3)/\
(-d2*d4*d1+d2*d3*d1+d1*d4**2-d3*d4*d1\
+d2*d4*d4-d4*d2*d3-d4*d4*d4\
+d3*d4*d4)/d4)[i]
return np.mat(D)
[docs]def D(x, n=2):
r"""Return D where D is the nth order finite difference
operator over the lattice x.
Parameters
----------
n : int
Order. Note only n=1, 2 are reliable.
Examples
--------
>>> x = np.cumsum(np.sin(np.linspace(0, 100, 500))**2)
>>> x = x/abs(x).max()
>>> d1 = D(x, 1)
>>> d2 = D(x, 2)
>>> f = np.sin(x)
>>> abs(np.dot(d1, f)-np.cos(x)).max() < 1e-5
True
>>> abs(np.dot(d2, f)+np.sin(x)).max() < 0.001
True
"""
if 0 == n:
return np.eye(len(x))
elif 1 == n:
return D1(x)
elif 2 == n:
return D2(x)
elif 2 < n:
d2 = np.mat(D2(x))
def f(n):
if n == 1:
return D1(x)
elif n == 2:
return d2
else:
return d2*f(n-2)
return f(n)
else:
raise ValueError("0 <= n must be an integer.")
[docs]def D_Hermitian(x, n=2, boundary='zero'):
r"""Return `(D, g, ginv)` where D is the nth order finite difference
operator over the lattice x with boundary conditions:
Parameters
----------
x : array of floats
Abscissa
n : int
Order of derivative
boundary : 'zero'
Specify boundary conditions:
`'zero'`: f(x[0]) = f(x[-1]) = 0
"""
if n != 2:
raise ValueError("Only order n=2 supported.")
if boundary is not 'zero':
raise ValueError("Only boundary conditions 'zero' supported.")
if x[0] == -np.inf:
x[0] = np.finfo(float).min
if x[-1] == np.inf:
x[-1] = np.finfo(float).max
# Compute finite diferences etc.
dx = x[1:]-x[0:-1]
d2x = x[2:]-x[0:-2]
# Form matrices
a = -2.0/dx[1:]/dx[0:-1]
b = 2.0/dx[1:-1]/d2x[:-1]
c = 2.0/dx[1:-1]/d2x[1:]
D2 = np.matrix(np.diag(a, 0)+np.diag(c, -1)+np.diag(b, 1))
g = np.matrix(np.diag(1./np.sqrt(d2x)))
ginv = np.matrix(np.diag(np.sqrt(d2x)))
return D2, g, ginv
def Hermitize(T):
r"""Return (Ginv*T*G, G, Ginv) where G is diagonal and G*T*Ginv is
Hermitian.
"""
(N, M) = T.shape
assert M==N
a = np.diag(T, 0) # diagonal
b = np.diag(T, 1) # super diagonal
c = np.diag(T, -1) # sub diagonal
assert (T==np.diag(a)+np.diag(b, 1)+np.diag(c, -1)).all(), \
"Matrix must be Tri-diagonal"
assert (0 <= c*b).all(), "Matrix must have only real eigenvalues!"
g = np.zeros(N, float)
g[0] = 1.0
for n in range(N-1):
if b[n] == 0 or c[n] == 0:
g[n+1] = g[n]
else:
g[n+1] = g[n]*np.sqrt(c[n].conjugate()/b[n])
g /= g[N/2]
G = np.diag(g)
Ginv = np.diag(1./g)
GinvTG = Ginv*T*G
GinvTG = GinvTG + GinvTG.conjugate().transpose()
GinvTG /= 2.0
return (GinvTG, G, Ginv)
def mesh(min, max, scales, centers=0.0, N=100):
r"""Return a mesh for x in [min, max] inclusive with selective
sampling about the centers with specified length scales. The last
argument N may be an array specifying the number of points in each
interval.
For example, mesh(-inf, inf, [1E-3, 1, 2], [0, 0, 10], N=100) will return
a union of tangent meshes, each with 100 point, with the following
intervals well sampled: [-1E-3, 1E-3], [-1, 1] and [8, 12].
"""
if np.iterable(scales):
Nmesh = len(scales)
elif np.iterable(centers):
Nmesh = len(centers)
else:
Nmesh = 1
if not np.iterable(scales):
scales = np.ones(Nmesh)*scales
if not np.iterable(centers):
centers = np.ones(Nmesh)*centers
if not np.iterable(N):
N = np.ones(Nmesh)*N
mesh = np.array([], float)
for n in range(Nmesh):
new_mesh = tanMesh(min=min,
max=max,
N=N[n],
lengthScale=scales[n],
origin=centers[n])
mesh = np.unique(np.union1d(mesh, new_mesh))
mesh.sort()
mesh[0] = min
mesh[-1] = max
return mesh
def tanMesh(min=0, max=np.inf, N=100, lengthScale=1.0, origin=0.0):
r"""Return an array of N+2 points for x \in [min, max] using
x = origin + lengthScale*tan(y)
where y are equally spaced. The endpoints are included.
"""
assert 0 < lengthScale
def x_f(y):
return origin + lengthScale*np.tan(y)
def y_f(x):
return np.arctan((x-origin)/lengthScale)
if min == -np.inf:
y_min = -np.pi/2.
else:
y_min = y_f(min)
if max == np.inf:
y_max = np.pi/2.
else:
y_max = y_f(max)
assert y_min < y_max
y = np.linspace(y_min, y_max, N+2)
x = x_f(y)
x[0] = min
x[-1] = max
return x