lorenz equations http://youtu.be/QloSRF7Nj_o

lorenz.py

# -*- coding: utf-8 -*-
from tomIntegrate import intor, randr
from matplotlib import pyplot as plt
from numpy import array, polyfit, log, exp

def dxdt(X, r=28.):
    s = 10.
    b = 8. / 3.
    x = X[0]
    y = X[1]
    z = X[2]
    F = lambda x: x ** 3 / 3 - x
    return array([
                    s * (y - x),
                    x * (r - z) - y,
                    x * y - b * z
                  ])

def demo():
    import mpl_toolkits.mplot3d.axes3d as p3
    fig=plt.figure()
    ax = p3.Axes3D(fig)
    for i in range(14):
        minp = -10
        maxp = 10
        x0 = randr(minp, maxp)
        y0 = randr(minp, maxp)
        z0 = randr(0, maxp)
        (x, y, z), times = intor((x0, y0, z0), dxdt, tmax=4, nstep=444)
        ax.plot(x,y,z)
        #plt.plot(x, y)
    plt.show()

def altrDemo(r=166.3):
    def wrapr(func, r):
        def wrappor(X, t=0):
            return func(X, r=r)
        return wrappor
    diff = wrapr(dxdt, r)
    
    import mpl_toolkits.mplot3d.axes3d as p3
    fig=plt.figure()
    ax = p3.Axes3D(fig)
    for i in range(1):
        minp = -10
        maxp = 10
        x0 = randr(minp, maxp)
        y0 = randr(minp, maxp)
        z0 = randr(0, maxp)
        (x, y, z), times = intor((x0, y0, z0), diff, tmax=48, nstep=16000)
        ax.plot(x,y,z)
        extra = "($r=%.1f$)" % r
        if r == 166.3:
            title = "Intermittent Chaos"
            fig.suptitle(title + '\n' + extra)
        if r == 212:
            title = "Noisy Periodicity"
            fig.suptitle(title + '\n' + extra)
        #plt.plot(x, y)
    fig.savefig(title + "-Phase-" + ".pdf")
    fig2 = plt.figure()
    ax2 = fig2.add_subplot(1, 1, 1)
    ax2.plot(times, x, label=r"$x$")
    ax2.plot(times, y, label=r"$y$")
    ax2.plot(times, z, label=r"$z$")
    ax2.set_xlabel(r"$t$")
    ax2.legend(loc="best")
    fig2.suptitle(title + '\n' + extra)
    fig2.savefig(title + "-Timecourses-" + ".pdf")
    plt.show()

def hw():
    import mpl_toolkits.mplot3d.axes3d as p3
    N = 44444
    tmax = 34
    fig=plt.figure(figsize=(11, 8.5))
    ax = p3.Axes3D(fig)
    
    #Start with a random IC, and integrate for a little bit to get a new IC "on" the attactor.
    (x, y, z), times = intor((-5.817, -10.203, 13.087), dxdt, tmax=1, nstep=N)
        # not this:
        #a0 = (-5.817, -10.203, 13.087)
        #b0 = (-5.818, -10.203, 13.087)
    a0 = (x[-1], y[-1], z[-1])
    b0 = list(a0)
    b0[0] += 1e-15
    b0 = tuple(b0)
    initials = [a0, b0]
    timecourses = [[],[]]
    colors = ['r', 'b']
    labels = [r"$\alpha$", r"$\beta$"]
    for i in range(len(initials)):
        x0 = initials[i][0]
        y0 = initials[i][1]
        z0 = initials[i][2]
        (x, y, z), times = intor((x0, y0, z0), dxdt, tmax=tmax, nstep=N)
        ax.plot(x,y,z, colors[i], label=labels[i])
        timecourses[i] = [x, y, z]
        ax.scatter([x[-1]],[y[-1]],[z[-1]], color=colors[i], label=labels[i])
    ax.scatter([x0],[y0],[z0], color='k', label="initial")
    ax.plot([],[],'ko', label="initial")
    ax.set_xlabel(r"$x$")
    ax.set_ylabel(r"$y$")
    ax.set_zlabel(r"$z$")
    fig.suptitle("Lorenz Phase Portrait" + '\n' + r"$\alpha_0=(%.3f, %.3f, %.3f)$" % a0
                                         + '\n' + r"$\beta_0=(%.3f, %.3f, %.3f)$" % b0)
    ax.legend(loc="best")
    fig.savefig("448-hw5-3-phasespace.pdf")
        
    
    def vecLength(vec):
        squared = [x**2 for x in vec]
        return (sum(squared)) ** 0.5
    
    def vecDiff(vec1, vec2):
        if type(vec1) != type(vec2) != tuple:
            print "vectors not both tuples:", type(vec1), type(vec2)
            return -1
        if len(vec1) != len(vec2):
            print "vectors of different size!"
            return -1
        else:
            diff = [x - y for x, y in zip(list(vec1), list(vec2))]
            return diff
    
    a = range(N)
    b = range(N)
    dcourse = range(N)
    for i in range(N):
        #print i
        a[i] = (timecourses[0][0][i], timecourses[0][1][i], timecourses[0][2][i])
        b[i] = (timecourses[1][0][i], timecourses[1][1][i], timecourses[1][2][i])
        dcourse[i] = vecLength(vecDiff(a[i], b[i]))
    #print dcourse
    fig2 = plt.figure(figsize=(11, 8.5))
    ax2 = fig2.add_subplot(1, 1, 1)
    ax2.plot(times, dcourse, 'k-')
    maxlintime = 22
    maxind = int(N * maxlintime / tmax)
    print "max time for fitting:", times[maxind]
    [m, b] = polyfit(times[:maxind], [log(abs(d)) for d in dcourse][:maxind], 1)
    print "line fit:", "ln(y) = %.3f x + %.3f" % (m, b)
    fit = [exp(m * t + b) for t in times]
    ax2.plot(times, fit, "k--", label=r"exponential fit $ln(t)=%.3f t + %.3f$ (for $0 \leq t \leq %d$)" % (m, b, maxlintime))
    ax2.legend(loc="best")
    ax2.set_yscale("log")
    ax2.set_title(r"exponential growth of $\delta(t)$" + '\n'
                  + r"largest Liapunov exponent is $\lambda=%.3f$" % m)
    ax2.set_xlabel(r"$t$")
    ax2.set_ylabel("$ln||\delta(times)||$")
    fig2.savefig("448-hw5-3-errorgrowth.pdf")
    fig2.savefig("448-hw5-3-errorgrowth.png")
    plt.show()

if __name__=="__main__":
    altrDemo(r=166.3)
    altrDemo(r=212)
    #hw()
    #demo()

tomIntegrate.py

# -*- coding: utf-8 -*-
from numpy import array, random, linspace, arange
import pylab as p
from scipy import integrate

###Integrating the ODE using scipy.integrate
def doint(initial, dxdt, tmin=0, tmax=800, nstep=1e5):
    def dxdtextra(X, t=0):
        return dxdt(X)
    t = linspace(tmin, tmax, nstep)
    X0 = array(initial)
    X, infodict = integrate.odeint(dxdtextra, X0, t, full_output=True)
    return (X, t)

def intor(initial, dxdt, **kwargs):
    X, t = doint(initial, dxdt, **kwargs)
    return X.T, t

def randr(a, b):
    if a > b:
        print "b > a, dumpass"
        return 4  # http://xkcd.com/221
    else:
        return random.rand() * (b - a) + a

if __name__=="__main__":
    # Example 7.5.1
    u = 4.0
    def F(x):
        return x ** 3 / 3 - x
    def vdp(X, t=0):
        x = X[0]
        y = X[1]
        return array([
                        u * (y - F(x)),
                        -x / u
                    ])

    numTrials = 10
    for i in range(numTrials):
        initial = (randr(-4, 4), randr(-4, 4))
        position, rate = doint(initial, vdp).T
        p.plot(position, rate)
        p.scatter(initial[0], initial[1])
    fullrange = list(arange(-2.5, 2.5, .01))
    p.plot(fullrange, [F(x) for x in fullrange], 'k--')
    p.show()