In [33]:

    

import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate
from time import time
%matplotlib inline


    

In [3]:

    

def Runge_Kutta(f, t, x0):
    n = len(x0)
    N = len(t)
    x = np.zeros((n, N))
    x[:, 0] = x0
    dt = t[1] - t[0]
    for i in range(N - 1):
        k1 = f(x[:, i])
        k2 = f(x[:, i] + dt * k1 / 2)
        k3 = f(x[:, i] + dt * k2 / 2)
        k4 = f(x[:, i] + dt * k3)
        x[:, i + 1] =  x[:, i] + dt * (k1 + 2 *  (k2 + k3) + k4) / 6
    return x


    

In [18]:

    

def Lorenz(X, t=0):
    return np.array([10 * (X[1] - X[0]),
                     28 * X[0] - X[1] - X[0] * X[2],
                     X[0] * X[1] - 8 * X[2] / 3])


    

In [19]:

    

N = 10 ** 4
t = np.linspace(0, 50, N)
x0 = np.array([1, 0, 0])
x = Runge_Kutta(Lorenz, t, x0)
plt.plot(t, x.T)
plt.xlabel('t')
plt.ylabel('x, y, z')


    

    
Out[19]:





[<matplotlib.lines.Line2D at 0x7f988bf87c18>,
 <matplotlib.lines.Line2D at 0x7f988bf87dd8>,
 <matplotlib.lines.Line2D at 0x7f988bf87f98>]






    

In [20]:

    

def Lorenz100(X, t=0):
    return np.hstack((10 * (X[100:200] - X[:100]),
                      28 * X[:100] - X[100:200] - X[:100] * X[200:300],
                      X[:100] * X[100:200] - 8 * X[200:300] / 3))


    

In [38]:

    

N = 10 ** 5
t = np.linspace(0, 50, N)
x0 = np.hstack((np.random.normal(1, .05, size=100), np.zeros(200)))
current_time = time()
x = Runge_Kutta(Lorenz100, t, x0)
print('It took', time() - current_time, 's')


    

    




It took 69.72971725463867 s

In [43]:

    

x1 = np.mean(x[:100, :], axis=0)
x2 = np.mean(x[100:200, :], axis=0)
x3 = np.mean(x[200:300, :], axis=0)


    

In [44]:

    

plt.plot(t, x1)
plt.plot(t, x2)
plt.plot(t, x3)
plt.xlabel('t')
plt.ylabel('mean of x, mean of y, mean of z')


    

    
Out[44]:





<matplotlib.text.Text at 0x7f988a588d30>






    

In [45]:

    

var_x1 = np.var(x[:100, :], axis=0)
var_x2 = np.var(x[100:200, :], axis=0)
var_x3 = np.var(x[200:300, :], axis=0)

plt.plot(t, var_x1)
plt.plot(t, var_x2)
plt.plot(t, var_x3)
plt.ylabel('var of x, var of y, var of z')


    

    
Out[45]:





<matplotlib.text.Text at 0x7f988fa6a0b8>






    

In [30]:

    

def trya():
    def Lorenz(x, y, z):
        return np.array([10 * (y - x),
                         28 * x - y - x * z,
                         x * y - 8 * z / 3])
    N = 10 ** 5
    t = np.linspace(0, 50, N)
    x0 = np.array([1, 0, 0])
    n = len(x0)
    N = len(t)
    x = np.zeros(N)
    y = np.zeros(N)
    z = np.zeros(N)
    #x = x0
    dt = t[1] - t[0]
    x[0] = 1
    for i in range(N - 1):
        k1 = Lorenz(x[i], y[i], z[i])
        k2 = Lorenz(x[i] + dt * k1[0] / 2, y[i] + dt * k1[1] / 2, z[i] + dt * k1[2] / 2)
        k3 = Lorenz(x[i] + dt * k2[0] / 2, y[i] + dt * k2[1] / 2, z[i] + dt * k2[2] / 2)
        k4 = Lorenz(x[i] + dt * k3[0], y[i] + dt * k3[1], z[i] + dt * k3[2])

        x[i + 1] = x[i] + dt * (k1[0] + 2 *  (k2[0] + k3[0]) + k4[0]) / 6
        y[i + 1] = y[i] + dt * (k1[1] + 2 *  (k2[1] + k3[1]) + k4[1]) / 6
        z[i + 1] = z[i] + dt * (k1[2] + 2 *  (k2[2] + k3[2]) + k4[2]) / 6
%timeit trya()


    

    




1 loop, best of 3: 11.2 s per loop

In [30]:

    

np.random.randn(100)


    

    
Out[30]:





array([-0.72716519, -0.02058777, -1.51594323, -1.18432998,  0.09269376,
       -0.4429285 , -0.64289323, -0.05178049, -0.66772672,  0.91812463,
       -0.95029808, -1.30591762, -1.46467793, -0.61783536, -1.34486354,
        0.54804895, -0.05145432, -1.89929532, -0.66263269,  0.19703072,
       -0.77757872,  0.14121466, -0.405104  ,  0.01251034,  0.67794063,
       -0.99397382,  0.56693904, -0.70627103, -0.33247358,  0.49343373,
       -0.00580554, -0.2744945 , -1.25305373,  0.41688131,  1.53116205,
        0.92867199, -2.93789449, -0.48727762, -0.29826876,  0.46173463,
       -0.49272694,  0.92998859,  0.54006121, -1.30415885, -0.59601325,
        0.22803473, -0.36864005,  2.85002914,  1.08893887, -1.31550696,
       -0.9639277 , -1.95050475, -0.39708963,  0.00557951,  0.73862203,
       -0.60658556,  0.84130735, -0.83125291, -0.75639711, -0.24626057,
        0.40253579, -0.134268  ,  0.02997676,  0.21124218, -0.63111478,
        0.02005445, -0.06615966,  1.77546336, -1.16328277, -1.29530291,
        0.36768733, -0.08235877,  0.17943627, -1.12435255,  0.23148413,
        0.16279953,  0.26016607,  0.02680808, -0.50076882, -1.07862074,
        0.78676639,  0.09585915,  0.55780738,  0.2187324 ,  0.39706541,
        0.6786982 ,  0.39823586, -0.65425623,  0.9073741 , -1.25465456,
        0.56644717, -0.18324247,  0.42634894,  0.15205491, -0.08515778,
        0.41454226, -0.77506747,  0.71302698, -2.59299147,  0.39412461])

In [31]:

    

10 ** 5


    

    
Out[31]:





100000

In [ ]: