In [62]:

    

from mpmath import mpf, mp
import pickle
import sys
%pylab inline --no-import-all
from pylab import *

# Constants for bits of significand precision
# TODO these should be for mp.prec, not mp.dps
SINGLE = 24
DOUBLE = 53

def vectorize(functions):
    """
    Takes a list of functions and returns a function that accepts a list of
    values and applies each function given to the list of values.
    """
    def _vectorized(x, ys):
        return tuple(f(x, *ys) for f in functions)

    return _vectorized

def inc_vec(vec, inc):
    return [x + inc for x in vec]

def add_vec(*vecs):
    return [sum(nums) for nums in zip(*vecs)]

def mul_vec(vec, mul):
    return [x * mul for x in vec]

def rk4_system(y, x0, y0, h, steps):
    """
    y: list of ODE
    y0: vector (tuple) of initial y values
    Finds y(x0 + h * steps)
    """
    x1 = x0
    y1 = list(y0)
    xs = [x0]
    ys = [[y0[i]] for i,f in enumerate(y)]

    f = vectorize(y)

    for i in range(steps):
        k1 = f(x1, y1)
        k2 = f(x1 + (h / 2), add_vec(y1, mul_vec(k1, h / 2)))
        k3 = f(x1 + (h / 2), add_vec(y1, mul_vec(k2, h / 2)))
        k4 = f(x1 + h, add_vec(y1, mul_vec(k3, h)))

        y1 = add_vec(y1, mul_vec(add_vec(k1, mul_vec(k2, 2), mul_vec(k3, 2), k4), h / 6))
        for i, val in enumerate(y1):
            ys[i].append(val)
        x1 += h
        xs.append(x1)

    return xs, ys


def runrange(file_suffix, a15rate, precision, steps):
    IS_MPMATH = isinstance(precision, int)
    print("Using mpmath:", IS_MPMATH)
    if IS_MPMATH:
        a1 = mpf('2')
        a2 = mpf('2')
        a3 = mpf('3')
        a4 = mpf('0.7')
        a5 = mpf('1')
        a6 = mpf('0.02')
        a7 = mpf('20')
        a8 = mpf('1')
        a9 = mpf('1')
        a10 = mpf('1')
        a11 = mpf('0.01')
        a12 = mpf('0.1')
        a13 = mpf('0.3')
        a14 = mpf('50')
        a16 = mpf('0.1')
        a17 = mpf('0.14')
        B = mpf('1')
        start = mpf('0')
        initialconds = (mpf('0.5'), mpf('0'), mpf('1'), mpf('0.1'))
        delta = mpf('0.1')
    
        mp.dps = precision
        print("Precision:", mp.dps)
    else:
        a1 = 2.0
        a2 = 2.0
        a3 = 3.0
        a4 = 0.7
        a5 = 1.0
        a6 = 0.02
        a7 = 20.0
        a8 = 1.0
        a9 = 1.0
        a10 = 1.0
        a11 = 0.01
        a12 = 0.1
        a13 = 0.3
        a14 = 50.0
        a16 = 0.1
        a17 = 0.14
        B = 1
        start = 0.0
        initialconds = (0.5, 0.0, 1.0, 0.1)
        delta = 0.1
        print("Precision: Python float (double)", sys.float_info)
    f1 = lambda p: p**a1 / (a2**a1 + p**a1)
    f2 = lambda p: a4 * a5**a3 / (a5**a3 + p**a3)
    p = lambda E, a15: (a14 * E * B) / a15
    
    y = [
        lambda t, A, M, E, a15: f1(p(E,a15))*(a6+a7*M)-a8*A-a9*A**2, #dA/dt
        lambda t, A, M, E, a15: a10*f2(p(E,a15))*A-f1(p(E,a15))*a16*a7*M-a11*M, #dM/dt
        lambda t, A, M, E, a15: a12*(1-f2(p(E,a15)))*A-a13*E, #dE/dt,
        lambda t, A, M, E, a15: a15rate
    ]
    xs, ys = rk4_system(y, start, initialconds, delta, steps)
    pickle.dump(list(zip(xs, *ys)), open('time_py_np_{}.dat'.format(file_suffix), 'wb'))


    

    




Populating the interactive namespace from numpy and matplotlib

In [42]:

    

runrange('slow_single', 0.0002, float, 150000)
#runrange('fast_single', mpf('0.01'), SINGLE)


    

    




Using mpmath: False
Precision: Python float (double) sys.float_info(max=1.7976931348623157e+308, max_exp=1024, max_10_exp=308, min=2.2250738585072014e-308, min_exp=-1021, min_10_exp=-307, dig=15, mant_dig=53, epsilon=2.220446049250313e-16, radix=2, rounds=1)

In [46]:

    

runrange('slow_100', mpf('0.0002'), 100, 150000)
#runrange('fast_double', mpf('0.01'), DOUBLE)


    

    




Using mpmath: True
Precision: 100

In [47]:

    

ssingle = pickle.load(open('time_py_np_slow_single.dat', 'rb'))
sdouble = pickle.load(open('time_py_np_slow_100.dat', 'rb'))


    

In [93]:

    

matplotlib.rcParams.update({'font.size': 14})
matplotlib.rcParams.update({'figure.autolayout': True})


    

In [94]:

    

figure()
title("Python float (double precision; 53 bits)")
plot([float(x[4]) for x in ssingle], [float(x[1]) for x in ssingle])
xlabel("delta_p")
ylabel("A")
axvline(0.5707, color='black')
axhline(color='black')
annotate("Hopf bifurcation", (0.5707, 2), xytext=(0.8,2.5), arrowprops=dict(facecolor='black',lw=0.01))
savefig("precision_double.png")


    

In [97]:

    

figure()
title("mpmath float (336 bits precision)")
ylim(0,3)
plot([float(x[4]) for x in sdouble], [float(x[1]) for x in sdouble])
xlabel("delta_p")
ylabel("A")
axvline(0.5707, color='black')
axhline(color='black')
annotate("Hopf bifurcation", (0.5707, 2), xytext=(0.8,2.5), arrowprops=dict(facecolor='black',lw=0.01))
savefig("precision_100.png")


    

In [ ]: