c++ Newton-Raphson problem

Discussion in 'C++' started by pauldepstein@att.net, Nov 15, 2008.

  1. Guest

    Let double NR( double x, double(*)(const double&) f ) be the
    signature of a Newton-Raphson function NR.

    Here, f is a function which returns a double and accepts a const
    double&. The aim of the game is to find a zero of this function f
    (the point at which f crosses the x-axis). This zero-of-f which
    solves our problem is the double which NR returns. It remains to
    explain what the "double x" represents. This is the starting-guess
    that is required in Newton-Raphson implementations.

    In my case, I have the following amended Newton-Raphson situation. I
    have a function of the form

    double MyFunc(double x1, double x2, double x3, double x4, double x5)

    I want to solve the following problem: Fix x1, x2, x3, and x4. Then
    use Newton Raphson to return the double y such that MyFunc(x1, x2, x3,
    x4, y) = 0.

    I was unable to find a way of using the ready-made function NR because
    it assumes f accepts 1 double and returns 1 double, whereas My Func
    accepts 5 doubles and returns 1 double.

    My very-inelegant solution was to copy-paste the NR code and adapt it
    so that the pointer-to-function parameter was of the type I needed.

    Is there a more elegant approach that calls on the NR function already
    present?

    Thank for your help.

    Paul Epstein
     
    , Nov 15, 2008
    #1
    1. Advertising

  2. Kai-Uwe Bux Guest

    wrote:

    > Let double NR( double x, double(*)(const double&) f ) be the
    > signature of a Newton-Raphson function NR.
    >
    > Here, f is a function which returns a double and accepts a const
    > double&. The aim of the game is to find a zero of this function f
    > (the point at which f crosses the x-axis). This zero-of-f which
    > solves our problem is the double which NR returns. It remains to
    > explain what the "double x" represents. This is the starting-guess
    > that is required in Newton-Raphson implementations.
    >
    > In my case, I have the following amended Newton-Raphson situation. I
    > have a function of the form
    >
    > double MyFunc(double x1, double x2, double x3, double x4, double x5)
    >
    > I want to solve the following problem: Fix x1, x2, x3, and x4. Then
    > use Newton Raphson to return the double y such that MyFunc(x1, x2, x3,
    > x4, y) = 0.
    >
    > I was unable to find a way of using the ready-made function NR because
    > it assumes f accepts 1 double and returns 1 double, whereas My Func
    > accepts 5 doubles and returns 1 double.
    >
    > My very-inelegant solution was to copy-paste the NR code and adapt it
    > so that the pointer-to-function parameter was of the type I needed.
    >
    > Is there a more elegant approach that calls on the NR function already
    > present?


    I would change NR into a template:

    template < typename Float, typename Func >
    Float find_zero ( Float initial_guess, Func f );

    Then, you could use bind() from c++0x or Boost to fix the first four
    arguments and pass the resulting function object into the template.


    Best

    Kai-Uwe Bux
     
    Kai-Uwe Bux, Nov 15, 2008
    #2
    1. Advertising

  3. James Kanze Guest

    On Nov 15, 9:53 am, wrote:
    > Let double NR( double x, double(*)(const double&) f )  be the
    > signature of a Newton-Raphson function NR.


    > Here, f is a function which returns a double and accepts a
    > const double&. The aim of the game is to find a zero of
    > this function f (the point at which f crosses the x-axis).
    > This zero-of-f which solves our problem is the double which NR
    > returns. It remains to explain what the  "double x"
    > represents. This is the starting-guess that is required in
    > Newton-Raphson implementations.


    > In my case, I have the following amended Newton-Raphson
    > situation.  I have a function of the form


    > double MyFunc(double x1, double x2, double x3, double x4, double x5)


    > I want to solve the following problem:  Fix x1, x2, x3, and
    > x4.  Then use Newton Raphson to return the double y such that
    > MyFunc(x1, x2, x3, x4, y) = 0.


    > I was unable to find a way of using the ready-made function NR
    > because it assumes f accepts 1 double and returns 1 double,
    > whereas My Func accepts 5 doubles and returns 1 double.


    That's because the interface to the existing NR function is very
    poorly designed. In C++, the "standard" solution for any
    callback would be:

    class NRCallBack
    {
    public:
    virtual ~NRCallBack() {}
    virtual double operator()( double ) const = 0 ;
    } ;

    So the signature of NR would be:

    double NR( double x, NRCallBack const& f ) ;

    Rather than providing a function, you then derive from
    NRCallBack, and define the appropriate operator.

    In your precise case, it's probably a bit wordy, because we
    don't have lambda classes, and you'd have to do something like:

    double
    NRforMyFunc( double x1, double x2, double x3, double x4 )
    {
    class F : public NRCallBack
    {
    public:
    NRCallBack( double x1, double x2, double x3, double x4 )
    : x1( x1 )
    , x2( x2 )
    , x3( x3 )
    , x4( x4 )
    {
    }
    virtual double operator()( double x ) const
    {
    return MyFunc( x1, x2, x3, x4, x ) ;
    }

    private:
    double x1 ;
    double x2 ;
    double x3 ;
    double x4 ;
    } ;
    return NR( 0.0, F() ) ;
    }

    If (as may be the case), NR is in fact a C function, and must be
    callable from C, the established convention is to pass an
    additional void* with user data, i.e.:
    double NR( double x, double (*f)( double, void* ), void* ) ;
    Again, you have to write a wrapper function which takes the
    additional, fixed values as a void*, move these values into an
    array, and pass the address of the array to NR.

    > My very-inelegant solution was to copy-paste the NR code and
    > adapt it so that the pointer-to-function parameter was of the
    > type I needed.


    You may end up having to do this, if it's interface is broken.

    > Is there a more elegant approach that calls on the NR function
    > already present?


    Depending on the context of what you're doing, you may be able
    to use static variables and a wrapper function. IMHO, it's
    playing with fire, however, and you'd be better off rewriting
    the function to use one of the above interfaces, depending on
    whether it is pure C++, or it must be callable from C as well.

    --
    James Kanze (GABI Software) email:
    Conseils en informatique orientée objet/
    Beratung in objektorientierter Datenverarbeitung
    9 place Sémard, 78210 St.-Cyr-l'École, France, +33 (0)1 30 23 00 34
     
    James Kanze, Nov 15, 2008
    #3
  4. James Kanze Guest

    On Nov 15, 10:26 am, Kai-Uwe Bux <> wrote:
    > wrote:
    > > Let double NR( double x, double(*)(const double&) f )  be the
    > > signature of a Newton-Raphson function NR.


    > > Here, f is a function which returns a double and accepts a
    > > const double&.      The aim of the game is to find a zero of
    > > this function f (the point at which f crosses the x-axis).
    > >  This zero-of-f which solves our problem is the double which
    > > NR returns.   It remains to explain what the  "double x"
    > > represents.   This is the starting-guess that is required in
    > > Newton-Raphson implementations.


    > > In my case, I have the following amended Newton-Raphson situation.  I
    > > have a function of the form


    > > double MyFunc(double x1, double x2, double x3, double x4, double x5)


    > > I want to solve the following problem:  Fix x1, x2, x3, and x4.  Then
    > > use Newton Raphson to return the double y such that MyFunc(x1, x2, x3,
    > > x4, y) = 0.


    > > I was unable to find a way of using the ready-made function NR because
    > > it assumes f accepts 1 double and returns 1 double, whereas My Func
    > > accepts 5 doubles and returns 1 double.


    > > My very-inelegant solution was to copy-paste the NR code and adapt it
    > > so that the pointer-to-function parameter was of the type I needed.


    > > Is there a more elegant approach that calls on the NR function already
    > > present?


    > I would change NR into a template:


    >   template < typename Float, typename Func >
    >   Float find_zero ( Float initial_guess, Func f );


    > Then, you could use bind() from c++0x or Boost to fix the
    > first four arguments and pass the resulting function object
    > into the template.


    This is a very elegant solution for a few special cases, but it
    results in an infection template; if the call to this function
    is in a function which receives the callback function as an
    argument, that function must be a template as well. And so on,
    add infitum; depending on the use pattern, you can very quickly
    end up with an unmanageable mess, where all of your functions
    are templates. (This might be workable if your compiler
    supports export, but not many do.)

    --
    James Kanze (GABI Software) email:
    Conseils en informatique orientée objet/
    Beratung in objektorientierter Datenverarbeitung
    9 place Sémard, 78210 St.-Cyr-l'École, France, +33 (0)1 30 23 00 34
     
    James Kanze, Nov 15, 2008
    #4
  5. Kai-Uwe Bux Guest

    James Kanze wrote:

    > On Nov 15, 10:26 am, Kai-Uwe Bux <> wrote:
    >> wrote:
    >> > Let double NR( double x, double(*)(const double&) f )  be the
    >> > signature of a Newton-Raphson function NR.

    >
    >> > Here, f is a function which returns a double and accepts a
    >> > const double&.      The aim of the game is to find a zero of
    >> > this function f (the point at which f crosses the x-axis).
    >> > This zero-of-f which solves our problem is the double which
    >> > NR returns.   It remains to explain what the  "double x"
    >> > represents.   This is the starting-guess that is required in
    >> > Newton-Raphson implementations.

    >
    >> > In my case, I have the following amended Newton-Raphson situation.  I
    >> > have a function of the form

    >
    >> > double MyFunc(double x1, double x2, double x3, double x4, double x5)

    >
    >> > I want to solve the following problem:  Fix x1, x2, x3, and x4.  Then
    >> > use Newton Raphson to return the double y such that MyFunc(x1, x2, x3,
    >> > x4, y) = 0.

    >
    >> > I was unable to find a way of using the ready-made function NR because
    >> > it assumes f accepts 1 double and returns 1 double, whereas My Func
    >> > accepts 5 doubles and returns 1 double.

    >
    >> > My very-inelegant solution was to copy-paste the NR code and adapt it
    >> > so that the pointer-to-function parameter was of the type I needed.

    >
    >> > Is there a more elegant approach that calls on the NR function already
    >> > present?

    >
    >> I would change NR into a template:

    >
    >> template < typename Float, typename Func >
    >> Float find_zero ( Float initial_guess, Func f );

    >
    >> Then, you could use bind() from c++0x or Boost to fix the
    >> first four arguments and pass the resulting function object
    >> into the template.

    >
    > This is a very elegant solution for a few special cases, but it
    > results in an infection template; if the call to this function
    > is in a function which receives the callback function as an
    > argument, that function must be a template as well.


    Huh? I admit that this sentence has too many "this" and "that" for me to get
    references straight. So, I do not really understand what you mean. Anyway,
    I also do not see any reason why this template could not be called from
    ordinary functions:

    template < typename Float, typename Func >
    Float find_zero ( Float initial_guess, Func f ) {
    return ( 1 );
    }

    double caller ( double x, double(*f)(double) ) {
    return ( find_zero( x, f ) );
    }

    double id ( double x ) {
    return (x);
    }

    #include <iostream>
    #include <ostream>

    int main ( void ) {
    std::cout << caller( 1.0, &id ) << '\n';
    }


    (BTW: a Google search for "infection template" only yield medical stuff.)


    > And so on,
    > add infitum; depending on the use pattern, you can very quickly
    > end up with an unmanageable mess, where all of your functions
    > are templates. (This might be workable if your compiler
    > supports export, but not many do.)


    I don't see that. Could you please elaborate?


    Best

    Kai-Uwe Bux
     
    Kai-Uwe Bux, Nov 15, 2008
    #5
  6. James Kanze Guest

    On Nov 15, 11:48 am, Kai-Uwe Bux <> wrote:
    > James Kanze wrote:
    > > On Nov 15, 10:26 am, Kai-Uwe Bux <> wrote:
    > >> wrote:
    > >> > Let double NR( double x, double(*)(const double&) f )  be the
    > >> > signature of a Newton-Raphson function NR.


    > >> > Here, f is a function which returns a double and accepts
    > >> > a const double&.      The aim of the game is to find a
    > >> > zero of this function f (the point at which f crosses the
    > >> > x-axis). This zero-of-f which solves our problem is the
    > >> > double which NR returns.   It remains to explain what the
    > >> >  "double x" represents.   This is the starting-guess that
    > >> > is required in Newton-Raphson implementations.


    > >> > In my case, I have the following amended Newton-Raphson
    > >> > situation.  I have a function of the form


    > >> > double MyFunc(double x1, double x2, double x3, double x4, double x5)


    > >> > I want to solve the following problem:  Fix x1, x2, x3,
    > >> > and x4.  Then use Newton Raphson to return the double y
    > >> > such that MyFunc(x1, x2, x3, x4, y) = 0.


    > >> > I was unable to find a way of using the ready-made
    > >> > function NR because it assumes f accepts 1 double and
    > >> > returns 1 double, whereas My Func accepts 5 doubles and
    > >> > returns 1 double.


    > >> > My very-inelegant solution was to copy-paste the NR code
    > >> > and adapt it so that the pointer-to-function parameter
    > >> > was of the type I needed.


    > >> > Is there a more elegant approach that calls on the NR
    > >> > function already present?


    > >> I would change NR into a template:


    > >> template < typename Float, typename Func >
    > >> Float find_zero ( Float initial_guess, Func f );


    > >> Then, you could use bind() from c++0x or Boost to fix the
    > >> first four arguments and pass the resulting function object
    > >> into the template.


    > > This is a very elegant solution for a few special cases, but
    > > it results in an infection template; if the call to this
    > > function is in a function which receives the callback
    > > function as an argument, that function must be a template as
    > > well.


    > Huh? I admit that this sentence has too many "this" and "that"
    > for me to get references straight. So, I do not really
    > understand what you mean. Anyway, I also do not see any reason
    > why this template could not be called from ordinary functions:


    The problem is that the fact that it is a template, and is not
    resolved dynamically, propagates. If the ordinary function
    calls it with a known callback, fine; the propagation stops
    there. But if the ordinary function calls it with a functor
    that is passed in, then if templates are used, the ordinary
    function has to be a template too. And so on; imagine what
    would happen if iostream used templates, rather than virtual
    functions, in streambuf.

    > template < typename Float, typename Func >
    > Float find_zero ( Float initial_guess, Func f ) {
    >   return ( 1 );
    > }


    > double caller ( double x, double(*f)(double) ) {
    >   return ( find_zero( x, f ) );
    > }


    In this case, I fail to see what you've gained with respect to
    the initial problem.

    > double id ( double x ) {
    >   return (x);
    > }


    > #include <iostream>
    > #include <ostream>


    > int main ( void ) {
    >   std::cout << caller( 1.0, &id ) << '\n';
    > }


    > (BTW: a Google search for "infection template" only yield
    > medical stuff.)


    Yes. I'm not really sure what the appropriate word should be.
    The closest analogy I can think of is the GNU license, but of
    course, templates aren't that infectious. The fact remains that
    any time you implement genericity with a template, it means that
    any client code which wants to maintain that genericity neads to
    be a template as well. Anytime you have a real choice, you
    should prefer inheritance with virtual functions to templates.
    (Of course, most of the time you don't have a real choice;
    templates are designed to solve a different set of problems than
    virtual functions.)

    --
    James Kanze (GABI Software) email:
    Conseils en informatique orientée objet/
    Beratung in objektorientierter Datenverarbeitung
    9 place Sémard, 78210 St.-Cyr-l'École, France, +33 (0)1 30 23 00 34
     
    James Kanze, Nov 15, 2008
    #6
  7. Kai-Uwe Bux Guest

    James Kanze wrote:

    > On Nov 15, 11:48 am, Kai-Uwe Bux <> wrote:
    >> James Kanze wrote:
    >> > On Nov 15, 10:26 am, Kai-Uwe Bux <> wrote:
    >> >> wrote:
    >> >> > Let double NR( double x, double(*)(const double&) f )  be the
    >> >> > signature of a Newton-Raphson function NR.

    >
    >> >> > Here, f is a function which returns a double and accepts
    >> >> > a const double&.      The aim of the game is to find a
    >> >> > zero of this function f (the point at which f crosses the
    >> >> > x-axis). This zero-of-f which solves our problem is the
    >> >> > double which NR returns.   It remains to explain what the
    >> >> > "double x" represents.   This is the starting-guess that
    >> >> > is required in Newton-Raphson implementations.

    >
    >> >> > In my case, I have the following amended Newton-Raphson
    >> >> > situation.  I have a function of the form

    >
    >> >> > double MyFunc(double x1, double x2, double x3, double x4, double x5)

    >
    >> >> > I want to solve the following problem:  Fix x1, x2, x3,
    >> >> > and x4.  Then use Newton Raphson to return the double y
    >> >> > such that MyFunc(x1, x2, x3, x4, y) = 0.

    >
    >> >> > I was unable to find a way of using the ready-made
    >> >> > function NR because it assumes f accepts 1 double and
    >> >> > returns 1 double, whereas My Func accepts 5 doubles and
    >> >> > returns 1 double.

    >
    >> >> > My very-inelegant solution was to copy-paste the NR code
    >> >> > and adapt it so that the pointer-to-function parameter
    >> >> > was of the type I needed.

    >
    >> >> > Is there a more elegant approach that calls on the NR
    >> >> > function already present?

    >
    >> >> I would change NR into a template:

    >
    >> >> template < typename Float, typename Func >
    >> >> Float find_zero ( Float initial_guess, Func f );

    >
    >> >> Then, you could use bind() from c++0x or Boost to fix the
    >> >> first four arguments and pass the resulting function object
    >> >> into the template.

    >
    >> > This is a very elegant solution for a few special cases, but
    >> > it results in an infection template; if the call to this
    >> > function is in a function which receives the callback
    >> > function as an argument, that function must be a template as
    >> > well.

    >
    >> Huh? I admit that this sentence has too many "this" and "that"
    >> for me to get references straight. So, I do not really
    >> understand what you mean. Anyway, I also do not see any reason
    >> why this template could not be called from ordinary functions:

    >
    > The problem is that the fact that it is a template, and is not
    > resolved dynamically, propagates. If the ordinary function
    > calls it with a known callback, fine; the propagation stops
    > there. But if the ordinary function calls it with a functor
    > that is passed in, then if templates are used, the ordinary
    > function has to be a template too. And so on; imagine what
    > would happen if iostream used templates, rather than virtual
    > functions, in streambuf.


    You keep using pronouns and terms like "ordinary function" whose reference
    is not clear to me. Maybe, I am just dense.

    So for concreteness, here is a simple template implementation of find_zero
    using the Newton method:

    template < typename Float, typename Func >
    Float find_zero ( Func f,
    Float initial_guess,
    Float eps = std::numeric_limits<Float>::epsilon() * 1000,
    Float h = std::numeric_limits<Float>::epsilon() * 1000,
    unsigned int iteration_limit = 50 ) {
    Float x = initial_guess;
    for ( unsigned int n = 0; n < iteration_limit; ++n ) {
    Float y = f(x);
    if ( std::abs(y) < eps ) {
    return (x);
    }
    x = x - y * h / ( f(x+h/2) - f(x-h/2) );
    }
    throw( std::invalid_argument( "iteration limit exceeded" ) );
    }

    I fail to see any problem with that. Clearly, there is no requirement that
    the call back function f has to be a template or that this function can
    only be called from templates. If one uses the signature

    template < typename Float, typename Func >
    Float find_zero ( Func const & f,
    Float initial_guess,
    Float eps = std::numeric_limits<Float>::epsilon() * 1000,
    Float h = std::numeric_limits<Float>::epsilon() * 1000,
    unsigned int iteration_limit = 50 );

    one can even support the inheritance based solution you suggested
    elsethread.


    If you could provide code that illustrates the propagation of templates with
    this example, it would be highly appreciated.


    >> template < typename Float, typename Func >
    >> Float find_zero ( Float initial_guess, Func f ) {
    >> return ( 1 );
    >> }

    >
    >> double caller ( double x, double(*f)(double) ) {
    >> return ( find_zero( x, f ) );
    >> }

    >
    > In this case, I fail to see what you've gained with respect to
    > the initial problem.


    Nothing: the code is just meant to illustrate why I don't see that something
    propagates. That why the function is called "caller". Clearly, I am still
    missing your point.


    >> double id ( double x ) {
    >> return (x);
    >> }

    >
    >> #include <iostream>
    >> #include <ostream>

    >
    >> int main ( void ) {
    >> std::cout << caller( 1.0, &id ) << '\n';
    >> }

    >
    >> (BTW: a Google search for "infection template" only yield
    >> medical stuff.)

    >
    > Yes. I'm not really sure what the appropriate word should be.
    > The closest analogy I can think of is the GNU license, but of
    > course, templates aren't that infectious. The fact remains that
    > any time you implement genericity with a template, it means that
    > any client code which wants to maintain that genericity neads to
    > be a template as well.


    Here the key phrase is "wants to maintain that genericity". There is no
    reason for client code to always be like that. All algorithms is the
    standard library are templates and they are used in non-templated code. I
    do not see the infectious trait you mention.


    > Anytime you have a real choice, you
    > should prefer inheritance with virtual functions to templates.


    I disagree. In my codebase, that would be utterly inappropriate. More to the
    point: finding roots of functions clearly seems like a task for a template
    and forcing clients to inherit from some call back just to fit the
    interface is clearly wrong.


    > (Of course, most of the time you don't have a real choice;
    > templates are designed to solve a different set of problems than
    > virtual functions.)


    Agreed.


    Best

    Kai-Uwe Bux
     
    Kai-Uwe Bux, Nov 16, 2008
    #7
    1. Advertising

Want to reply to this thread or ask your own question?

It takes just 2 minutes to sign up (and it's free!). Just click the sign up button to choose a username and then you can ask your own questions on the forum.
Similar Threads
  1. Andreas Suurkuusk
    Replies:
    0
    Views:
    3,984
    Andreas Suurkuusk
    Jul 27, 2003
  2. Ted Miller
    Replies:
    0
    Views:
    5,158
    Ted Miller
    Sep 13, 2003
  3. moi

    newton-raphson problem

    moi, Nov 5, 2003, in forum: C++
    Replies:
    2
    Views:
    870
    Rajeev Ayyagari
    Nov 5, 2003
  4. Mike

    Problem problem problem :( Need Help

    Mike, May 7, 2004, in forum: ASP General
    Replies:
    2
    Views:
    551
    Bullschmidt
    May 11, 2004
  5. Theodore Knab

    Ruby math: Newton's Law of Cooling

    Theodore Knab, Feb 1, 2004, in forum: Ruby
    Replies:
    5
    Views:
    553
    Josef 'Jupp' SCHUGT
    Feb 9, 2004
Loading...

Share This Page