Monday, May 29, 2017

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


Read More

Saturday, May 27, 2017

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


Read More

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



Read More

Thursday, May 25, 2017

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
Read More

Tuesday, May 16, 2017

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
Read More

Saturday, May 13, 2017

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


 

 

 

 

 

 
Read More