t.ipynb

Kernel: Python 2 (SageMath)

import numpy

def rk2a( f, x0, t ):
    n = len( t )
    x = numpy.array( [ x0 ] * n )
    for i in xrange( n - 1 ):
        h = t[i+1] - t[i]
        k1 = h * f( x[i], t[i] ) / 2.0
        x[i+1] = x[i] + h * f( x[i] + k1, t[i] + h / 2.0 )

    return x

rk2a(f,x0,t)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-2-a905ef492582> in <module>()
----> 1 rk2a(f,x0,t)

NameError: name 'f' is not defined

def rk2a( f, x0, t ):
    n = len( t )
    x = numpy.array( [ x0 ] * n )
    for i in xrange( n - 1 ):
        h = t[i+1] - t[i]
        k1 = h * f( x[i], t[i] ) / 2.0
        x[i+1] = x[i] + h * f( x[i] + k1, t[i] + h / 2.0 )

    return x

o = rk2a([1,2,3], [4,5,6], [7,8,9])

def rk2 (y , time , dt , derivs) :
    k0 = dt * derivs(y , time)
    k1 = dt * derivs (y + k0 , time + dt)
    y.next = y + 0.5 * (k0 + k1)
    return y.next

o = rk2(3, 1000, .08, 2)

from math import sqrt
 
def rk4(f, x0, y0, x1, n):
    vx = [0]*(n + 1)
    vy = [0]*(n + 1)
    h = (x1 - x0)/n
    vx[0] = x = x0
    vy[0] = y = y0
    for i in range(1, n + 1):
        k1 = h*f(x, y)
        k2 = h*f(x + 0.5*h, y + 0.5*k1)
        k3 = h*f(x + 0.5*h, y + 0.5*k2)
        k4 = h*f(x + h, y + k3)
        vx[i] = x = x0 + i*h
        vy[i] = y = y + (k1 + k2 + k2 + k3 + k3 + k4)/6
    return vx, vy
 
def f(x, y):
    return x*sqrt(y)
 
vx, vy = rk4(f, 0, 1, 10, 100)
for x, y in list(zip(vx, vy))[::10]:
    print(x, y, y - (4 + x*x)**2/16)

def f(t, u):
    return -u

def euler(u, v, t, dt):
    v_new = v + dt * f(t, u)
    u_new = u + dt * v
    return u_new, v_new

def rk(u, v, t, dt):
    k1 = f(t,u)
    k2 = f(t+dt*0.5,u+k1*0.5*dt)
    k3 = f(t+dt*0.5,u+k2*0.5*dt)
    k4 = f(t+dt,u+k3*dt)

    v += dt * (k1+2*k2+2*k3+k4)/6

    # v doesn't explicitly depend on other variables
    k1 = k2 = k3 = k4 = v

    u += dt * (k1+2*k2+2*k3+k4)/6

    return k4

o1 = euler(5,12,1000,.08)
o1

o2 = rk(5,12,1000,.08)
o2

def Euler(f, xa, xb, ya, n):
    h = (xb - xa) / float(n)
    x = xa
    y = ya
    for i in range(n):
        y += h * f(x, y)
        x += h
    return y

o3 = Euler(lambda x,y: math.cos(x) + math.sin(y),1,1,1000)

def euler(f,y0,a,b,h):
    t,y = a,y0
    while t <= b:
        print "%6.3f %6.3f" % (t,y)
        t += h
        y += h * f(t,y)
        
def newtoncooling(time, temp):
    return -0.07 * (temp - 20)
 
euler(newtoncooling,100,0,100,10)

import numpy as np

def f(t, u):
    return -u

#def euler(u, v, t, dt):
def euler(arr = np.array([x,y,z], dt):
    #arr = np.array([x,y,z])
    self.arr = arr
    x = arr[0]
    y = arr[1]
    z = arr[2]
    y = x + dt * f(z, x)
    x = x + dt * y
    #v_new = v + dt * f(t, u)
    #u_new = u + dt * v
    return x, y

def rk(u, v, t, dt):
    k1 = f(t,u)
    k2 = f(t+dt*0.5,u+k1*0.5*dt)
    k3 = f(t+dt*0.5,u+k2*0.5*dt)
    k4 = f(t+dt,u+k3*dt)

    v += dt * (k1+2*k2+2*k3+k4)/6

    # v doesn't explicitly depend on other variables
    k1 = k2 = k3 = k4 = v

    u += dt * (k1+2*k2+2*k3+k4)/6

    return k4

import math
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
import bokeh.plotting as blt
blt.output_notebook()

class Attractor(object):
    def __init__(self, s = 10.0, p = 28.0, b = 8.0/3.0, start = 0.0, end = 80.0, points = 10000):
        inarr = np.array([s,p,b])
        self.s = s
        self.p = p
        self.b = b
        self.params = inarr
        self.start = start
        self.end = end
        self.points = points
        self.dt = ((self.end - self.start) / self.points)
        self.t = np.linspace(self.start, self.end, self.points)
        
    #def f1(self, gs, t, xco, yco):
        #return (gs * (yco - xco))
    
    #def f2(self, gp, t, xco, yco, zco):
        #return (xco * (gp - zco) - yco)
    
    #def f3(self, gb, t, xco, yco, zco):
        #return (xco * yco - gb * zco)
    
    def given(self, arr):
        
        x0,y0,z0 = arr
        s,p,b = self.params
        x = s * (y0 - x0)
        y = x0 * (p - z0) - y0
        z = x0 * y0 - b * z0
        
        return np.array([x,y,z])
    
    #def euler(self, r = np.array([])):
    def euler(self, arr):
        
        k1 = arr + self.given(arr) * self.dt
        
        return k1
        #dx = self.f1(self.s, self.start, arr[0], arr[1]) * self.dt
        #dy = self.f2(self.p, self.start, arr[0], arr[1], arr[2]) * self.dt
        #dz = self.f3(self.b, self.start, arr[0], arr[1], arr[2]) * self.dt
        
        #x,y,z = r
        #kx = self.params[0] * (y - x)
        #ky = x * (self.params[1] - z) - y
        #kz = x * y - self.params[2] * z
        
        #print(dx, dy, dz)
        
    def rk2(self, arr):
        
        k1f = self.given(arr)
        k2f = self.given(arr + k1f * self.dt / 2.0)
        
        k2 = arr + k2f * self.dt
        
        return k2
    
    
        #k1x = self.euler = self.kx * self.dt
        #k2x = self.euler = self.ky * self.dt
        #k3x = self.euler = self.kz * self.dt
              
              
        #k1x = self.f1(self.s, self.start, arr[0], arr[1]) * self.dt
        #k1y = self.f2(self.p, self.start, arr[0], arr[1], arr[2]) * self.dt
        #k1z = self.f3(self.b, self.start, arr[0], arr[1], arr[2]) * self.dt
        
        #k2x = self.f1(self.s, self.start + 0.5 * self.dt, arr[0] + 0.5 * self.dt * k1x, arr[1] + 0.5 * self.dt * k1y) * self.dt
        #k2y = self.f2(self.p, self.start + 0.5 * self.dt, arr[0] + 0.5 * self.dt * k1x, arr[1] + 0.5 * self.dt * k1y, arr[2] + 0.5 * self.dt * k1z) * self.dt
        #k2z = self.f3(self.b, self.start + 0.5 * self.dt, arr[0] + 0.5 * self.dt * k1x, arr[1] + 0.5 * self.dt * k1y, arr[2] + 0.5 * self.dt * k1z) * self.dt
        
        #dx = self.dt * k2x
        #dy = self.dt * k2y
        #dz = self.dt * k2z
        
        #print(dx, dy, dz)
        
    def rk4(self, arr):
        
        k1f = self.given(arr)
        k2f = self.given(arr + k1f * self.dt / 2.0)
        k3f = self.given(arr + k2f * self.dt / 2.0)
        k4f = self.given(arr + k3f * self.dt)
        
        k4 = arr + self.dt / 6.0 * (k1f + 2 * k2f + 2 * k3f + k4f)
        
        return k4
                                     
        
        #k1x = self.f1(self.s, self.start, arr[0], arr[1]) * self.dt
        #k1y = self.f2(self.p, self.start, arr[0], arr[1], arr[2]) * self.dt
        #k1z = self.f3(self.b, self.start, arr[0], arr[1], arr[2]) * self.dt
        
        #k2x = self.f1(self.s, self.start + 0.5 * self.dt, arr[0] + 0.5 * self.dt * k1x, arr[1] + 0.5 * self.dt * k1y) * self.dt
        #k2y = self.f2(self.p, self.start + 0.5 * self.dt, arr[0] + 0.5 * self.dt * k1x, arr[1] + 0.5 * self.dt * k1y, arr[2] + 0.5 * self.dt * k1z) * self.dt
        #k2z = self.f3(self.b, self.start + 0.5 * self.dt, arr[0] + 0.5 * self.dt * k1x, arr[1] + 0.5 * self.dt * k1y, arr[2] + 0.5 * self.dt * k1z) * self.dt
        
        #k3x = self.f1(self.s, self.start + 0.5 * self.dt, arr[0] + 0.5 * self.dt * k2x, arr[1] + 0.5 * self.dt * k2y) * self.dt
        #k3y = self.f2(self.p, self.start + 0.5 * self.dt, arr[0] + 0.5 * self.dt * k2x, arr[1] + 0.5 * self.dt * k2y, arr[2] + 0.5 * self.dt * k2z) * self.dt
        #k3z = self.f3(self.b, self.start + 0.5 * self.dt, arr[0] + 0.5 * self.dt * k2x, arr[1] + 0.5 * self.dt * k2y, arr[2] + 0.5 * self.dt * k2z) * self.dt
        
        #k4x = self.f1(self.s, self.start + self.dt, arr[0] + self.dt * k3x, arr[1] + self.dt * k3y)
        #k4y = self.f2(self.p, self.start + self.dt, arr[0] + self.dt * k3x, arr[1] + self.dt * k3y, arr[2] + self.dt * k3z)
        #k4z = self.f3(self.b, self.start + self.dt, arr[0] + self.dt * k3x, arr[1] + self.dt * k3y, arr[2] + self.dt * k3z)
        
        #dx = (self.dt / 6) * (k1x + 2 * k2x + 2 * k3x + k4x)
        #dy = (self.dt / 6) * (k1y + 2 * k2y + 2 * k3y + k4y)
        #dz = (self.dt / 6) * (k1z + 2 * k2z + 2 * k3z + k4z)
        
        #print(dx, dy, dz)
    """"
    def evolve(self, r0, order = 4):
        x0 = 0.1
        y0 = 0.0
        z0 = 0.0
        r0 = np.array([x0,y0,z0], dtype=float)
        #self.x0 = x
        #self.y0 = y
        #self.z0 = z
        for order in xrange(4):
            if order == 1:
                self.euler(r0)
                #self.solution = pd.meshgrid[t,x,y,z]
            elif order == 2:
                self.rk2(r0)
                #self.solution = pd.meshgrid[t,x,y,z]
            elif order == 4:
                self.rk4(r0)
                #self.solution = pd.meshgrid[t,x,y,z]
            else:
                break
        #return self.solution
        """
    def evolve(self, r0 = np.array([0.1, 0.0, 0.0]), order = 4):
        """print("Please enter a 1, 2, or 4.")"""
        if order == 1:
            a = self.euler(r0)
            #self.solution = pd.meshgrid[t, x, y, z]
            return a
        elif order == 2:
            a = self.rk2(r0)
            #self.solution = pd.meshgrid[t, x, y, z]
            return a
        elif order == 4:
            a = self.rk4(r0)
            #self.solution = pd.meshgrid[t, x, y, z]
            return a
        
    def save(self)
    
    def plotx()
    
    def ploty()
    
    def plotz()
    
    def plotxy()
    
    def plotyz()
    
    def plotzx()
    
    def plot3d()

obj = Attractor()
obj.euler([1,0,1])

array([ 0.92      ,  0.216     ,  0.97866667])

obj = Attractor()
obj.rk2([1,0,1])

array([ 0.93184   ,  0.20657792,  0.97972366])

obj = Attractor()
obj.rk4([1,0,1])

array([ 0.93130186,  0.20741988,  0.97968745])

obj = Attractor()
obj.evolve([1,0,1],1)

array([ 0.92      ,  0.216     ,  0.97866667])

obj = Attractor()
obj.evolve([1,0,1],2)

array([ 0.93184   ,  0.20657792,  0.97972366])

obj = Attractor()
obj.evolve([1,0,1],4)

array([ 0.93130186,  0.20741988,  0.97968745])