conceptor.chaos.Lyapunov.java Source code

Introduction

Here is the source code for conceptor.chaos.Lyapunov.java

Source

/*
 * To change this license header, choose License Headers in Project Properties.
 * To change this template file, choose Tools | Templates
 * and open the template in the editor.
 */

package conceptor.chaos;

import org.apache.commons.math3.util.MathArrays;

/**
   A class for computing the maximal lyapunov exponent of
   a dynamical system
 */
public class Lyapunov {

    /**
       Given a DynamicalSystem, returns the maximul lyapunov exponent
       using the method outline by J.C. Sprott see:
       http://sprott.physics.wisc.edu/chaos/lyapexp.htm
     */
    public static double computeLargestExponent(DynamicalSystem system) {

        double d0 = 10e-8;

        // First, evolve the system such that it's very likely to be on
        // the attractor (1000 iterations is a decent guess)
        evolveSystem(system, 10000);

        // Next, pick a test point a distance d0 away
        double[] testPoint = getTestPoint(system, d0);
        double[] x = system.getState();

        double d1 = 0.0;
        double l = 0.0;

        // Evolve both the test point and the original
        for (int i = 0; i < 64100; i++) {

            testPoint = system.evolveOther(testPoint);
            x = system.evolve();

            // Test for unbounded orbits
            if (isUnbounded(x)) {
                return Double.NaN;
            }

            // Compute distance
            d1 = MathArrays.distance(testPoint, x);

            // Compute lyapunov exponent this round
            if (i > 100) {
                l = l + Math.log(d1 / d0);
            }

            // Renormalize test point
            for (int j = 0; j < x.length; j++) {
                testPoint[j] = x[j] + (d0 / d1) * (testPoint[j] - x[j]);
            }
        }
        return l / (63999 * system.getStepSize());
    }

    private static void evolveSystem(DynamicalSystem system, int times) {
        for (int i = 0; i < times; i++) {
            system.evolve();
        }
    }

    private static double[] getTestPoint(DynamicalSystem system, double d0) {
        double[] testPoint = new double[system.getDimension()];
        double denom = Math.sqrt((new Integer(system.getDimension())).doubleValue());
        double[] x = system.getState();

        for (int i = 0; i < x.length; i++) {
            testPoint[i] = x[i] + d0 / denom;
        }

        return testPoint;
    }

    private static boolean isUnbounded(double[] x) {
        double absSum = 0.0;
        for (int i = 0; i < x.length; i++) {
            absSum += Math.abs(x[i]);
        }
        if (absSum >= 10e6) {
            return true;
        } else {
            return false;
        }
    }
}