Numerical Solutions to Lotka Volterra and Lorens Equations

Lorenz Equations with Python

import matplotlib.pyplot as plt
import numpy as np

vSigma = 10
vBeta = 8/3
vRho = 28


def f1(t,x,y,z):
    f1 = vSigma*y - vSigma*x
    return f1

def f2(t,x,y,z):
    f2 = vRho*x - y - x*z
    return f2

def f3(t,x,y,z):
    f3 = -vBeta*z + x*y
    return f3


print("Sample input: lorenzEquation(0,20,200,5,5,5)")

def lorenzEquation(first,second,N,alpha1,alpha2, alpha3):
    h = (second-first)/N
    t = first

    w1 = alpha1
    w2 = alpha2
    w3 = alpha3

    A = [w1]
    B = [w2]
    C = [w3]
    tempTime = alpha1
    time = [alpha1]
    
    A1 = [w1]
    B1 = [w2]
    C1 = [w3]


    for num in range(1,N):
        k1x = h*f1(t,w1,w2,w3)
        k1y = h*f2(t,w1,w2,w3)
        k1z = h*f3(t,w1,w2,w3)
        
        k2x = h*f1(t + (h/2), w1 + (k1x/2), w2 + (k1y/2), w3 + (k1z/2))
        k2y = h*f2(t + (h/2), w1 + (k1x/2), w2 + (k1y/2), w3 + (k1z/2))
        k2z = h*f3(t + (h/2), w1 + (k1x/2), w2 + (k1y/2), w3 + (k1z/2))
        
        k3x = h*f1(t + (h/2), w1 + (k2x/2), w2 + (k2y/2), w3 + (k2z/2))
        k3y = h*f2(t + (h/2), w1 + (k2x/2), w2 + (k2y/2), w3 + (k2z/2))
        k3z = h*f3(t + (h/2), w1 + (k2x/2), w2 + (k2y/2), w3 + (k2z/2))
        
        
        k4x = h*f1(t + h, w1 + k3x, w2 + k3y, w3 + k3z)
        k4y = h*f2(t + h, w1 + k3x, w2 + k3y, w3 + k3z)
        k4z = h*f3(t + h, w1 + k3x, w2 + k3y, w3 + k3z)
        
        w1 = w1 + (1/6)*(k1x + 2*k2x + 2*k3x + k4x)
        w2 = w2 + (1/6)*(k1y + 2*k2y + 2*k3y + k4y)
        w3 = w3 + (1/6)*(k1z + 2*k2z + 2*k3z + k4z)

        A.append(w1)
        B.append(w2)
        C.append(w3)
        tempTime = alpha1 + num*h
        time.append(tempTime)

    w1 = alpha1 + 0.001
    w2 = alpha2
    w3 = alpha3


    for num in range(1,N):
        k1x = h*f1(t,w1,w2,w3)
        k1y = h*f2(t,w1,w2,w3)
        k1z = h*f3(t,w1,w2,w3)
        
        k2x = h*f1(t + (h/2), w1 + (k1x/2), w2 + (k1y/2), w3 + (k1z/2))
        k2y = h*f2(t + (h/2), w1 + (k1x/2), w2 + (k1y/2), w3 + (k1z/2))
        k2z = h*f3(t + (h/2), w1 + (k1x/2), w2 + (k1y/2), w3 + (k1z/2))
        
        k3x = h*f1(t + (h/2), w1 + (k2x/2), w2 + (k2y/2), w3 + (k2z/2))
        k3y = h*f2(t + (h/2), w1 + (k2x/2), w2 + (k2y/2), w3 + (k2z/2))
        k3z = h*f3(t + (h/2), w1 + (k2x/2), w2 + (k2y/2), w3 + (k2z/2))
        
        
        k4x = h*f1(t + h, w1 + k3x, w2 + k3y, w3 + k3z)
        k4y = h*f2(t + h, w1 + k3x, w2 + k3y, w3 + k3z)
        k4z = h*f3(t + h, w1 + k3x, w2 + k3y, w3 + k3z)
        
        w1 = w1 + (1/6)*(k1x + 2*k2x + 2*k3x + k4x)
        w2 = w2 + (1/6)*(k1y + 2*k2y + 2*k3y + k4y)
        w3 = w3 + (1/6)*(k1z + 2*k2z + 2*k3z + k4z)

        A1.append(w1)
        B1.append(w2)
        C1.append(w3)
    #
    plt.plot(time, A)
    plt.plot(time, A1)
    plt.legend(['Initial condition [5,5,5]', 'Initial condition [5.001,5,5]'], loc='upper left')
    plt.show()

    plt.plot(A, B)
    plt.show()
    
    plt.plot(A, C)
    plt.show()

    return



Lotka Voltera Equations with Python

import matplotlib.pyplot as plt
import numpy as np

a = 1.2
b = 0.6
c = 0.8
d = 0.3

def f1(t,x,y):
    f1 = a*x - b*x*y
    return f1

def f2(t,x,y):
    f2 = -c*y + d*x*y
    return f2


print("Sample input:  predatorPrey(0, 30, 300, 2, 1)")


def predatorPrey(first,second,N,alpha1,alpha2):
    h = (second-first)/N
    t = first
    
    w1 = alpha1
    w2 = alpha2
    
    A = [w1]
    B = [w2]
    tempTime = alpha1
    time = [alpha1]
    
    for num in range(1,N):
        k1x = h*f1(t,w1,w2)
        k1y = h*f2(t,w1,w2)
        
        k2x = h*f1(t + (h/2), w1 + (k1x/2), w2 + (k1y/2))
        k2y = h*f2(t + (h/2), w1 + (k1x/2), w2 + (k1y/2))
        
        k3x = h*f1(t + (h/2), w1 + (k2x/2), w2 + (k2y/2))
        k3y = h*f2(t + (h/2), w1 + (k2x/2), w2 + (k2y/2))
        
        
        k4x = h*f1(t + h, w1 + k3x, w2 + k3y)
        k4y = h*f2(t + h, w1 + k3x, w2 + k3y)
        
        w1 = w1 + (1/6)*(k1x + 2*k2x + 2*k3x + k4x)
        w2 = w2 + (1/6)*(k1y + 2*k2y + 2*k3y + k4y)
        
        A.append(w1)
        B.append(w2)
        tempTime = alpha1 + num*h
        time.append(tempTime)
        
    
    plt.plot(time, A)
    plt.plot(time, B)
    plt.legend(['x, prey', 'y, predator'], loc='upper left')
    ax.grid()
    ax.set_xlabel("Time (h)")
    plt.show()
    
    plt.plot(A, B)
    plt.show()
    
    return