Homec4science
Diffusion Lyapunov (master) Browse

Lyapunov/lyapunov/lyapunov.py
0e7b2e34bcfdmaster

Last Change
Raw File

lyapunov.py
View Options

import numpy as np
def lyapunov(X, r=2.0**(1/16), dt=1.0, isotropic=True, d0=None, radts=False,
interp=True):
"""
Compute the finite-size Lyapunov exponent of a set of trajectories.
X: Array-like of dimensions (L, K, M, N), with:
L the number of repeated experiment,
K the number of time snapshots,
M the number of spatial dimensions,
N the number of points.
r: Coefficient defining the sequence of thresholds, based on the initial
scale of the points cloud.
dt: Time sampling rate, by default 1.0.
isotropic: Consider the isotropic spreading of the cloud.
radts: Return also the time series of the radiuses (default: False).
interp: Use linear interpolation to compute doubling time (default: True)
See Artale et al., Physics of Fluids, 1997
Note that Aurell et al. 1996 used a different definition of FSLE!
"""
X = np.asanyarray(X)
if X.ndim != 4:
raise ValueError("The input array must be 4-D.")
nexp, ntime, ndim, npart = X.shape
if isotropic:
radius = np.mean(np.power(
np.linalg.norm(X - X.mean(axis=3)[:, :, :, np.newaxis], axis=2),
2),
axis=2)
else:
raise NotImplementedError
radius = np.mean(np.power(X - X.mean(axis=3)[:, :, :, np.newaxis], 2),
axis=3)
radius = np.sqrt(radius)
# estimate of the initial size
if d0 is None:
d0 = radius[..., 0].mean()
nmax = np.int32(np.floor(np.log(radius.max() / d0) / np.log(r)))
if nmax < 1:
raise ValueError("Insufficient dynamic range of the spreading, "
"try with longer time series or more particles.")
delta_n = d0 * np.power(r, np.arange(0, nmax + 1))
if isotropic:
Ln = np.zeros(nmax)
Lnstd = np.zeros(nmax)
Tn_old = np.zeros(nexp)
logr = np.log(r)
if not interp:
Dn_old = np.ones(nexp) * d0
else:
Ln = np.zeros((nmax, ndim))
Lnstd = np.zeros((nmax, ndim))
Tn_old = np.zeros((nexp, ndim))
logr = np.log(r)
if not interp:
Dn_old = np.ones((nexp, ndim)) * d0
for nn, dn in enumerate(delta_n[1:]):
# argmax returns zero if not found
tind = np.argmax(radius >= dn, axis=1)
notfound = tind == 0
# nn is shifted with respect to delta_n
if isotropic:
if not interp:
Tn = tind * dt
Dn = radius[range(nexp), tind]
Tn[notfound] = np.nan
Dn[notfound] = np.nan
else:
Tn1 = tind * dt
Tn0 = Tn1 - dt
Dn1 = radius[range(nexp), tind]
Dn0 = radius[range(nexp), tind-1]
# linear interpolation
Tn = Tn0 + (dn - Dn0) * (Tn1 - Tn0) / (Dn1 - Dn0)
Tn[notfound] = np.nan
else:
Dn = np.log(radius[range(nexp), tind] / delta_n[nn])
if not interp:
# for formulas, see Ku 1966
XX = Tn - Tn_old
X = 1/np.nanmean(XX)
vX = np.nanvar(XX) / np.power(X, -4)
YY = np.log(Dn / Dn_old)
Y = np.nanmean(YY)
vY = np.nanvar(YY)
Ln[nn] = X * Y
Lnstd[nn] = X * Y * np.sqrt(vX/X/X + vY/Y/Y)
else:
XX = Tn - Tn_old
X = 1/np.nanmean(XX)
sX = np.nanstd(XX) / np.power(X, -2)
Ln[nn] = X * logr
Lnstd[nn] = sX * logr
# store current results for next cycle
Tn_old[:] = Tn
if not interp:
Dn_old[:] = Dn
# end cycle on delta_n
if radts:
return delta_n[:-1], Ln, Lnstd, radius
else:
return delta_n[:-1], Ln, Lnstd
def standard_map(X0, K, n):
nprt = X0.shape[1]
x = np.empty((n, nprt))
y = np.empty((n, nprt))
x[0] = X0[0, :]
y[0] = X0[1, :] % (2*np.pi)
for nn in range(1, n):
x[nn, :] = x[nn-1, :] + K * np.sin(y[nn-1, :])
y[nn, :] = (y[nn-1, :] + x[nn, :]) % (2*np.pi)
return x,y
def standard_map_reflecting(X0, K, l, L, n):
nprt = X0.shape[1]
from scipy.optimize import minimize_scalar
def f(x):
if np.abs(x)<l:
return 1
else:
return (L-np.abs(x)) / (L - l)
def fprime(x):
if np.abs(x) < l:
return 0
else:
return -1*np.sign(x)
def xnp1(x, xn, yn):
return xn + K * f(x) * np.sin(yn)
x = np.empty((n, nprt))
y = np.empty((n, nprt))
x[0] = X0[0, :]
y[0] = X0[1, :] % (2*np.pi)
for nn in range(1, n):
for kk in range(nprt):
res = minimize_scalar(xnp1, args=(x[nn-1, kk], y[nn-1, kk]))
x[nn, kk] = res.x
y[nn, kk] = (y[nn-1, kk] + x[nn, kk] - K * fprime(x[nn-1, kk] * np.cos(y[nn-1, kk]))) % (2*np.pi)
return x,y
c4science · Help