diff --git a/README.md b/README.md index d756b2a..30f68d0 100644 --- a/README.md +++ b/README.md @@ -40,7 +40,7 @@ data = l63_obj.simulate(5000) # -> data.shape = (5000,3) iterator = l63_obj.iterate lle = lpy.measures.largest_lyapunov_exponent( iterator_func=iterator, - starting_point=l63_obj.default_starting_point, + starting_point=l63_obj.get_default_starting_pnt(), dt=l63_obj.dt ) # -> lle = 0.905144329... diff --git a/src/lorenzpy/measures.py b/src/lorenzpy/measures.py index 69850fe..dfb0a1c 100644 --- a/src/lorenzpy/measures.py +++ b/src/lorenzpy/measures.py @@ -74,3 +74,95 @@ def largest_lyapunov_exponent( ) else: return float(np.average(log_divergence) / (dt * part_time_steps)) + + +def lyapunov_exponent_spectrum( + iterator_func: Callable[[np.ndarray], np.ndarray], + starting_point: np.ndarray, + deviation_scale: float = 1e-10, + steps: int = int(1e3), + part_time_steps: int = 15, + steps_skip: int = 50, + dt: float = 1.0, + m: int | None = None, + initial_pert_directions: np.ndarray | None = None, + return_convergence: bool = False, +) -> np.ndarray: + """Calculate the spectrum of m largest lyapunov exponent given an iterator function. + + A mixture of: + - The algorithm for the largest lyapunov exponent: Sprott, Julien Clinton, and + Julien C. Sprott. Chaos and time-series analysis. Vol. 69. Oxford: Oxford + university press, 2003. + - The algorithm for the spectrum given in 1902.09651 "LYAPUNOV EXPONENTS of the + KURAMOTO-SIVASHINSKY PDE". + + Args: + iterator_func: Function to iterate the system to the next time step: + x(i+1) = F(x(i)) + starting_point: The starting_point of the main trajectory. + deviation_scale: The L2-norm of the initial perturbation. + steps: Number of renormalization steps. + part_time_steps: Time steps between renormalization steps. + steps_skip: Number of renormalization steps to perform, before tracking the log + divergence. + Avoid transients by using steps_skip. + dt: Size of time step. + m: Number of Lyapunov exponents to compute. If None: take all (m = x_dim). + initial_pert_directions: + - If np.ndarray: The direction of the initial perturbations of shape + (m, x_dim). Will be qr decomposed + and the q matrix will be used. + - If None: The direction of the initial perturbation is assumed to be + np.eye(x_dim)[:m, :]. + return_convergence: If True, return the convergence of the largest LE; a + numpy array of the shape (N, m). + + Returns: + The Lyapunov exponent spectrum of largest m values. If return_convergence is + True: The convergence (2D N x m np.ndarray), else a (1D m-size np.ndarray), + which holds the last values in the convergence. + """ + x_dim = starting_point.size + + if m is None: + m = x_dim + + if initial_pert_directions is None: + initial_perturbations = ( + np.eye(x_dim)[:m, :] * deviation_scale + ) # each vector is of length deviation_scale + else: + q, _ = np.linalg.qr(initial_pert_directions) + initial_perturbations = q * deviation_scale + + log_divergence_spec = np.zeros((steps, m)) + + x = starting_point + x_pert = starting_point + initial_perturbations + + for i_n in range(steps + steps_skip): + for i_t in range(part_time_steps): + x = iterator_func(x) + for i_m in range(m): + x_pert[i_m, :] = iterator_func(x_pert[i_m, :]) + dx = x_pert - x + q, r = np.linalg.qr(dx.T, mode="reduced") + x_pert = x + q.T * deviation_scale + + if i_n >= steps_skip: + log_divergence_spec[i_n - steps_skip, :] = np.log( + np.abs(np.diag(r)) / deviation_scale + ) + + if return_convergence: + return np.array( + np.cumsum(log_divergence_spec, axis=0) + / ( + np.repeat(np.arange(1, steps + 1)[:, np.newaxis], m, axis=1) + * dt + * part_time_steps + ) + ) + else: + return np.average(log_divergence_spec, axis=0) / (dt * part_time_steps) diff --git a/tests/test_measures.py b/tests/test_measures.py index ee01b30..309531a 100644 --- a/tests/test_measures.py +++ b/tests/test_measures.py @@ -1,19 +1,35 @@ """Testing simulations.""" + import numpy as np +import pytest -from lorenzpy import measures +from lorenzpy import measures, simulations -def test_largest_lyapunov_one_d_linear_time_dependent(): - """Testing measures.largest_lyapunov_exponent.""" +@pytest.fixture +def get_demo_exponential_decay_system(): + """Fixture to get a demo iterator func and starting point. + + Returns the iterator_func, starting_point, dt and a for exponential decay. + """ dt = 0.1 a = -0.1 def iterator_func(x): return x + dt * a * x - return_convergence = False starting_point = np.array(1) + + return iterator_func, starting_point, dt, a + + +def test_largest_lyapunov_one_d_linear_time_dependent( + get_demo_exponential_decay_system, +): + """Testing measures.largest_lyapunov_exponent.""" + iterator_func, starting_point, dt, a = get_demo_exponential_decay_system + + return_convergence = False deviation_scale = 1e-10 steps = 10 steps_skip = 1 @@ -33,16 +49,13 @@ def iterator_func(x): np.testing.assert_almost_equal(actual, desired, 2) -def test_largest_lyapunov_one_d_linear_time_dependent_return_conv(): +def test_largest_lyapunov_one_d_linear_time_dependent_return_conv( + get_demo_exponential_decay_system, +): """Testing measures.largest_lyapunov_exponent with convergence.""" - dt = 0.1 - a = -0.1 - - def iterator_func(x): - return x + dt * a * x + iterator_func, starting_point, dt, a = get_demo_exponential_decay_system return_convergence = True - starting_point = np.array(1) deviation_scale = 1e-10 steps = 10 steps_skip = 1 @@ -65,3 +78,93 @@ def iterator_func(x): * steps ) np.testing.assert_almost_equal(actual, desired, 2) + + +def test_lyapunov_spectrum_lorenz63(): + """Testing lyapunov_exponent_spectrum with initial_pert_directions.""" + dt = 0.05 + lorenz63_obj = simulations.Lorenz63(dt=dt) + iterator_func = lorenz63_obj.iterate + starting_point = lorenz63_obj.get_default_starting_pnt() + + m = 3 + deviation_scale = 1e-10 + steps = 1000 + part_time_steps = 15 + steps_skip = 50 + initial_pert_directions = None + return_convergence = False + + lyap_spec_conv = measures.lyapunov_exponent_spectrum( + iterator_func=iterator_func, + starting_point=starting_point, + deviation_scale=deviation_scale, + steps=steps, + part_time_steps=part_time_steps, + steps_skip=steps_skip, + dt=dt, + m=m, + initial_pert_directions=initial_pert_directions, + return_convergence=return_convergence, + ) + + expected_lyapunovs = np.array([9.05e-01, 1.70e-03, -1.44e01]) + np.testing.assert_almost_equal(lyap_spec_conv, expected_lyapunovs, decimal=1) + + +def test_lyapunov_spectrum_vs_largest_le_lorenz63(): + """Testing lyapunov_exponent_spectrum with initial_pert_directions.""" + dt = 0.05 + lorenz63_obj = simulations.Lorenz63(dt=dt) + iterator_func = lorenz63_obj.iterate + starting_point = lorenz63_obj.get_default_starting_pnt() + + m = 1 + deviation_scale = 1e-10 + steps = 500 + part_time_steps = 10 + steps_skip = 50 + return_convergence = True + + lyap_spec_conv = measures.lyapunov_exponent_spectrum( + iterator_func=iterator_func, + starting_point=starting_point, + deviation_scale=deviation_scale, + steps=steps, + part_time_steps=part_time_steps, + steps_skip=steps_skip, + dt=dt, + m=m, + initial_pert_directions=None, + return_convergence=return_convergence, + ) + + lyap_largest_conv = measures.largest_lyapunov_exponent( + iterator_func=iterator_func, + starting_point=starting_point, + deviation_scale=deviation_scale, + steps=steps, + part_time_steps=part_time_steps, + steps_skip=steps_skip, + dt=dt, + initial_pert_direction=None, + return_convergence=return_convergence, + ) + + np.testing.assert_almost_equal(lyap_spec_conv, lyap_largest_conv[:, np.newaxis], 3) + + +def test_lyapunov_spectrum_custom_pert(get_demo_exponential_decay_system): + """Test if it works if custom initial_pert_directions is given.""" + iterator_func, starting_point, dt, a = get_demo_exponential_decay_system + + lle = measures.lyapunov_exponent_spectrum( + iterator_func=iterator_func, + starting_point=starting_point, + dt=dt, + steps=10, + part_time_steps=1, + steps_skip=0, + initial_pert_directions=np.array([[1]]), + ) + np.testing.assert_almost_equal(lle, a, decimal=2)