Jump to content

Using Newton's method to find the root of an equation! Need verification!

- - - - -

  • Please log in to reply
17 replies to this topic

#1
Yuriy M

Yuriy M

    Learning Programmer

  • Members
  • PipPipPip
  • 70 posts
Hi guys. :)

I've been working on a programming problem that makes use of Newton's method to help find the root of an equation. The procedure starts with an initial guess of "c / 2" and c is equal to the nth power of x. The program uses test values consisting of the square root of 2, the cube root of 7 and the cube root of -1. I'm also making use of the iterative formula to help find the root consisting of the new guess of xj + 1 being equal to the current guess of xj minus the power rule function of xj divided by the derivative of xj.

The program makes 100 guesses before quitting and moving on to the next value in the input file. Here is what I have:

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

double power_rule(double, int);                         /* Function that calculates the power rule of the equation */
double derivative(double, int);                         /* Function that calculates the derivative of the equation */

int main(void)
{
    /* Declare variables */
    
    FILE   *inp;                                        /* The input file */
    double curr_guess,                                  /* The current guess value for the root of the equation */
           new_guess,                                   /* The new guess value for the root of the equation */
           c;                                           /* The current value to be used to determine the initial guess */
    int    trial_counter,                               /* Counts the guess trials for determining the root of the equation */
           inp_value,                                   /* The input value from the file */
           inp_counter = 0,                             /* Counts the input values from the file */
           input_status;                                /* The status of the file input */
    
    /* Open file */
    
    inp = fopen("root_test_values.dat", "r");
    input_status = fscanf(inp, "%d", &inp_value);
    
    /* Guess the root */
    
    while (input_status != EOF)
    {
     if (inp_value < 0)
     {
      printf("\n%d is a negative value. Negative values cannot be used for root operations.\n", inp_value);
      inp_counter++;
      input_status = fscanf(inp, "%d", &inp_value);
     }
     else
     {
      if (inp_counter == 0)
       c = pow(inp_value, 0.5);
      else
       c = pow(inp_value, 0.33);
      
      curr_guess = c / 2;
      printf("\nStarting guess for root is %.6f", curr_guess);
     
      for (trial_counter = 0; trial_counter < 100; trial_counter++)
      {
       new_guess = curr_guess - (power_rule(curr_guess, inp_counter) / derivative(curr_guess, inp_counter));
       printf("\nNew guess for root is %.6f", new_guess);
      
       if (inp_counter == 0)
       {
        if (curr_guess == pow(c, 0.5))
        {
         printf("\nRoot found at %.6f", curr_guess);
         break;
        }
        else if (new_guess == pow(c, 0.5))
        {
         printf("\nRoot found at %.6f", new_guess);
         break;
        }
       }
       else
       {
        if (curr_guess == pow(c, 0.33))
        {
         printf("\nRoot found at %.6f", curr_guess);
         break;
        }
        else if (new_guess == pow(c, 0.33))
        {
         printf("\nRoot found at %.6f", new_guess);
         break;
        }
       }
       curr_guess = new_guess;
      } 
     
      if (inp_counter == 0 && trial_counter == 100)
       printf("\nNo root was found for the square root of %d.\n", inp_value);
      if (inp_counter != 0 && trial_counter == 100)
       printf("\nNo root was found for the cube root of %d.\n", inp_value);
      
      inp_counter++; 
      input_status = fscanf(inp, "%d", &inp_value);
     }
    }
    
    /* Close file */
    
    fclose(inp);
    
    return 0;   
}

double power_rule(double x, int inp_count)
{
    double n;                                           /* The denomination of the fractional power */
    
    if (inp_count == 0)
    {
     n = 1 / 0.5;
     return pow(x, n);
    }
    else 
    {
     n = 1 / 0.33;
     return pow(x, n);
    }
}

double derivative(double x, int inp_count)
{
    double n,                                           /* The denomination of the fractional power */ 
           exp_sub;                                     /* The subtracted exponent */
    
    if (inp_count == 0)
    {
     n = 1 / 0.5;
     exp_sub = n - 1;
     return n * pow(x, exp_sub);            
    }
    else
    {
     n = 1 / 0.33;
     exp_sub = n - 1;
     return n * pow(x, exp_sub);  
    }
}

The problem in my textbook seems to indicate that f(x) = nth root of c but I'm not sure. I am also not sure if Newton's method applies to the square or cube root of negative values. Some help would be appreciated. Thanks. :)
For $1000: Something that is a miserable pile of secrets.

#2
WingedPanther

WingedPanther

    A spammer's worst nightmare

  • Moderators
  • 16,831 posts
  • Location:Upstate, South Carolina
  • Programming Language:C, C++, PL/SQL, Delphi/Object Pascal, Pascal, Transact-SQL, Others
  • Learning:Java, C#, PHP, JavaScript, Lisp, Fortran, Haskell, Others
it applies to cube roots of negative values, but not square roots.
Programming is a branch of mathematics.
My CodeCall Blog | My Personal Blog

#3
Yuriy M

Yuriy M

    Learning Programmer

  • Members
  • PipPipPip
  • 70 posts
So the cube root of -1 is still okay for this type of problem?
For $1000: Something that is a miserable pile of secrets.

#4
WingedPanther

WingedPanther

    A spammer's worst nightmare

  • Moderators
  • 16,831 posts
  • Location:Upstate, South Carolina
  • Programming Language:C, C++, PL/SQL, Delphi/Object Pascal, Pascal, Transact-SQL, Others
  • Learning:Java, C#, PHP, JavaScript, Lisp, Fortran, Haskell, Others
Yes. cube root (-1) = -1. With that said, your pow implementation may have issues with it :)
Programming is a branch of mathematics.
My CodeCall Blog | My Personal Blog

#5
Yuriy M

Yuriy M

    Learning Programmer

  • Members
  • PipPipPip
  • 70 posts
Thanks for clarifying that, Panther. :) Now, I hope you can clarify one other thing.

After looking over the requirements for the question, I'm not exactly sure what f(x) is really supposed to be. Here is what the question originally asked:

Textbook said:

In this chapter we studied the bisection method for finding a root of an equation. Another method for finding a root, Newton's method, usually converges to a solution even faster than the bisection method, if it converges at all. Newton's method starts with an initial guess for a root x0, and then generates successive approximate root x1, x2, ..., xj, xj+1, ..., using the iterative formula

xj+1 = xj - f(xj) / f`(xj)

where f`(xj) is the derivative of function f evaluated at x = xj. The formula generates a new guess, xj + 1, from a previous one, xj. Sometimes Newton's method will fail to converge to a root. In this case, the program should terminate after many trials, perhaps 100.

Write a program that uses Newton's method to approximate the nth root of a number to six decimal places. If x^n = c, then x^n - c = 0. Finding a root of the second equation will give you n√c. Test your program on √2, 3√7, and 3√-1. Your program could use c/2 as its initial guess.
In this case, would f(x) = xj? :confused:
For $1000: Something that is a miserable pile of secrets.

#6
WingedPanther

WingedPanther

    A spammer's worst nightmare

  • Moderators
  • 16,831 posts
  • Location:Upstate, South Carolina
  • Programming Language:C, C++, PL/SQL, Delphi/Object Pascal, Pascal, Transact-SQL, Others
  • Learning:Java, C#, PHP, JavaScript, Lisp, Fortran, Haskell, Others
The idea is that xj is the j'th estimate of the true 0. f(x) is just the function (you could use a function pointer to make your code flexible). f(x) is whatever function you want, such as sqrt(x) or cubert(x).
Programming is a branch of mathematics.
My CodeCall Blog | My Personal Blog

#7
Yuriy M

Yuriy M

    Learning Programmer

  • Members
  • PipPipPip
  • 70 posts
Okay, here's the revised code:

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

double power_rule(double);                              /* Function that calculates the power rule of the equation */
double derivative(double);                              /* Function that calculates the derivative of the equation */

int main(void)
{
    /* Declare variables */
    
    FILE   *inp;                                        /* The input file */
    double curr_guess,                                  /* The current guess value for the root of the equation */
           new_guess,                                   /* The new guess value for the root of the equation */
           root_comp,                                   /* The comparison equation for a root match */
           c;                                           /* The current value to be used to determine the initial guess */
    int    trial_counter,                               /* Counts the guess trials for determining the root of the equation */
           inp_value,                                   /* The input value from the file */
           inp_counter = 0,                             /* Counts the input values from the file */
           input_status;                                /* The status of the file input */
    
    /* Open file */
    
    inp = fopen("root_test_values.dat", "r");
    input_status = fscanf(inp, "%d", &inp_value);
    
    /* Guess the root */
    
    while (input_status != EOF)
    {
     /* Check input value and assign appropriate square or cube root */
     if (inp_counter == 0)
     {
      c = pow(inp_value, 0.5);
      root_comp = pow(c, 0.5);
     }
     else if (inp_counter == 1)
     {
      c = pow(inp_value, 0.33);
      root_comp = pow(c, 0.33);
     }
     else
     {
      inp_value = inp_value * -1;
      c = pow(inp_value, 0.33);
      root_comp = -1 * pow(c, 0.33);
     }
     
     /* Provide starting guess */
     curr_guess = c / 2;
     printf("\nStarting guess for root is %.6f", curr_guess);
     
     if (curr_guess == root_comp)
     {
      printf("\nRoot found at %.6f\n", curr_guess);
      inp_counter++; 
      input_status = fscanf(inp, "%d", &inp_value);
     }
     else
     {
      /* Perform trials */
      for (trial_counter = 0; trial_counter < 100; trial_counter++)
      {
       new_guess = curr_guess - (power_rule(c) / derivative(c));
       printf("\nNew guess for root is %.6f", new_guess);
      
       if (new_guess == root_comp)
       {
        printf("\nRoot found at %.6f\n", new_guess);
        break;
       }
       else
        curr_guess = new_guess;
      } 
    
      if (inp_counter == 0 && trial_counter == 100)
       printf("\nNo root was found for the square root of %d\n", inp_value);
      if (inp_counter == 1 && trial_counter == 100)
       printf("\nNo root was found for the cube root of %d\n", inp_value);
      if (inp_counter == 2 && trial_counter == 100)
       printf("\nNo root was found for the cube root of %d\n", inp_value * -1);
      
      inp_counter++; 
      input_status = fscanf(inp, "%d", &inp_value);
     }
    }
     
    /* Close file */
    
    fclose(inp);
    
    return 0;   
}

double power_rule(double x)
{
    return (2 * pow(x, 2)) + (4 * pow(x, 3));
}

double derivative(double x)
{
    return (4 * x) + (12 * pow(x, 2));
}
The results appear to be looking good but I'm not sure If I'm supposed to be using "c" or "inp_value" in the function arguments for "power_rule" and "derivative". The code provided is using the argument "c" as part of the "new_guess" calculation.

P.S. For f(x), I'm using the power rule formula: 2x^2 + 4x^3 and for f`(x), the derivative: 4x + 12x^2. :)
For $1000: Something that is a miserable pile of secrets.

#8
WingedPanther

WingedPanther

    A spammer's worst nightmare

  • Moderators
  • 16,831 posts
  • Location:Upstate, South Carolina
  • Programming Language:C, C++, PL/SQL, Delphi/Object Pascal, Pascal, Transact-SQL, Others
  • Learning:Java, C#, PHP, JavaScript, Lisp, Fortran, Haskell, Others
Why are you testing inp_counter to decide what to do with inp_value?
Programming is a branch of mathematics.
My CodeCall Blog | My Personal Blog

#9
Yuriy M

Yuriy M

    Learning Programmer

  • Members
  • PipPipPip
  • 70 posts
Okay, how about this then?
#include <stdio.h>
#include <math.h>

double power_rule(double);                              /* Function that calculates the power rule of the equation */
double derivative(double);                              /* Function that calculates the derivative of the equation */

int main(void)
{
    /* Declare variables */
    
    FILE   *inp;                                        /* The input file */
    double curr_guess,                                  /* The current guess value for the root of the equation */
           new_guess,                                   /* The new guess value for the root of the equation */
           root_comp,                                   /* The comparison equation for a root match */
           c;                                           /* The current value to be used to determine the initial guess */
    int    trial_counter,                               /* Counts the guess trials for determining the root of the equation */
           inp_radicand,                                /* Input radicand of the considered root */
           calc_radicand,                               /* The radicand being calculated */
           degree;                                      /* The degree of the considered root */
    char   response;                                    /* The response to the availability of another radicand */
    
    /* Prompt user for availability of radicand */
    printf("\nProcess root? (Y/N): ");
    scanf("%c", &response);
    
    /* Guess the root */
    while (response == 'y' || response == 'Y')
    {
     /* Collect radicand and degree */
     printf("\nEnter radicand: ");
     scanf("%d", &inp_radicand);
     printf("\nEnter degree: ");
     scanf("%d", °ree);
     
     /* Error check */
     while (degree != 3 && inp_radicand < 0)
     {
      printf("\nNegative root values that are not cube roots cannot be calculated.\n"); 
      printf("\nEnter radicand: ");
      scanf("%d", &inp_radicand);
      printf("\nEnter root: ");
      scanf("%d", °ree); 
     }
     
     calc_radicand = inp_radicand;
     
     /* Check input value and assign appropriate root */
     if (degree == 3 && calc_radicand < 0)
     {
      calc_radicand = calc_radicand * -1;
      c = pow(calc_radicand, (float)1 / degree);
      root_comp = -1 * pow(c, (float)1 / degree);
     }
     else
     {
      c = pow(calc_radicand, (float)1 / degree);
      root_comp = pow(c, (float)1 / degree);
     }
      
     /* Provide starting guess */
     curr_guess = c / 2;
     printf("\nStarting guess for root is %.6f", curr_guess);
     
     if (curr_guess == root_comp)
     {
      printf("\nRoot found at %.6f\n", curr_guess);
      
      /* Prompt user for availability of radicand */
      printf("\nProcess root? (Y/N): ");
      scanf("%c", &response);
     }
     else
     {
      /* Perform trials */
      for (trial_counter = 0; trial_counter < 100; trial_counter++)
      {
       new_guess = curr_guess - (power_rule(c) / derivative(c));
       printf("\nNew guess for root is %.6f", new_guess);
      
       if (new_guess == root_comp)
       {
        printf("\nRoot found at %.6f\n", new_guess);
        break;
       }
       else
        curr_guess = new_guess;
      } 
    
      if (degree == 2 && trial_counter == 100)
       printf("\nNo root was found for the square root of %d\n", inp_radicand);
      else if (degree == 3 && trial_counter == 100 && calc_radicand < 0)
       printf("\nNo root was found for the cube root of %d\n", inp_radicand);
      else if (degree == 3 && trial_counter == 100 && calc_radicand >= 0)
       printf("\nNo root was found for the cube root of %d\n", inp_radicand);
      else if (degree > 3 && trial_counter == 100)
       printf("\nNo root was found for the %dth root of %d\n", degree, inp_radicand);
     }
     /* Prompt user for availability of radicand */
     printf("\nProcess root? (Y/N): ");
     scanf(" %c", &response);
    }
    
    return 0;   
}

double power_rule(double x)
{
    return (2 * pow(x, 2)) + (4 * pow(x, 3));
}

double derivative(double x)
{
    return (4 * x) + (12 * pow(x, 2));
}
Is this more like it?
For $1000: Something that is a miserable pile of secrets.

#10
WingedPanther

WingedPanther

    A spammer's worst nightmare

  • Moderators
  • 16,831 posts
  • Location:Upstate, South Carolina
  • Programming Language:C, C++, PL/SQL, Delphi/Object Pascal, Pascal, Transact-SQL, Others
  • Learning:Java, C#, PHP, JavaScript, Lisp, Fortran, Haskell, Others
Assuming it works, my only complaint with it right now is general structure. I found my eyes cross as I scanned through the code. However, it's pretty short, and I'm tired, so I wouldn't stress too much :)
Programming is a branch of mathematics.
My CodeCall Blog | My Personal Blog

#11
Yuriy M

Yuriy M

    Learning Programmer

  • Members
  • PipPipPip
  • 70 posts
Well, I'll give you a breakdown of the changes I made to the code:

  • Removed file input altogether. Radicands and degrees are now manually inputted.
  • inp_counter and inp_value removed as they defeat the program's purpose of universally finding the root of all potential values and not just the examples that are provided in the textbook.
  • The program now searches for root values with degrees greater than 3.
  • The program checks to see that negative radicand values can be calculated with a degree of 3 and not with any others.
I love status updates, don't you? ^^
For $1000: Something that is a miserable pile of secrets.

#12
WingedPanther

WingedPanther

    A spammer's worst nightmare

  • Moderators
  • 16,831 posts
  • Location:Upstate, South Carolina
  • Programming Language:C, C++, PL/SQL, Delphi/Object Pascal, Pascal, Transact-SQL, Others
  • Learning:Java, C#, PHP, JavaScript, Lisp, Fortran, Haskell, Others
negative radicand values work with odd degrees, not just 3 :)
Programming is a branch of mathematics.
My CodeCall Blog | My Personal Blog




1 user(s) are reading this topic

0 members, 1 guests, 0 anonymous users