Newton's Method

Hey everyone. I am trying to write a program using Newtons Method. I know there is no code here, but that's my problem... I don't know where to start. I know the method is similar to bisection so if anyone could explain either to me or just help me get started, please let me know. Here is my problem:

I am trying to write a program in C++ to find a root of an equation. It starts with an initial guess, x(0), and generates successive approximate roots x(1), x(2), ... x(j), x(j+1) using the formula:
x(j+1) = x(j) - ((fx(j)) / f'x(j))
{x sub j+1 equals x sub j minus f(x sub j) divided by f prime (x sub j)} where f'(x sub j) is the derivative of function f evaluated at x = x(j). The formula generates a new guess, x(j=1), from a previous one, f(j). The program should terminate after 100 trials.
From geometry we get the equation:
(y(j+1) - y(j))/(x(j+1) - x(j)) = m
where m is the slope of the line between points (x(j+1), y(j+1)) and (x(j), y(j)). y(j+1) is zero, y(j) is f(x(j)), and m is f'(x(j)). After rearranging we get:
-f(x(j)) = f'(x(j)) * (x(j+1) - x(j))
Write a program that uses Newtons method to approximate the nth root of a number to 6 decimal places.
If x^n = c thenx^n - c = 0.
Finding the rood of the second equation you get nroot(c).
Test the program on sqroot(2), 3root(7), and 3root(-1).
The program can use c/2 as its initial guess.

Code:

---

#include<iostream>
#include<cmath>
using namespace std;
int main(){

      return 0;
}

---

That is how I started to set it up. I just don't know where to start... Thanks for any help. :)

you need to define some tolerance: 10^-6 as the stopping condition. I would do it in the preprocessor.
you need a c++ function that take in double value x, and return double value for the derivative f'(x).
you need a c++ function that take in double value of x, and return double value for the function f(x);
then, just for-loop the heck out of it;

Code:

---

for ( int i = 0 ; i < 100 ; ++ i)
{
      x -= f(x)/f'(x);
      if ( abs(x) <= tolerance) break;
}

---

lol alright thanks. :) Ill let ya know how it goes.. workin on it now.

How can you calculate a derivative in c++? I know how on paper but don't know how to type it up...

there is no built-in symbolic derivative for c++. you have the option of doing the derivative on paper then write it into your f'(x) function. like if f(x) = x^2; i would write

Code:

---

double f_prime(double x)
{
 return 2.0*x;
}

---

you can calculate the derivative numerically from the function f(x), but i don't think that's part of your assignment.

Yeah. I'm starting to feel overwhelmed... :( lol I usually get these program fast but I'm stuck on this one... I know f(x) is the function of the curve, I think, so how do I even get that equation? would it be the slope equation y - y1 = m (x - x1) and solve for y? I need the derivative of f(x) but what excactly is that?

I'm a math major so I can understand the concept of Newton's method but am so confused with the variables and equations..

Newton's method - Wikipedia, the free encyclopedia
look at the examples.

Quote:

---

would it be the slope equation y - y1 = m (x - x1) and solve for y

---

not quite, what you think of must be local linear approximation of sqrt(2), that could be done in one step without all the iterations.
for your test cases; your funtion f(x) should be

Code:

---

f(x) = x^n - c;

---

such that f( nroot(c) ) = 0 .

I'll give you a hint here

Code:

---

double f(double x, int n, int c)

---

How will I know what f(x) is if the user only enters the number to find the nth root of and n (the degree of the root)?

What do you mean?
In this type of program, the user must provide the functions for f(x) and f'(x).

Edit: i think i know what you meant. That's why I gave you that function prototype as a hint. that function take in the value of x, the value of the root n, and the value of c. there is only one version of f(x) that you need for this problem

Code:

---

f(x) = x^n - c;

---

Thats what I thought but the end of the instructions said test your program on sqroot(2) 3rd root(7) and 3rd root(-1)... I just asked my professor because that didn't make sense to me.

so for example, 3rd root(7)
n = 3
c = 7
your function in the Newton's method become:

Code:

---

x -= -f(x,3,7)/f_prime(x,3,7);

---

So, maybe you didn't get what Newton's method does. Using the example above, you want to estimate x=7^(1/3) (or 3rd root of 7).
I can cube both sides to get x^3 = 7. Then if I substract both side by 7, we'll get x^3-7=0. So you need to find an x that satisfis that condition x^3-7=0, which is equilvalent to find a root of a function. So for any x that satifises x^3-7=0, x=7^1/3. So instead of using local linear approx on 3rdroot of 7, you iterater on f(x)=x^3-7.

Alright I think I get it. I'll work on it some more... I'll let ya know how it goes tomorrow prob.

Thanks. :)

Okay I think I have everything done, but how do you set a tolerance? Like I have to set a tolerance of 10^-6 (6 decimal places). Would it be like:

Code:

---

#define TOLERANCE 10^-6

---

? That's what I tried but I couldn't get it to work. I'll post the code I have now...

Code:

---

#include<iostream>
#include<cmath>

void getInitGuess(double&);
double f(double, int, int);
double fPrime(double, int, int);
void findRoots(double&, int, int);
using namespace std;
int main(){
        double root, guess;
        int c, n;

        cout << "Enter the number to find the nth root of: ";
        cin >> c;
        cout << "Enter n (the degree of the root): ";
        cin >> n;

        findRoots(guess, n, c);

        root = guess;
        cout << "The " << n << " root of " << c << " is: " << root << endl;
 
        return 0;
}

void getInitGuess(double& guess){
        cout << "Enter your initial guess: ";
        cin >> guess;
}

double f(double x, int n, int c){
        double fx;

        fx = pow(x, n) - c;

        return fx;
}

double fPrime(double x, int n, int c){
        double fp;

        fp = n * pow(x, (n-1));

        return fp;
}

void findRoots(double& x, int n, int c){
        double TOLERANCE = .000001;                      //Here is where I defined tolerance, but is this right?
        x = 0;
        for(int i=0; i < TOLERANCE; i++){
                x -= (f(x, n, c) / fPrime(x, n, c));
                if(abs(x) <= TOLERANCE)
                        break;
        }
}

---

When I run it as it is, this happens:
-------
Enter the number to find the nth root of: 2
Enter n (the degree of the root): 2
The 2 root of 2 is 1.#INF

The answer is the same for 3rd root of 7 and 3rd root of -1

What is #INF? How could I fix this? I think this is all I'm lacking now.

Code:

---

for(int i=0; i < TOLERANCE; i++

---

should be

Code:

---

for(int i=0; i < 100; i++

---

Also;

Code:

---

void findRoots(double& x, int n, int c){
        double TOLERANCE = .000001;                      //Here is where I defined tolerance, but is this right?
        x = 0; // GET RID OF THIS LINE;

---

That didn't change anything that I can see. The answer still comes out to -1.#INF. Even when I try sqroot(16) I get the same answer (-1.#INF not 4). Why would I change TOLERANCE to 100? Since the program has to have a tolerance of .000001, would I have to use #include<iomanip> and setprecision(6)? I forgot how to define a tolerance..

Quote:

---

Originally Posted by hottiefee

That didn't change anything that I can see. The answer still comes out to -1.#INF. Even when I try sqroot(16) I get the same answer (-1.#INF not 4). Why would I change TOLERANCE to 100? Since the program has to have a tolerance of .000001, would I have to use #include<iomanip> and setprecision(6)? I forgot how to define a tolerance..

---

I edited the other post, you need to get rid of the line "x=0" in findRoots function.

Also you need

Code:

---

guess = c/2;

---

in main() before you call findRoot();

Awesome. It's working now :) Thanks! The only thing is how do I set it to where it always shows the answer with 6 decimal places? I tried setpresicion(7) but that displays 7 digits instead of 6 decimal places.

Use fixed along with setpricision.

Code:

---

cout<< fixed << setprecison(6).

---

Oh yeah!! Thanks. I haven't used that since the begining of last semester so I forgot it. Thank you so much for the help :)

No problem.

Quote:

---

Originally Posted by nimitzhunter

Use fixed along with setpricision.

Code:

---

cout<< fixed << setprecison(6).

---

---

I believe you will also need the library,

Code:

---

#include <iomanip>

---

to use setprecision.

#define TOLERANCE 10^-6
this is incorrect. "^" is not the power operator.
I believe there is a syntax such as
const float g_Tolerance = 10E-6;
Can't remember exactly.
Also, use const instead of define in C++.

This is what i came up with, not the most efficient but it works

Code:

---

double Sqrt(double num)
{
    return Guess(num/2, num, 0, 256);
}
double Guess(double x, double y, int depth, int maxDepth)
{
    if (++depth < maxDepth)
        return Guess(x - (((x * x) - y) / (2 * x)), y, depth, maxDepth);
    return x;
}

---

What? You invented your own sqrt or something?

Quote:

---

Originally Posted by Elysia

What? You invented your own sqrt or something?

---

Using Newton's Method, yes i did.

Quote:

---

Originally Posted by Elysia

#define TOLERANCE 10^-6
this is incorrect. "^" is not the power operator.
I believe there is a syntax such as
const float g_Tolerance = 10E-6;
Can't remember exactly.
Also, use const instead of define in C++.

---

I just can't shake MATLAB out of my system sometimes.

For a bit of fun, I thought I'd write my version of this program. Criticisms encouraged :-).

Code:

---

#include <cassert>
#include <cmath>
#include <cstdlib>
#include <iostream>

#include "boost/bind.hpp"
#include "boost/cstdint.hpp"

template <typename Function1, typename Function2>
double NewtonsMethod(
    Function1 f,
    Function2 fPrime,
    double x,
    double tolerance,
    boost::uint32_t maxIterations)
{
    assert(tolerance > 0);
    assert(maxIterations != 0);

    double last;
    boost::uint32_t iterations = 0;

    do
    {
        last = x;
        x -= f(x) / fPrime(x);
        ++iterations;
    }
    while (std::abs(x - last) > tolerance && iterations != maxIterations);

    return x;
}

double XPowNMinusC(double x, boost::uint32_t n, double c) { return std::pow(x, n) - c; }
double XPowNMinusCPrime(double x, boost::uint32_t n, double c) { return n * std::pow(x, n - 1); }

int main(int argc, char** argv)
{
    if (argc != 3)
    {
        std::cerr << "Usage: " << argv[0] << " n x" << std::endl;
        std::cerr << "(Prints the nth root of x)" << std::endl;
        return 1;
    }

    boost::uint32_t n;
    double x;

    n = atoi(argv[1]);
    x = atof(argv[2]);

    double nthRootOfX = NewtonsMethod(
        boost::bind(XPowNMinusC, _1, n, x), // f
        boost::bind(XPowNMinusCPrime, _1, n, x), // f'
        x / 2, // x[0]
        10E-6, // tolerance
        100); // max iterations

    std::cout << nthRootOfX << std::endl;

    return 0;
}

---

The only thing I can say is avoid atoi (use boost::lexical_cast instead) and I would prefer lambdas over boost::bind (error messages are horrible as usual).
Otherwise it looks like a very nice C++-like solution.

Quote:

---

Originally Posted by Elysia

The only thing I can say is avoid atoi (use boost::lexical_cast instead) and I would prefer lambdas over boost::bind (error messages are horrible as usual).
Otherwise it looks like a very nice C++-like solution.

---

Ah yes boost::lexical_cast would be better. I'm not very familiar with C++0x's lambda functions, can they be used in that program without changing NewtonsMethod() ?

Yes. Since lambdas behave just like functors, and you are using template parameters, it will be fine.

Updated with lambda functions and boost::lexical_cast. I also decided to go with fp for the derivative name instead of fPrime, and that x is a less arbitrary initial guess than x / 2.

Code:

---

#include <cassert>
#include <cmath>
#include <cstdlib>
#include <iostream>

#include "boost/cstdint.hpp"
#include "boost/lexical_cast.hpp"

template <typename Function1, typename Function2>
double NewtonsMethod(
    Function1 f,
    Function2 fp,
    double x,
    double tolerance,
    boost::uint32_t maxIterations)
{
    assert(tolerance > 0);
    assert(maxIterations != 0);
   
    double last;
    boost::uint32_t iterations = 0;

    do 
    { 
        last = x;
        x -= f(x) / fp(x);
        ++iterations;
    } 
    while (std::abs(x - last) > tolerance && iterations != maxIterations);

    return x;
}

int main(int argc, char** argv)
{
    if (argc != 3)
    { 
        std::cerr << "Usage: " << argv[0] << " n x" << std::endl;
        std::cerr << "(Prints the nth root of x)" << std::endl;
        return 1;
    } 

    boost::uint32_t n;
    double x;

    n = boost::lexical_cast<boost::uint32_t>(argv[1]);
    x = boost::lexical_cast<double>(argv[2]);

    double nthRootOfX = NewtonsMethod(
        [n, x] (double v) -> double { return std::pow(v, n) - x; }, // f
        [n, x] (double v) -> double { return n * std::pow(v, n - 1); }, // f'
        x, // x[0]
        10E-6, // tolerance
        100); // max iterations

    std::cout << nthRootOfX << std::endl;

    return 0;
}

---

Looks nice. The only complaint is that lexical_cast throws if it fails. You should probably intercept this and display an error.
Are you sure the derivative is right, though? Looks like you've missed a "-1".

Ha, every solution brings another problem :-P. I don't know quite what it is about exception handling, but I'm resistant to doing it. I guess I'm more concerned with how my code looks than how the program it produces looks! Maybe I should think twice about those feelings. Here's the new code:

Code:

---

#include <cassert>
#include <cmath>
#include <cstdlib>
#include <iostream>

#include "boost/cstdint.hpp"
#include "boost/lexical_cast.hpp"

template <typename Function1, typename Function2>
double NewtonsMethod(
    Function1 f,
    Function2 fp,
    double x,
    double tolerance,
    boost::uint32_t maxIterations)
{
    assert(tolerance > 0);
    assert(maxIterations != 0);
   
    double last;
    boost::uint32_t iterations = 0;
   
    do
    {
        last = x;
        x -= f(x) / fp(x);
        ++iterations;
    } 
    while (std::abs(x - last) > tolerance && iterations != maxIterations);
   
    return x;
}

void UsageErrorMessage(const char* programName)
{
    std::cerr << "Usage: " << programName << " n x" << std::endl;
    std::cerr << "(Prints the nth root of x)" << std::endl;
}

int main(int argc, char** argv)
{
    if (argc != 3)
    { 
        UsageErrorMessage(argv[0]);
        return 1;
    } 
   
    boost::uint32_t n;
    double x;
   
    try
    {
        n = boost::lexical_cast<boost::uint32_t>(argv[1]);
        x = boost::lexical_cast<double>(argv[2]);
    }
    catch (boost::bad_lexical_cast)
    {
        std::cerr << "Error: Could not convert arguments to numbers" << std::endl;
        UsageErrorMessage(argv[0]);
        return 1;
    }
   
    double nthRootOfX = NewtonsMethod(
        [n, x] (double v) -> double { return std::pow(v, n) - x; }, // f
        [n, x] (double v) -> double { return n * std::pow(v, n - 1); }, // f'
        x, // x[0]
        10E-6, // tolerance
        100); // max iterations
   
    std::cout << nthRootOfX << std::endl;
   
    return 0;
}

---