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