Infinite loop problem in Newton-Raphson program

[goldenllama]

I wrote this program to calculate f(x) and I'm running into a problem that is causing an infinite loop. The output is fine until the first time that the if-loop catches f_x <= 0.0. I've been looking endlessly for what might be causing it to no avail. Maybe fresher eyes can help. Anyone?

#include <stdio.h>
#include <math.h>
double func1(double x);
double func1_d(double x);
double Newton_Raphson(double x0);
main(void)
{
        double x, x_begin=-1.0, del_x=0.25, x_old, x0, root, f_x, f_x_old;
        int k;
        char sign_change;
        printf("----------------------------------\n");
        printf("      x         f(x)   sign change\n");
        printf("----------------------------------\n");
        x = x_begin;
        f_x = func1(x);
                printf("%8.2f %12.4f\n",x, f_x);
        for (k=1; k<25; k++) {
                x_old = x;
                f_x_old = f_x;
                sign_change = ' ';
               x = x_begin + (double)k * del_x;
                f_x = func1(x);
                if(f_x*f_x_old <= 0.0) {
                        sign_change = 'Y';
                        printf("%8.2f %12.4f    %c\n",x ,f_x ,sign_change);
                        x0 = 0.5 * (x + x_old);
                        root = Newton_Raphson(x0);
                        printf("     A refined root is %-12.6e\n",root);
                }
                else {
                        printf("%8.2f %12.4f    %c\n",x ,f_x ,sign_change);
                }
        }
        printf("\n");
 
        exit(0);/* normal termination */
}
double func1(double x)
{
/* f(x)= x*x*x - x*x -9.0*x + 8.9 */
        double f_x;
        f_x = x * x * x - x * x - 9.0 * x + 8.9;
        return f_x;
}
}
double func1_d(double x)
{
/* f'(x) = 3.0*x*x -2.0 * x - 9.0 */
        double fd_x;
        fd_x + 3.0 * x * x - 2.0 * x - 9.0;
        return fd_x;
}
double Newton_Raphson(double x0)
{
        int debug = 1;
        double tolx, tolf, x1, del_x;
        double f0, f1, f_d0;
        f0 = func1(x0);
        if(debug != 0) printf("     f(%g) = %e \n", x0, f0);
/* define tolerances */
        tolx = 1.e-8 * fabs(x0); /* tolerance for |x1 - x0| */
        tolf = 1.e-6 * fabs(f0); /* tolerance for |f(x1)| */
        do{
                f_d0 = func1_d(x0);
                x1 = x0 - f0/f_d0;
                f1 = func1(x1);
 
                if(debug != 0) printf("     f(%g) = %e\n", x1, f1);
                del_x = fabs(x1 - x0);
/* update x0 and f0 for the next iteration */
                x0 = x1;
                f0 = f1;
        } while(del_x > tolx && fabs(f1) > tolf);
        return x1;
}

Much, much thanks in advance.

[goldenllama]

Some more information. Here is my output:


----------------------------------
x f(x) sign change
----------------------------------
-1.00 15.9000
-0.75 14.6656
-0.50 13.0250
-0.25 11.0719
0.00 8.9000
0.25 6.6031
0.50 4.2750
0.75 2.0094
1.00 -0.1000 Y
f(0.875) = 9.292969e-01
f(-0.125) = 1.000742e+01
f(-1.125) = 1.633555e+01
f(-2.125) = 1.391367e+01
f(-3.125) = -3.258203e+00
f(-4.125) = -4.118008e+01
f(-5.125) = -1.058520e+02
...

When it should be:

----------------------------------
x f(x) sign change
----------------------------------
-1.00 15.9000
-0.75 14.6656
-0.50 13.0250
-0.25 11.0719
0.00 8.9000
0.25 6.6031
0.50 4.2750
0.75 2.0094
1.00 -0.1000 Y
f(0.875) = 9.292969e-01
f(0.984935) = 2.096803e-02
f(0.987537) = 1.324866e-05
f(0.987539) = 5.316636e-12
A refined root is 9.875386e-01
1.25 -1.9594
1.50 -3.4750
1.75 -4.5531
.... .......

So it looks like the code is working as it should, but something is making x1 in Newton_Raphson calculate incorrectly. Still not sure what it is, though.

[aziarani]

You are most probably runing into Oscilation problem. It happens when the solution fluctuate around a local maximum or minimum.
Try some other initial values for your independent variable. This might resolve the issue but there is no guarantee.

[Ganon11]

This code

double func1_d(double x)
{
/* f'(x) = 3.0*x*x -2.0 * x - 9.0 */
        double fd_x;
        fd_x + 3.0 * x * x - 2.0 * x - 9.0;
        return fd_x;
}

isn't doing anything. You aren't assigning any values to fd_x, just using it in an addition expression.