lyap.cpp

#include <iostream>
#include <cmath>
#include <vector>

#define _USE_MATH_DEFINES 
#define NEQ 2
#define G 1
#define L 1

using namespace std;

double q = 0.;
double f = 0.;
double ohm = 0.;
double theta0 = 0.001;
double omega0 = 0.;
int num_steps = 65536;

void dYdt(double t, double* Y, double* R) {
    double theta = Y[0];
    double omega = Y[1];
    R[0] = omega;
    R[1] = -(G/L)*sin(theta) - q*omega + f*sin(ohm*t);
}

void Jacobian(double t, double* Y, double J[NEQ][NEQ]) {
    double theta = Y[0];
    J[0][0] = 0;
    J[0][1] = 1;
    J[1][0] = -(G/L)*cos(theta);
    J[1][1] = -q;
}

void RungeKutta4(double t, double* Y, void (*RHSFunc)(double, double*, double*), double dt, int neq) {
    double k1[NEQ], k2[NEQ], k3[NEQ], k4[NEQ], Ytemp[NEQ];

    RHSFunc(t, Y, k1);
    for (int k = 0; k < neq; k++) {
        Ytemp[k] = Y[k] + 0.5 * dt * k1[k];
    }
    RHSFunc(t + 0.5 * dt, Ytemp, k2);
    for (int k = 0; k < neq; k++) {
        Ytemp[k] = Y[k] + 0.5 * dt * k2[k];
    }
    RHSFunc(t + 0.5 * dt, Ytemp, k3);
    for (int k = 0; k < neq; k++) {
        Ytemp[k] = Y[k] + dt * k3[k];
    }
    RHSFunc(t + dt, Ytemp, k4);
    for (int k = 0; k < neq; k++) {
        Y[k] += (1.0/6.0) * dt * (k1[k] + 2.0*k2[k] + 2.0*k3[k] + k4[k]);
    }
}

void GramSchmidt(double* v1, double* v2) {
    double dot = v1[0] * v2[0] + v1[1] * v2[1];
    v2[0] -= dot * v1[0];
    v2[1] -= dot * v1[1];
    
    double norm1 = sqrt(v1[0] * v1[0] + v1[1] * v1[1]);
    double norm2 = sqrt(v2[0] * v2[0] + v2[1] * v2[1]);
    
    v1[0] /= norm1;
    v1[1] /= norm1;
    v2[0] /= norm2;
    v2[1] /= norm2;
}

void LyapunovExponents(double* Y0, int num_steps, double dt) {
    double Y[NEQ] = {Y0[0], Y0[1]};
    double v1[NEQ] = {1.0, 0.0};
    double v2[NEQ] = {0.0, 1.0};

    double t = 0.0;
    double sumLogNorm1 = 0.0, sumLogNorm2 = 0.0;

    for (int step = 0; step < num_steps; ++step) {
        // Integrate the main system
        RungeKutta4(t, Y, dYdt, dt, NEQ);

        // Integrate tangent vectors and apply Gram-Schmidt
        double J[NEQ][NEQ], Vtemp[NEQ];
        Jacobian(t, Y, J);
        for (int i = 0; i < NEQ; i++) {
            Vtemp[i] = 0;
            for (int j = 0; j < NEQ; j++) {
                Vtemp[i] += J[i][j] * v1[j];
            }
        }
        for (int i = 0; i < NEQ; i++) v1[i] = Vtemp[i];

        for (int i = 0; i < NEQ; i++) {
            Vtemp[i] = 0;
            for (int j = 0; j < NEQ; j++) {
                Vtemp[i] += J[i][j] * v2[j];
            }
        }
        for (int i = 0; i < NEQ; i++) v2[i] = Vtemp[i];

        GramSchmidt(v1, v2);

        // Compute norms for Lyapunov exponents
        sumLogNorm1 += log(sqrt(v1[0]*v1[0] + v1[1]*v1[1]));
        sumLogNorm2 += log(sqrt(v2[0]*v2[0] + v2[1]*v2[1]));

        // Normalize vectors
        double norm1 = sqrt(v1[0]*v1[0] + v1[1]*v1[1]);
        double norm2 = sqrt(v2[0]*v2[0] + v2[1]*v2[1]);
        v1[0] /= norm1;
        v1[1] /= norm1;
        v2[0] /= norm2;
        v2[1] /= norm2;

        t += dt;
    }

    double lambda1 = sumLogNorm1 / (num_steps * dt);
    double lambda2 = sumLogNorm2 / (num_steps * dt);

    cout << "Lyapunov Exponent 1: " << lambda1 << endl;
    cout << "Lyapunov Exponent 2: " << lambda2 << endl;
}

int main() {
    double t_start = 0.;
    double t_end = 50 * M_PI;
    double dt = (t_end - t_start) / num_steps;

    double Y0[NEQ] = {0.02, 0.1};

    f = 1.2;
    q = 0.5;
    ohm = 2. / 3.;

    LyapunovExponents(Y0, num_steps, dt);

    return 0;
}