/* simple2_nl.c  handle reciprocal terms    */
/*  f(x) = x + 2 x^2 + 3/x - 11.5    thus:  */
/*  f'(x) = 1 + 4 x - 6/x^2                 */
/*  Newton iteration:                       */
/*  x_i+1 = x_i - f(x)/f'(x)                */
#include <stdio.h>

double f(double x)
{
  return x + 2.0*x*x + 3.0/x - 11.5;
}

double fp(double x)
{
  return 1.0 + 4.0*x - 6.0/(x*x);
}

int main(int argc, char * argv[])
{
  double x, x_next, fx, fpx;
  int i;
  int n = 6;

  printf("simple2_nl.c running \n");
  printf("solve x + 2 x^2 + 3/x = 11.5 \n");
  printf("define f(x) = x + 2 x^2 + 3/x - 11.5 \n");
  printf("derivative fp(x) = 1.0 + 4.0*x - 6.0/(x*x) \n");
  printf("guess initial x \n");
  printf("iterate x_next = x - f(x)/fp(x) \n");
  printf("no guarentee of solution or unique solution \n");
  printf("f(x) should go to zero \n");

  x = 3.0;
  for(i=0; i<n; i++)
  {
    fx = f(x);
    fpx = fp(x);
    x_next = x - fx/fpx;
    printf("x=%8.4f, f(x)=%8.4f, fp(x)=%8.4f, f/fp=%8.4f, x_next=%8.4f \n",
           x, fx, fpx, fx/fpx, x_next);
    x = x_next;
  }
  printf("\nnew initial guess\n");
  x = 1.5;
  for(i=0; i<n; i++)
  {
    fx = f(x);
    fpx = fp(x);
    x_next = x - fx/fpx;
    printf("x=%8.4f, f(x)=%8.4f, fp(x)=%8.4f, f/fp=%8.4f, x_next=%8.4f \n",
           x, fx, fpx, fx/fpx, x_next);
    x = x_next;
  }
  printf("\nnew initial guess\n");
  x = 2.1;
  for(i=0; i<n; i++)
  {
    fx = f(x);
    fpx = fp(x);
    x_next = x - fx/fpx;
    printf("x=%8.4f, f(x)=%8.4f, fp(x)=%8.4f, f/fp=%8.4f, x_next=%8.4f \n",
           x, fx, fpx, fx/fpx, x_next);
    x = x_next;
  }
  printf("\nnew initial guess\n");
  x = 2.0;
  for(i=0; i<4; i++)
  {
    fx = f(x);
    fpx = fp(x);
    x_next = x - fx/fpx;
    printf("x=%8.4f, f(x)=%8.4f, fp(x)=%8.4f, f/fp=%8.4f, x_next=%8.4f \n",
           x, fx, fpx, fx/fpx, x_next);
    x = x_next;
  }

  printf("simple2_nl.c running \n");
  return 0;
} /* end simple_nl.c */