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?
Much, much thanks in advance.
Code:
#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;
}
Comment