Difference between two roots using bisection method in C?

Question

I have written this C code to compare the root of two functions using bisection method. My first function (g(x)) executes correctly but the second one (h(x)) outputs "#1QO" on screen. I can't find what I did wrong in the code.

Can you please explain? Any kind of help would be highly appreciated. Thanks!

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

typedef double (*DFD) (double);

double g (double x) 
{
   double y;
   y = pow (x,3) - pow (x,2)-1;
   return y;
}

double h (double x)
{
  double k;
  k = 1.0 + (1.0/x) + (1.0 /pow (x,2));
  return k;
}

double bisection (DFD f, double x0, double x1,double tol)
{
   int i;
   double middle;
     for (i=1;i<=50;) {
     middle = (x0+x1)/2.0;
     if (fabs (middle - x0) <tol) return middle;
     if (f(middle)* f(x0) <0.0) x1 = middle;
     else x0 = middle;
     }
 }

 int main () {
 double root_gx, root_hx = 0.0;

 root_gx = bisection (g,0,2,0.0005);
 printf ("Root found using g(x) = %.3lf\n",root_gx);

 root_hx = bisection (h,1,2,0.0005);
 printf ("Root found using h(x) = %.3lf\n",root_hx);

 printf ("Difference between the two roots = %.3lf\n", (fabs (root_gx- root_hx)));

return 0;

} 

EDIT: Initialised i=1 in bisection and changed bisection (h,0,2,0.0005) to bisection (h,1,2,0.0005) and it works Thanks everyone!

Answer 1

Value not initialized leading to undefined behavior. Set to some value.

// int i;
int i = 0;
...
for (;i<=50;) {

---

i is never incremented.
50 is arbitrary.

Recommend using the binary precision of double as bisection() halves each iteration.

for (i=0; i <= DBL_MANT_DIG; i++) {
  ...
}
return middle;

---

Also recommend an algorithm change to allow a change of |difference| of 0.0 to meet a tolerance of 0.0.

// if (fabs (middle - x0) <tol) return middle;
if (fabs (middle - x0) <= tol) return middle;

Answer 2

double h (double x)
{
  double k;
  k = 1.0 + (1.0/x) + (1.0 /pow (x,2));
  return k;
}

gives division by zero if x==0.0. You do this when you call bisection (h,0,2,0.0005);, which calls f(x0), which is h(0.0); in this case.

If you divide through an input parameter in a function, then you should always check, if the actual value of the input parameter is not equal to 0.0:

double h (double x)
{
  if ( x == 0.0 ) return 1e99; // infinite
  double k;
  k = 1.0 + (1.0/x) + (1.0 /pow (x,2));
  return k;
}