Rcjp's Weblog

December 30, 2006

Lyaponov Exponent

Filed under: c — rcjp @ 6:42 pm

Most introductions to chaos theory will talk about the logistic mapping x_{n+1}=r\, x_n(1-x_n) and then calculate the lyapunov exponent which is a measure of how rapidly small deviations in the initial conditions give exponential changes in results giving an indication how chaotic the system is. It is calculated from \lambda = \frac{1}{N} \sum^N_{n=1} \log(\frac{d x_{n+1}}{d x_n}) The program below calculates this for 2.95<r<4.0 after iterating for 1000 times to stabalise things then using N=30.

The code is slightly contrived; I wrote it to play with passing functions through C++ templates. The idea being that logis and deriv (its derivative) since they are templated, can get inlined by the compiler; though Id’ guess a much more straightforward C program would probably optimise just as well. In addition the std::bind1st(std::mem_fun(&IterMap::lyaponov_func),this) didn’t exactly trip off my fingers.

#include <iostream>
#include <fstream>
#include <algorithm>
#include <vector>
#include <numeric>
#include <functional>
#include <cmath>

typedef double (*CALC)(double, double);

class IterMap
{
    std::vector<double> data;
    const int nstabalise;
    const int nrecord;
    double r;
 public:
    IterMap(int rec, int stab=1000) : nstabalise(stab), nrecord(rec) {}
    template<CALC derivfunc>double lyaponov_func(double x);
    template<CALC mapfunc> void stabalise(double x);
    template<CALC mapfunc, CALC derivfunc> double lyapunov_exp(double r);
};

template<CALC mapfunc> void IterMap::stabalise(double x)
{
    for(int i = 0; i < nstabalise; ++i)
        x = mapfunc(r, x);

    for(int i = 0; i < nrecord; ++i) {
        x = mapfunc(r, x);
        data.push_back(x);
    }
}

template<CALC derivfunc> double IterMap::lyaponov_func(double x)
{
    return log(fabs( derivfunc(r,x) ));
}

template<CALC mapfunc, CALC derivfunc> double IterMap::lyapunov_exp(double ra)
{
    r = ra;
    stabalise<mapfunc>(0.1);

    // used mem_func, else have to declare func static

    transform(data.begin(), data.end(), data.begin(),
            std::bind1st(std::mem_fun(&IterMap::lyaponov_func<derivfunc>),this));
    double avg = accumulate(data.begin(), data.end(), 0.0) / data.size();
    data.clear();
    return avg;
}

////////////////////////////////////////////////////////////////////////////////

double logis(double r, double x)
{
    return r * x * (1.0 - x);
}

double deriv(double r, double x)
{
    return r - 2.0*r*x;
}

int main(int argc, char* argv[])
{
    std::ofstream out("ldata", std::ios::out);
    IterMap bir(30);
    for (double r = 2.95; r <= 4.0; r+=0.0005)
        out << r << " " << bir.lyapunov_exp<logis,deriv>(r) << std::endl;

    return 0;
}

Advertisements

Leave a Comment »

No comments yet.

RSS feed for comments on this post.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Create a free website or blog at WordPress.com.

%d bloggers like this: