Newton-Raphson method with approximate derivative

[Claudius7]

I have written the following code for the Newton Raphson method, which takes as input two function pointers, one that points to the function f(x), the other points to its derivative f'(x). Once I dereference f(x) to a specific function, and f'(x) to its derivative (example: f(x) = x3 - 10, f'(x) = 3x^2), I can call the Newton routine on that function. So far, everything works smoothly:

#include <iostream>
#include <cmath>

using namespace std;

double Newton(double (*pt_funct)(double), double (*pt_deriv)(double), double &x, const double& err)
{
    while (fabs((*pt_funct)(x)) > err)
    {
        x -= (*pt_funct)(x) / (*pt_deriv)(x);
    }
    return x;
}

double funct(double x)
{
    return x * x * x - 10.0;    //Take for example the function y=x^3-10
}

double deriv(double x)          //Its derivative is the function y=3x^2
{
    return 3 * x * x;
}


int main()
{
    double initial = 1.0;

    cout << "Newton method applied to y = x^3 - 10, initial point x_0 = 1.0, tolerance = 0.01: " << endl;
    
    cout << Newton(funct, deriv, initial, 0.1) << endl;

    return 0;
}

Now, I want to extend my code in the following way: suppose I have a very complex expression for f(x), so that f'(x) cannot be evaluated with a simple algebraic expression as in the case above, but can be only approximated, for instance by f(x+0.001) - f(x-0.001)/0.002.

My question is: how should I modify my code in order to allow for this more general situation? I would like to be able to write at some point something like:

double(*pt_deriv)(double(x))
{
return 
((*pt_funct)(x+0.001) - (*pt_funct)(x+0.001)) / 0.002;
}

but it fails horribly. What do you think can be done?
Thanks in advance for your ideas.

[jonnin]

are you asking how to replace the function for your equation with a function pointer? There are tons of explains online on how to do this, though its generally used more in C than C++.

Also, you can just use the derivative formula once generically .. ( (f(x)- (fx2))/stuff ) .. if you get tired of coding up both the function and its dx.

another approach would be to just parse an equation format from a string.

[lastchance]

Like so?

#include <iostream>
#include <cmath>
using namespace std;

double Newton( double f(double), double dfdx(double fn(double), double), double &x, const double &err)
{
    while ( fabs( f( x ) ) > err )
    {
        x -= f( x ) / dfdx( f, x );
    }
    return x;
}


double funct( double x )
{
    return x * x * x - 10.0;    //Take for example the function y=x^3-10
}


double deriv( double f(double), double x )
{
   return ( f( x + 0.001 ) - f( x - 0.001 ) ) / 0.002;
}


int main()
{
    double initial = 1.0;
    double tolerance = 0.001;

    cout << "Newton method, initial point = " << initial << "  , tolerance = " << tolerance << endl;
    
    cout << Newton( funct, deriv, initial, tolerance ) << endl;

    return 0;
}

I've turned down your tolerance (a lot). Also, you probably need to put in some check that it doesn't just go on looping endlessly: N-R isn't guaranteed to converge.

EDITED: edited to avoid multiple deriv() name clashes.

[Claudius7]

@lastchance
thank you. But are you using function pointers in your code?? Your 'deriv' is a function, a pointer to a function, a composite of two functions.... what is exactly?

[lastchance]

As far as its effect goes, it's the function being passed as an argument. I simply find it easier to write down like this. It does exactly the same as with the (*f) notation. Feel free to put the *'s back in your own code if you prefer.

And yes, the C++ book that I learnt from told me to write them your way in the function declaration ... but then admitted it wasn't necessary.

EDITED: As you can see from the deriv() declaration in my code on line 21, it takes two arguments: a function and a value of x. This needs to be reflected in the Newton declaration on line 5. (I've edited the name here as I realised I was over-using the deriv() name.

[Claudius7]

Hi thanks a lot.
However, admittedly my code has problems now: I have now the following

#include <iostream>
#include <cmath>

using namespace std;

double Newton(double (*pt_funct)(double), double (*pt_deriv)(double (*f)(double), double),
              double &x, const double& err)
{
    while (fabs((*pt_funct)(x)) > err)
    {
        x -= (*pt_funct)(x) / (*pt_deriv)((*pt_funct), x);
    }
    return x;
}

double funct(double x)
{
    return sqrt(x) - 2.0;    //Want: approximate the square root of 2 using Newton method.
}

double deriv(double (*f)(double), double x)          //Approximate derivative.
{
    return  ((*f)( x + 0.001 ) - (*f)( x - 0.001 ) ) / 0.002;
}


int main()
{
    double initial = 4.0;

    cout << "Newton method applied to y = sqrt(x) -2, initial point x_0 = 4.0, 
tolerance = 0.001: " << endl;

    cout << Newton(funct, deriv, initial, 0.001) << endl;

    return 0;
}

Code compiles nicely, but when I run it, instead of giving me an approximate value for sqrt(2), all I get is 4.0... there should be something wrong, but I don't see exactly what and why. Can you spot the problem??

[integralfx]

Code compiles nicely, but when I run it, instead of giving me an approximate value for sqrt(2), all I get is 4.0... there should be something wrong, but I don't see exactly what and why. Can you spot the problem??

double funct(double x)
{
    return sqrt(x) - 2.0;    //Want: approximate the square root of 2 using Newton method.
}

That's the incorrect function to be using. Likewise with your derivative.
x² = 2
f(x) = x² - 2
f'(x) = 2x

Also, similar post:
http://www.cplusplus.com/forum/beginner/209898/

[Claudius7]

Oooops! Sorry, how silly of me! Thanks, integralfix. Now it works perfectly indeed.