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

Regula Falsi Or Method of False Position with Python

Regula Falsi or Method of False Position

     The regula falsi method iteratively determines a sequence of root enclosing intervals, $(a_n, b_n)$, and a sequence of approximations, which shall be denoted by $p_n$. Similar to the bisection method, the root should be in ther interval being considered. During each iteration, a single point is selected from $(a_n, b_n)$ to approximate the location of the root and serve as $p_n$. If $p_n$ is an accurate enough approximation, the iterative process is terminated. Otherwise, the Intermediate Value Theorem is used to determine whether the root lies on the subinterval $(a_n, p_n)$ or the subinterval $(p_n, b_n)$. The entire process is then repeated on that subinterval. It was developed because the Bisection method converges at a fairly slow rate.

Let f be a continuous function on the interval $[a,b]$ s.t. $f(a) \cdot f(b) < 0$, locate the point $(p1,0)$ where the line joining the points $(a, f(a))$ and $(b,f(b))$ crosses the x-axis. Hence,
       $$p_1 = b -  \frac{f(b)(b-a)}{f(b)-f(a)} = \frac{af(b) - bf(a)}{f(b) - f(a)}$$

Algorithm

To find a solution to $f(x) = 0$ given the continuous function $f$ on the interval $[a, b]$, where $f(a)$ and $f(b)$ have opposite signs:

INPUT endpoints a, b; tolerance TOL; maximum number of iterations $N_0$.

STEP 1 Set $i = 1$
                     $FA = f(a)$.

STEP 2 While $i \le N_0$ do Steps 3-6.

        STEP 3 Set $p = \frac{af(b) - bf(a)}{f(b) - f(a)}$
                               $FP = f(p)$

        STEP 4 If $FP = 0$ or |f(p)| < TOL  then
                         STOP
                     else OUTPUT(P)
                 
        STEP 5 Set $i = i + 1$

       STEP 6 If $FA \times FP > 0$ then set $a = p$;
                          $FA = FP$
                     else set $b = p$.

STEP 7 OUTPUT("Method failed after $N_0$")

Sample Problem:

Use Regula Falsi method to approximate the solution of $f(x) = x^3 + 2x^2 - 3x - 1 = 0$ within $[1, 2]$ that is accurate to at least within $10^-4$.


For the approximation, see the outpout below:

    n                    $a_n$                                $b_n$                        $p_n$                                  $f(p_n)$
         
   1                     1                                   2                       1.1                                  -0.549            

   2                    1.1                                 2                       1.1517436                      -0.27440072      

   3                    1.1517436                     2                       1.1768409                      -0.13074253      

   4                    1.1768409                     2                       1.1886277                      -0.060875863      

  5                    1.1886277                      2                       1.1940789                      -0.028040938      

  6                    1.1940789                      2                       1.1965821                      -0.01285224      

  7                    1.1965821                      2                       1.1977278                       -0.0058772415    

  8                    1.1977278                      2                       1.1982513                       -0.0026848163    

  9                    1.1982513                      2                       1.1984904                       -0.001225881      

 10                   1.1984904                     2                        1.1985996                       -0.0005596125    

 11                   1.1985996                     2                        1.1986494                        -0.00025543669    

 12                   1.1986494                     2                        1.1986721                        -0.0001165895    

Python Code:

import math
import numpy as np



def f(x):
    f = math.pow(x,3) + 2*math.pow(x,2) - 3*x - 1
    return f
 
 
print("Sample input: regulaFalsi(1,2,10**-4, 100)")
 
def regulaFalsi(a,b,TOL,N):
    i = 1
    FA = f(a)
    
    print("%-20s %-20s %-20s %-20s %-20s" % ("n","a_n","b_n","p_n","f(p_n)"))
     
    while(i <= N):
        p = (a*f(b)-b*f(a))/(f(b) - f(a))
        FP = f(p)
         
        if(FP == 0 or np.abs(f(p)) < TOL):
            break
        else:
             print("%-20.8g %-20.8g %-20.8g %-20.8g %-20.8g\n" % (i, a, b, p, f(p)))
        
         
        i = i + 1
         
        if(FA*FP > 0):
            a = p
        else:
            b = p
     
    return

Regula Falsi or Method of False Position with Scilab

Regula Falsi or Method of False Position

     The regula falsi method iteratively determines a sequence of root enclosing intervals, $(a_n, b_n)$, and a sequence of approximations, which shall be denoted by $p_n$. Similar to the bisection method, the root should be in ther interval being considered. During each iteration, a single point is selected from $(a_n, b_n)$ to approximate the location of the root and serve as $p_n$. If $p_n$ is an accurate enough approximation, the iterative process is terminated. Otherwise, the Intermediate Value Theorem is used to determine whether the root lies on the subinterval $(a_n, p_n)$ or the subinterval $(p_n, b_n)$. The entire process is then repeated on that subinterval. It was developed because the Bisection method converges at a fairly slow rate.

Let f be a continuous function on the interval $[a,b]$ s.t. $f(a) \cdot f(b) < 0$, locate the point $(p1,0)$ where the line joining the points $(a, f(a))$ and $(b,f(b))$ crosses the x-axis. Hence,
       $$p_1 = b -  \frac{f(b)(b-a)}{f(b)-f(a)} = \frac{af(b) - bf(a)}{f(b) - f(a)}$$

Algorithm

To find a solution to $f(x) = 0$ given the continuous function $f$ on the interval $[a, b]$, where $f(a)$ and $f(b)$ have opposite signs:

INPUT endpoints a, b; tolerance TOL; maximum number of iterations $N_0$.

STEP 1 Set $i = 1$
                     $FA = f(a)$.

STEP 2 While $i \le N_0$ do Steps 3-6.

        STEP 3 Set $p = \frac{af(b) - bf(a)}{f(b) - f(a)}$
                               $FP = f(p)$

        STEP 4 If $FP = 0$ or |f(p)| < TOL  then
                         STOP
                     else OUTPUT(P)
                 
        STEP 5 Set $i = i + 1$

       STEP 6 If $FA \times FP > 0$ then set $a = p$;
                          $FA = FP$
                     else set $b = p$.

STEP 7 OUTPUT("Method failed after $N_0$")

Sample Problem:

Use Regula Falsi method to approximate the solution of $f(x) = x^3 + 2x^2 - 3x - 1 = 0$ within $[1, 2]$ that is accurate to at least within $10^-4$.


For the approximation, see the outpout below:

    n                    $a_n$                              $b_n$                    $p_n$                        $f(p_n)$
         
    1                    1                                 2                    1.1                        -0.549            
    2                    1.1                              2                    1.1517436            -0.27440072      
    3                    1.1517436                  2                    1.1768409            -0.13074253      
    4                    1.1768409                  2                    1.1886277            -0.060875863      
    5                    1.1886277                  2                    1.1940789            -0.028040938      
    6                    1.1940789                  2                    1.1965821            -0.01285224      
    7                    1.1965821                  2                    1.1977278            -0.0058772415    
    8                    1.1977278                  2                    1.1982513            -0.0026848163    
    9                    1.1982513                  2                    1.1984904            -0.001225881      
   10                   1.1984904                  2                    1.1985996            -0.0005596125    
   11                   1.1985996                  2                    1.1986494            -0.00025543669    
   12                   1.1986494                  2                    1.1986721            -0.0001165895

Scilab Code:

clear
clc
 
function f = f(x)
    f = x^3 + 2*x^2 - 3*x -1 
endfunction
 
 
disp("Sample input: regulaFalsi(1,2,10^-4, 100)")
 
function regulaFalsi(a,b,TOL,N)
    i = 1
    FA = f(a)
    finalOutput = [i, a, b, a + (b-a)/2, f(a + (b-a)/2)]
     
    printf("%-20s %-20s %-20s %-20s %-20s \n","n","a_n","b_n","p_n","f(p_n)")
    
    while(i <= N),
        p = (a*f(b)-b*f(a))/(f(b) - f(a))
        FP = f(p)
         
         
        if(FP == 0 | abs(f(p)) < TOL) then
            break
        else
             printf("%-20.8g %-20.8g %-20.8g %-20.8g %-20.8g\n", i, a, b, p, f(p))
        end
         
        i = i + 1
         
        if(FA*FP > 0) then
            a = p
        else
            b = p
        end
    end
     
    //disp(finalOutput)
     
endfunction

Bisection Method with Python

The Bisection Method

     Suppose $f$ is a continuous function defined on the interval $[a, b]$, with $f(a)$ and $f(b)$ of opposite sign. The Intermediate Value Theorem implies that a number p exists in (a, b) with $f( p) = 0$. Although the procedure will work when there is more than one root in the interval $(a, b)$, we assume for simplicity that the root in this interval is unique. The method calls for a repeated halving (or bisecting) of subintervals of $[a, b]$ and, at each step, locating the half containing p.

Algorithm

To find a solution to $f(x) = 0$ given the continuous function $f$ on the interval $[a, b]$, where $f(a)$ and $f(b)$ have opposite signs:

INPUT endpoints a, b; tolerance TOL; maximum number of iterations $N_0$.

STEP 1 Set $i = 1$
                     $FA = f(a)$.

STEP 2 While $i \le N_0$ do Steps 3-6.

        STEP 3 Set $p = a + (b - a)/2$
                               $FP = f(p)$

        STEP 4 If $FP = 0$ or $(b-a)/2 < TOL$ then
                         STOP
                     else OUTPUT(P)
                   
        STEP 5 Set $i = i + 1$

       STEP 6 If $FA \times FP > 0$ then set $a = p$;
                          $FA = FP$
                     else set $b = p$.

STEP 7 OUTPUT("Method failed after $N_0$")

Sample Problem:

Show that $f(x) = x^3 + 4x^2 - 10 = 0$ has a root in $[1, 2]$, and use the Bisection method to determine an approximation to the root that is accurate to at least within $10^-4$.

Solution: Because $f(1) = -5$ and $f(2) = 14$ the Intermediate Value Theorem ensures that this continuous function has a root in $[1, 2]$.

For the approximation, see the outpout below:

    n                    $a_n$                        $b_n$                          $p_n$                        $f(p_n)$
           
    1                    1                           2                           1.5                       2.375                    

    2                    1                           1.5                        1.25                     -1.796875        

    3                    1.25                      1.5                        1.375                   0.16210938        

    4                    1.25                      1.375                    1.3125                 -0.84838867      

    5                    1.3125                  1.375                    1.34375               -0.35098267      

    6                    1.34375                1.375                    1.359375             -0.096408844      

    7                    1.359375              1.375                    1.3671875            0.032355785      

    8                    1.359375              1.3671875            1.3632812            -0.032149971      

    9                    1.3632812            1.3671875            1.3652344            7.2024763e-05    

    10                  1.3632812            1.3652344            1.3642578            -0.016046691      

    11                  1.3642578            1.3652344            1.3647461            -0.0079892628    

    12                  1.3647461            1.3652344            1.3649902            -0.0039591015    

    13                  1.3649902            1.3652344            1.3651123            -0.001943659      

Python Code:

def f(x):
    f = x**3 + 4*x**2 - 10
    return f
 
 
print("Sample input: bisectionMethod(1,2,10**-4, 100)")
 
def bisectionMethod(a,b,TOL,N):
    i = 1
    FA = f(a)
    
    print("%-20s %-20s %-20s %-20s %-20s" % ("n","a_n","b_n","p_n","f(p_n)"))
    print("%-20.8g %-20.8g %-20.8g %-20.8g %-20.8g\n" % (i, a, b, a + (b-a)/2, f(a + (b-a)/2) ))
    
     
    while(i <= N):
        p = a + (b-a)/2
        FP = f(p)
         
        if(FP == 0 or (b-a)/2 < TOL):
            break
        else:
             print("%-20.8g %-20.8g %-20.8g %-20.8g %-20.8g\n" % (i, a, b, p, f(p)))
        
         
        i = i + 1
         
        if(FA*FP > 0):
            a = p
        else:
            b = p
     
    return

Final Note: 

The Bisection method, though conceptually clear, has significant drawbacks. It is relatively
slow to converge (that is, N may become quite large before $| p − p_N|$ is sufficiently
small), and a good intermediate approximation might be inadvertently discarded. However,
the method has the important property that it always converges to a solution, and for that
reason it is often used as a starter for the more efficient methods

Bisection Method with Scilab

The Bisection Method

     Suppose $f$ is a continuous function defined on the interval $[a, b]$, with $f(a)$ and $f(b)$ of opposite sign. The Intermediate Value Theorem implies that a number p exists in (a, b) with $f( p) = 0$. Although the procedure will work when there is more than one root in the interval $(a, b)$, we assume for simplicity that the root in this interval is unique. The method calls for a repeated halving (or bisecting) of subintervals of $[a, b]$ and, at each step, locating the half containing p.

Algorithm

To find a solution to $f(x) = 0$ given the continuous function $f$ on the interval $[a, b]$, where $f(a)$ and $f(b)$ have opposite signs:

INPUT endpoints a, b; tolerance TOL; maximum number of iterations $N_0$.

STEP 1 Set $i = 1$
                     $FA = f(a)$.

STEP 2 While $i \le N_0$ do Steps 3-6.

        STEP 3 Set $p = a + (b - a)/2$
                               $FP = f(p)$

        STEP 4 If $FP = 0$ or $(b-a)/2 < TOL$ then
                         STOP
                     else OUTPUT(P)
                   
        STEP 5 Set $i = i + 1$

       STEP 6 If $FA \times FP > 0$ then set $a = p$;
                          $FA = FP$
                     else set $b = p$.

STEP 7 Display("Method failed after $N_0$")

Sample Problem:

Show that $f(x) = x^3 + 4x^2 - 10 = 0$ has a root in $[1, 2]$, and use the Bisection method to determine an approximation to the root that is accurate to at least within $10^-4$.

Solution: Because $f(1) = -5$ and $f(2) = 14$ the Intermediate Value Theorem ensures that this continuous function has a root in $[1, 2]$.

For the approximation, see the outpout below:

    n       $a_n$                   $b_n$                $p_n$                   $f(p_n)$

    1.      1.                   2.                  1.5                  2.375    
    1.      1.                   2.                  1.5                  2.375    
    2.      1.                   1.5                1.25                - 1.796875
    3.      1.25               1.5                1.375              0.1621094
    4.      1.25               1.375            1.3125            - 0.8483887
    5.      1.3125           1.375            1.34375          - 0.3509827
    6.      1.34375         1.375            1.359375        - 0.0964088
    7.      1.359375       1.375            1.3671875      0.0323558
    8.      1.359375       1.3671875    1.3632812      - 0.0321500
    9.      1.3632812     1.3671875    1.3652344      0.0000720
    10.    1.3632812     1.3652344    1.3642578      - 0.0160467
    11.    1.3642578     1.3652344    1.3647461      - 0.0079893
    12.    1.3647461     1.3652344    1.3649902      - 0.0039591
    13.    1.3649902     1.3652344    1.3651123      - 0.0019437

Scilab Code:

clear
clc

function f = f(x)
    f = x^3 + 4*x^2 - 10
endfunction


disp("Sample input: bisectionMethod(1,2,10^-4, 100)")

function bisectionMethod(a,b,TOL,N)
    i = 1
    FA = f(a)
    finalOutput = [i, a, b, a + (b-a)/2, f(a + (b-a)/2)]
    
    disp("   n      a_n          b_n          p_n         f(p_n)")
    
    while(i <= N), 
        p = a + (b-a)/2
        FP = f(p)
        
        if(FP == 0 | (b-a)/2 < TOL) then
            break
        else
             finalOutput = [finalOutput; i, a, b, p, f(p)]
        end
        
        i = i + 1
        
        if(FA*FP > 0) then
            a = p
        else
            b = p
        end
    end
    
    disp(finalOutput)
    
endfunction

Final Note: 

The Bisection method, though conceptually clear, has significant drawbacks. It is relatively
slow to converge (that is, N may become quite large before $| p − p_N|$ is sufficiently
small), and a good intermediate approximation might be inadvertently discarded. However,
the method has the important property that it always converges to a solution, and for that
reason it is often used as a starter for the more efficient methods

Integrating PrimeNG with Angular CLI

1.) Install the dependency, and save it to your project:

 npm install primeng --save
 npm install font-awesome --save

Once that is installed, you’ll have a ./node_modules/primeng/ directory (and a ./node_modules/font-awesome/ directory)

2.) In your project folder there is a file called styles.css, add the following lines:

 @import url('../node_modules/primeng/resources/themes/omega/theme.css');
 @import url('../node_modules/primeng/resources/primeng.min.css');
 @import url('../node_modules/font-awesome/css/font-awesome.min.css');

3.) Install the animations

 npm install @angular/animations

With the animation installed, update the app.module.ts to load the module

Replace this line of code: 

 import {BrowserModule} from '@angular/platform-browser';

with

 import { BrowserAnimationsModule } from '@angular/platform-browser/animations';

followed by the imports.

 BrowserAnimationsModule

Your app.module.ts will have something like:

 @NgModule({
 imports: [
      BrowserAnimationsModule,
      //...
 ],
      //...
 })

With the styles in place, it’s time to add the modules you’ll need. For this example, we'll just use the basic P-Button

4.) Add the following line to the top of your app.module.ts

 import {ButtonModule} from 'primeng/primeng';

Then finally, include the component in your imports: section of app.module.ts towards the bottom of the file. Your app.module.ts will have something like:

 @NgModule({
 imports: [
      BrowserAnimationsModule,
      //...,
      ButtonModule
 ],
      //...
 })

With everything in place, you can now start using p-Button tag in your actual markup.

5.) Add the following lines to your app.component.html