5 Replies - 7214 Views - Last Post: 30 October 2006 - 09:09 AM Rate Topic: -----

#1 goldenllama  Icon User is offline

  • New D.I.C Head

Reputation: 0
  • View blog
  • Posts: 8
  • Joined: 21-October 06

Infinite loop problem in Newton Raphson program

Posted 29 October 2006 - 06:20 PM

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.

Is This A Good Question/Topic? 0
  • +

Replies To: Infinite loop problem in Newton Raphson program

#2 Amadeus  Icon User is offline

  • g+ + -o drink whiskey.cpp
  • member icon

Reputation: 248
  • View blog
  • Posts: 13,506
  • Joined: 12-July 02

Re: Infinite loop problem in Newton Raphson program

Posted 29 October 2006 - 06:23 PM

Which loop is giving you the problem...if the f_x <= 0.0 catches, then is it the loop in the function that is infinite?
Was This Post Helpful? 0
  • +
  • -

#3 goldenllama  Icon User is offline

  • New D.I.C Head

Reputation: 0
  • View blog
  • Posts: 8
  • Joined: 21-October 06

Re: Infinite loop problem in Newton Raphson program

Posted 29 October 2006 - 09:54 PM

I'm fairly sure the infinite loop is in double Newton_Raphson(double x0). If I comment out root = Newton_Raphson(x0) in main(), the infinite loop doesn't come.
Was This Post Helpful? 0
  • +
  • -

#4 goldenllama  Icon User is offline

  • New D.I.C Head

Reputation: 0
  • View blog
  • Posts: 8
  • Joined: 21-October 06

Re: Infinite loop problem in Newton Raphson program

Posted 30 October 2006 - 07:56 AM

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.
Was This Post Helpful? 0
  • +
  • -

#5 horace  Icon User is offline

  • D.I.C Lover
  • member icon

Reputation: 291
  • View blog
  • Posts: 1,900
  • Joined: 25-October 06

Re: Infinite loop problem in Newton Raphson program

Posted 30 October 2006 - 08:17 AM

one problem is that you have a + where you should have an = in
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;
}



i.e.
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;   //** change = for +
		return fd_x;
}



Was This Post Helpful? 0
  • +
  • -

#6 goldenllama  Icon User is offline

  • New D.I.C Head

Reputation: 0
  • View blog
  • Posts: 8
  • Joined: 21-October 06

Re: Infinite loop problem in Newton Raphson program

Posted 30 October 2006 - 09:09 AM

That'll do it. Goes to show that the hardest problems are always the easiest ones. Thanks so much.
Was This Post Helpful? 0
  • +
  • -

Page 1 of 1