README.rst

Simulation of the Kaplan–Yorke map using ctypes and numpy

Result:

Kaplan–Yorke map

Definition:

x_{n+1} = T(x_{n})
y_{n+1} = lambda y_{n} + h(x_{n})
T(x) = 1 - mu x^2
h(x) = x

See:

• Thermodynamics of chaotic systems ... - Google Books
• Kaplan–Yorke map - Wikipedia, the free encyclopedia

ctypes

• Cookbook/Ctypes -
• ctypes tutorial
• 15.16. ctypes — A foreign function library for Python — Python v2.7rc1 documentation

kaplan_yorke.c

typedef struct kaplanyorkemap_{
  int step;
  double *xt;
  double *yt;
  double mu;
  double lambda;
} KaplanYorkeMap;

int kym_gene_seq(KaplanYorkeMap *self)
{
  int st;
  for (st = 1; st < self->step; ++st){
    self->xt[st] = 1 - self->mu * self->xt[st-1] * self->xt[st-1];
    self->yt[st] = self->lambda * self->yt[st-1] + self->xt[st-1];
  }
  return 0;
}

Makefile

objs=kaplan_yorke.o
slib=kaplan_yorke.so

all: $(slib)

kaplan_yorke.so: kaplan_yorke.o
	gcc -shared -o $@ $<

kaplan_yorke.o: kaplan_yorke.c
	gcc -c $< -fPIC -O3 -ftree-vectorizer-verbose=1

clean:
	rm $(objs) $(slib)

.PHONY: all clean

sim.py

import ctypes
import numpy
import pylab

double1d = ctypes.POINTER(ctypes.c_double)

class KaplanYorkeMapStruct(ctypes.Structure):
    _fields_ = [
        ("step", ctypes.c_int),
        ("xt", double1d),
        ("yt", double1d),
        ("mu", ctypes.c_double),
        ("ld", ctypes.c_double),
        ]

_kaplan_yorke = numpy.ctypeslib.load_library('kaplan_yorke.so', '.')
_kaplan_yorke.kym_gene_seq.restype = ctypes.c_int
_kaplan_yorke.kym_gene_seq.argtypes = [ctypes.POINTER(KaplanYorkeMapStruct)]

def kym_gene_seq(kyms):
    return _kaplan_yorke.kym_gene_seq(ctypes.pointer(kyms))

if __name__ == '__main__':
    step = 10**5
    xt = numpy.zeros(step, dtype=float)
    yt = numpy.zeros(step, dtype=float)
    mu = 2
    ld = 0.4

    xt[0] = 0.1
    yt[0] = 0.1
    
    kyms = KaplanYorkeMapStruct(
        step = xt.shape[0],
        xt = xt.ctypes.data_as(double1d),
        yt = yt.ctypes.data_as(double1d),
        mu = mu,
        ld = ld)
    kym_gene_seq(kyms)

    pylab.figure(1, figsize=(4,3))
    pylab.plot(xt, yt, 'k.', markersize=0.1)
    pylab.savefig('kym.png')