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