# 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 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; }