Regula-Falsi Algorithm?

Question

I'm trying to implement the Regula-Falsi Algorithm to solve the equation of 2(x^3)-x-2 but the problem is that the variable c value is staying constant and it is not changing, even that my code should change it .

#include<math.h>
#include<stdio.h>


float fonc(float x)
{
    int result;
    result=2*(pow(x,3))-x-2;
    return result;
}

int main(void)
{
    float eps=pow(10,-4);
    int i=0;
    float a,b,c;
    a=1;
    b=2;
    do
    {
        c=((a*fonc(b))-(b*fonc(a)))/((fonc(b)-fonc(a)));
        if(fonc(c)*fonc(a)<0)
        {
            b=c;
        }
        else
        {
            a=c;
        }
        i++;    
        printf("\n%f",c);
    }
    while(fabs(b-c)>eps);

    printf("le nombre d'itération %d",i);
}

Answer 1

Okay, just for the sake of correcting the algorithm to actually perform what it was meant to perform:

Along with changing the type of result from int to float, you also need to change the loop condition. Currently it is:

    while ( fabs( b - c ) > eps );

Which means that the loop will continue happening until the distance between the values b and c become lower than 0.0001, and under this condition, currently, at least on my end, the code runs forever.

We aren't after reducing the difference between b and c anyway. What we really are after is that the fonc( c ) = 2*c*c*c - c - 2 to be less than our eps. In the end, we want it to be as close as possible to zero, so that c becomes approximately a root of the function. So simply:

    while ( fabs( fonc( c ) ) > eps );

is what our condition should be. This way, along with the int --> float change, it doesn't fall into an infinite loop, completing the job in 14 iterations.

Answer 2

What can go wrong with the algorithm as presented, even if all data types were appropriate?

In contrast to the bisection method, the pure regula falsi will not force the interval length to go to zero. If a monotone and convex function is used, the iteration will stall at changing only one side of the interval with geometric convergence. Any sufficiently smooth function will eventually have these properties over the remaining bracketing interval.

To correctly catch this behavior, the midpoint c should be compared to both interval ends a, b immediately after computation. As bonus point, check if the value at c is small enough and if true, also break the iteration regardless of how far from the ends of the interval this is.

---

There exist many simple tricks to force the interval length to be zero. The complicated tricks lead to the Brent's method. Of simple tricks one is the Illinois variant. In these variants, the midpoint is considered as a convex sum

    c = |f(b)|/(|f(a)|+|f(b)|) * a + |f(a)|/(|f(a)|+|f(b)|) * b

Because of the opposite signs of f(a) and f(b), this is equivalent to the original formula. If the side b does not change, its importance in this convex sum is increased by decreasing the function value f(b), i.e., multiplying it with additional weight factors. This moves the midpoint c towards b, which will in very few steps find a midpoint that will replace b.

---

The following is an implementation of the Illinois variant of regula falsi (or false position method). The algorithm finds a solution with function value 2.2e-6 and an enclosing interval of length 6e-7 in 6 iterations.

#include<math.h>
#include<stdio.h>


float fonc(float x)
{
   return (2*x*x-1)*x-2;
}

int main(void)
{
    float eps=1e-6;
    int i=0;
    float a=1, fa = fonc(a);
    float  b=2, fb = fonc(b);
    
    printf("\na=%10.7f b=%10.7f  fa=%10.7f  fb=%10.7f\n------\n",a,b, fa,fb);
    
    if(signbit(fb)==signbit(fa)) {
        printf("Attention, les valeurs initiales de 'fonc' n'ont pas de signe opposeés!\n");
        
    }
    do
    {
        float c=(a*fb-b*fa)/(fb-fa), fc = fonc(c);
        
        if( signbit(fc)!=signbit(fa) )
        {
            b=a; fb=fa;
        }
        else
        {
            fb *= 0.5;
        }
        a=c; fa=fc;
        i++;  
        printf("\na=c=%10.7f b=%10.7f  fa=fc=%10.7f  fb=%10.7f",c,b, fc,fb);
        
        if(fabs(fc)<eps) break;
    }
    while(fabs(b-a)>eps);

    printf("\nle nombre d'itération %d\n",i);
    
    return 0;
}

The output is

     a= 1.0000000 b= 2.0000000  fa=-1.0000000  fb=12.0000000
     ------
     
     a=c= 1.0769231 b= 2.0000000  fa=fc=-0.5789710  fb= 6.0000000
     a=c= 1.1581569 b= 2.0000000  fa=fc=-0.0512219  fb= 3.0000000
     a=c= 1.1722891 b= 1.1581569  fa=fc= 0.0497752  fb=-0.0512219
     a=c= 1.1653242 b= 1.1722891  fa=fc=-0.0003491  fb= 0.0497752
     a=c= 1.1653727 b= 1.1722891  fa=fc=-0.0000022  fb= 0.0248876
     a=c= 1.1653733 b= 1.1653727  fa=fc= 0.0000020  fb=-0.0000022
     le nombre d'itération 6

Answer 3

One problem is that result is int. You almost certainly want to be float or double, as otherwise the result of fonc() is getting truncated to an integer.