How to Solve a Root...

Discussion in 'C++' started by Mark, Dec 2, 2005.

  1. Mark

    Mark Guest

    always wondered how a computer might go about doing this...
    i'm taking calculus at university... and we went over newton's method.
    i found this interesting...

    ------------------------------
    #include <cstdlib>
    #include <iostream>
    #include <cmath>

    using namespace std;

    double f1(double x);
    double f2(double x, double y);

    int main(int argc, char *argv[])
    {
    double x = 5;
    cout << f1(x) << endl;
    cout << sqrt(x) << endl;

    system("PAUSE");
    return EXIT_SUCCESS;
    }

    double f1(double x)
    {
    return f2(x,x/2);
    }

    double f2(double x, double y)
    {
    double z = y - (y*y-x)/(2*x);
    if( y == z ) return z;
    else return f2(x,z);
    }
    ---------------------------

    will work for cube roots and others if you just tweak the formula.
     
    Mark, Dec 2, 2005
    #1
    1. Advertising

  2. Mark

    Mark Guest

    (just thought i'd share. no question involved)
     
    Mark, Dec 2, 2005
    #2
    1. Advertising

  3. Mark

    Guest

    Mark wrote:
    > always wondered how a computer might go about doing this...
    > i'm taking calculus at university... and we went over newton's method.
    > i found this interesting...
    >
    > ------------------------------
    > #include <cstdlib>
    > #include <iostream>
    > #include <cmath>
    >
    > using namespace std;
    >
    > double f1(double x);
    > double f2(double x, double y);
    >
    > int main(int argc, char *argv[])
    > {
    > double x = 5;
    > cout << f1(x) << endl;
    > cout << sqrt(x) << endl;
    >
    > system("PAUSE");
    > return EXIT_SUCCESS;
    > }
    >
    > double f1(double x)
    > {
    > return f2(x,x/2);
    > }
    >
    > double f2(double x, double y)
    > {
    > double z = y - (y*y-x)/(2*x);
    > if( y == z ) return z;


    I thought it was bad to directly compare floats?

    > else return f2(x,z);
    > }
    > ---------------------------
    >
    > will work for cube roots and others if you just tweak the formula.


    I love recursion as much as the next guy, but in this mathematical case
    I see huge potential for a stack overflow. Maybe this should be rolled
    into a loop instead. Might actually look cleaner.

    Cheers,
    Andre
     
    , Dec 2, 2005
    #3
  4. Mark

    Guest

    Mark wrote:
    > always wondered how a computer might go about doing this...
    > i'm taking calculus at university... and we went over newton's method.


    > will work for cube roots and others if you just tweak the formula.


    You might want to look into the Runge-Kutta method. Quite robust. I
    coded it up once many years ago...
     
    , Dec 2, 2005
    #4
  5. Mark

    Mark Guest

    wrote:

    > I thought it was bad to directly compare floats?


    i have no idea... never heard that before.

    > I love recursion as much as the next guy, but in this mathematical case
    > I see huge potential for a stack overflow. Maybe this should be rolled
    > into a loop instead. Might actually look cleaner.


    well. when i've worked these out on paper... they only took about 6 or
    8 "loops" to get it accurate to about 8 decimal places.... it really
    shouldn't stack that much at all. *i think*

    but yes, maybe a for loop would be more efficient...but... oh well.


    wrote:
    > You might want to look into the Runge-Kutta method. Quite robust. I
    > coded it up once many years ago...


    Runge-kutta eh? perhaps i shall look into it. you wouldnt happen to
    know what method cmath uses, would you?

    i suppose i could open the library and check myself.... if i can
    understand it.
     
    Mark, Dec 2, 2005
    #5
  6. Mark

    Mark P Guest

    Mark wrote:
    > wrote:
    >
    >
    >>I thought it was bad to directly compare floats?

    >
    >
    > i have no idea... never heard that before.
    >


    "Bad" may be too strong, but it should be done with care. A computer
    will compare 2.0 and 1.999999 as unequal; sometimes this is undesirable.
    Often rather than comparing for equality it makes more sense to see if
    the absolute value of the difference is less than some small number (say
    0.0001).

    >
    > wrote:
    >
    >>You might want to look into the Runge-Kutta method. Quite robust. I
    >>coded it up once many years ago...

    >


    Runge-Kutta is a numerical method for solving differential equations.
    Quite different from root finding.

    -Mark
     
    Mark P, Dec 2, 2005
    #6
  7. Mark wrote:

    > well. when i've worked these out on paper... they only took about 6 or
    > 8 "loops" to get it accurate to about 8 decimal places.... it really
    > shouldn't stack that much at all. *i think*


    Seems not fast enough. Some years ago, talking about fixed point
    calculations for writing graphics demos, I recommended to try the Hero
    method (a Newton's variant specific for square roots), and the result was
    that just 4 iterations were enough for 3-d calculus.

    Buy will be better to talk about this things in group about graphics
    programming.

    --
    Salu2
     
    =?ISO-8859-15?Q?Juli=E1n?= Albo, Dec 2, 2005
    #7
  8. Mark

    Guest

    Mark P wrote:

    > Runge-Kutta is a numerical method for solving differential equations.
    > Quite different from root finding.


    Not really. Often these forms can be converted into each other,
    certainly
    for the trivial forms used here. And sqr(x)=y implies y*y-x = 0, which
    is
    truly trivial. Runga-Kutta takes three points, IIRC, which means it's a
    single iteration for second-degree functions like that. Newton uses two
    points, which means it's only a single iteration for first-degree
    functions
    (lines).

    HTH,
    Michiel.
     
    , Dec 2, 2005
    #8
  9. Mark

    Guest

    Mark wrote:
    > wrote:
    > > I love recursion as much as the next guy, but in this mathematical case
    > > I see huge potential for a stack overflow. Maybe this should be rolled
    > > into a loop instead. Might actually look cleaner.

    >
    > well. when i've worked these out on paper... they only took about 6 or
    > 8 "loops" to get it accurate to about 8 decimal places.... it really
    > shouldn't stack that much at all. *i think*


    I've measured it with a non-trivial number like '12345.6789' and it
    took 3669 iterations. That's quite a bit of stack pressure.

    >
    > but yes, maybe a for loop would be more efficient...but... oh well.
    >


    Here's my entry :D (based on your formula):

    double my_sqrt( const double & x )
    {
    double y = 0;
    double z = 1;

    while ( y != z )
    {
    z = y;
    y = y - (y * y - x) / (2 * x);
    }

    return y;
    }

    Cheers,
    Andre
     
    , Dec 2, 2005
    #9
  10. Mark

    mlimber Guest

    Mark P wrote:
    > Mark wrote:
    > > wrote:
    > >
    > >
    > >>I thought it was bad to directly compare floats?

    > >
    > >
    > > i have no idea... never heard that before.
    > >

    >
    > "Bad" may be too strong, but it should be done with care. A computer
    > will compare 2.0 and 1.999999 as unequal; sometimes this is undesirable.
    > Often rather than comparing for equality it makes more sense to see if
    > the absolute value of the difference is less than some small number (say
    > 0.0001).


    See the FAQ:

    http://www.parashift.com/c -faq-lite/newbie.html#faq-29.17

    Cheers! --M
     
    mlimber, Dec 2, 2005
    #10
  11. Mark

    Neil Cerutti Guest

    On 2005-12-02, mlimber <> wrote:
    > Mark P wrote:
    >> Mark wrote:
    >> > wrote:
    >> >
    >> >
    >> >>I thought it was bad to directly compare floats?
    >> >
    >> >
    >> > i have no idea... never heard that before.
    >> >

    >>
    >> "Bad" may be too strong, but it should be done with care. A
    >> computer will compare 2.0 and 1.999999 as unequal; sometimes
    >> this is undesirable. Often rather than comparing for equality
    >> it makes more sense to see if the absolute value of the
    >> difference is less than some small number (say 0.0001).

    >
    > See the FAQ:
    >
    > http://www.parashift.com/c -faq-lite/newbie.html#faq-29.17


    This algorithm will not be helped much by a better equality
    operation. Lots of representable positive real numbers don't
    converge with this algorithm, e.g., 0.05.

    A better algorith is what's needed, or a bunch of asserts.

    --
    Neil Cerutti
     
    Neil Cerutti, Dec 2, 2005
    #11
  12. Mark

    mlimber Guest

    Neil Cerutti wrote:
    > On 2005-12-02, mlimber <> wrote:
    > > Mark P wrote:
    > >> Mark wrote:
    > >> > wrote:
    > >> >
    > >> >
    > >> >>I thought it was bad to directly compare floats?
    > >> >
    > >> >
    > >> > i have no idea... never heard that before.
    > >> >
    > >>
    > >> "Bad" may be too strong, but it should be done with care. A
    > >> computer will compare 2.0 and 1.999999 as unequal; sometimes
    > >> this is undesirable. Often rather than comparing for equality
    > >> it makes more sense to see if the absolute value of the
    > >> difference is less than some small number (say 0.0001).

    > >
    > > See the FAQ:
    > >
    > > http://www.parashift.com/c -faq-lite/newbie.html#faq-29.17

    >
    > This algorithm will not be helped much by a better equality
    > operation.

    [snip]

    Sure, but it still does need a better "equality" operation.

    Cheers! --M
     
    mlimber, Dec 2, 2005
    #12
  13. Mark

    Greg Guest

    wrote:
    > Mark wrote:
    > > always wondered how a computer might go about doing this...
    > > i'm taking calculus at university... and we went over newton's method.
    > > i found this interesting...
    > >
    > > ------------------------------
    > > #include <cstdlib>
    > > #include <iostream>
    > > #include <cmath>
    > >
    > > using namespace std;
    > >
    > > double f1(double x);
    > > double f2(double x, double y);
    > >
    > > int main(int argc, char *argv[])
    > > {
    > > double x = 5;
    > > cout << f1(x) << endl;
    > > cout << sqrt(x) << endl;
    > >
    > > system("PAUSE");
    > > return EXIT_SUCCESS;
    > > }
    > >
    > > double f1(double x)
    > > {
    > > return f2(x,x/2);
    > > }
    > >
    > > double f2(double x, double y)
    > > {
    > > double z = y - (y*y-x)/(2*x);
    > > if( y == z ) return z;

    >
    > I thought it was bad to directly compare floats?


    It would depend on why the floats are being compared. Testing two
    different floating point values for equality is a bad idea since the
    values are compared exactly even though the values may not be
    represented exactly. So two values that would be considered equal by
    the programmer may not always compare as equal in the program.

    In this case, however, the comparison is fine. One floating value is
    being compared against itself, or more precisely, the result of one
    iteration is being compared against the result of the previous
    iteration. Since the difference between the values produced at each
    iteration becomes smaller and smaller, at some point the difference
    will have become so minute that the floating point variable can no
    longer detect it.

    At that point, the two floating point values will compare equal (since
    the new value appears unchanged from the previous one). Since there is
    no point in any further iterations when that happens, the function
    stops iterating and returns with its result.

    Greg
     
    Greg, Dec 3, 2005
    #13
  14. Mark

    Kai-Uwe Bux Guest

    Greg wrote:

    > wrote:
    >> Mark wrote:
    >> > always wondered how a computer might go about doing this...
    >> > i'm taking calculus at university... and we went over newton's method.
    >> > i found this interesting...
    >> >
    >> > ------------------------------
    >> > #include <cstdlib>
    >> > #include <iostream>
    >> > #include <cmath>
    >> >
    >> > using namespace std;
    >> >
    >> > double f1(double x);
    >> > double f2(double x, double y);
    >> >
    >> > int main(int argc, char *argv[])
    >> > {
    >> > double x = 5;
    >> > cout << f1(x) << endl;
    >> > cout << sqrt(x) << endl;
    >> >
    >> > system("PAUSE");
    >> > return EXIT_SUCCESS;
    >> > }
    >> >
    >> > double f1(double x)
    >> > {
    >> > return f2(x,x/2);
    >> > }
    >> >
    >> > double f2(double x, double y)
    >> > {
    >> > double z = y - (y*y-x)/(2*x);
    >> > if( y == z ) return z;

    >>
    >> I thought it was bad to directly compare floats?

    >
    > It would depend on why the floats are being compared. Testing two
    > different floating point values for equality is a bad idea since the
    > values are compared exactly even though the values may not be
    > represented exactly. So two values that would be considered equal by
    > the programmer may not always compare as equal in the program.
    >
    > In this case, however, the comparison is fine. One floating value is
    > being compared against itself, or more precisely, the result of one
    > iteration is being compared against the result of the previous
    > iteration. Since the difference between the values produced at each
    > iteration becomes smaller and smaller, at some point the difference
    > will have become so minute that the floating point variable can no
    > longer detect it.


    This reasoning is based on calculus and does not apply to floating point
    numbers of limitedprecision.

    > At that point, the two floating point values will compare equal (since
    > the new value appears unchanged from the previous one). Since there is
    > no point in any further iterations when that happens, the function
    > stops iterating and returns with its result.


    Nope:

    #include <iostream>
    #include <iomanip>
    #include <limits>
    #include <cmath>

    double my_sqrt( const double & x )
    {
    double y = 0;
    double z = 1;

    while ( y != z )
    {
    z = y;
    y = y - (y * y - x) / (2 * x);
    std::cout << std::setprecision(17)
    << y << " "
    << z << '\n';
    }

    return y;
    }


    double my_sqrt_b ( const double & x )
    {
    double root = 0;
    double quot = x;

    while ( std::abs( root - quot ) >
    std::numeric_limits<double>::epsilon() ) {
    root = ( root + quot ) / 2;
    quot = x / root;
    std::cout << root << " " << quot << '\n';
    }

    return ( root );
    }

    int main ( void ) {
    std::cout << my_sqrt( 0.5 ) << '\n';
    }

    prints on my machine:

    0.5 0
    0.75 0.5
    0.6875 0.75
    0.71484375 0.6875
    0.7038421630859375 0.71484375
    0.70844837254844606 0.7038421630859375
    0.7065492759819042 0.70844837254844606
    0.7073373965913512 0.7065492759819042
    0.70701120397472073 0.7073373965913512
    0.70714636142893661 0.70701120397472073
    0.70709038494675236 0.70714636142893661
    0.70711357246260598 0.70709038494675236
    0.70710396810177689 0.70711357246260598
    0.70710794639649821 0.70710396810177689
    0.70710629853942519 0.70710794639649821
    0.70710698110529846 0.70710629853942519
    0.70710669837744955 0.70710698110529846
    0.70710681548719212 0.70710669837744955
    0.70710676697875419 0.70710681548719212
    0.70710678707160801 0.70710676697875419
    0.70710677874887562 0.70710678707160801
    0.70710678219626433 0.70710677874887562
    0.70710678076830913 0.70710678219626433
    0.70710678135978755 0.70710678076830913
    0.7071067811147892 0.70710678135978755
    0.7071067812162708 0.7071067811147892
    0.70710678117423575 0.7071067812162708
    0.70710678119164727 0.70710678117423575
    0.70710678118443515 0.70710678119164727
    0.70710678118742254 0.70710678118443515
    0.70710678118618508 0.70710678118742254
    0.70710678118669767 0.70710678118618508
    0.70710678118648529 0.70710678118669767
    0.70710678118657333 0.70710678118648529
    0.7071067811865368 0.70710678118657333
    0.70710678118655201 0.7071067811865368
    0.70710678118654569 0.70710678118655201
    0.70710678118654824 0.70710678118654569
    0.70710678118654724 0.70710678118654824
    0.70710678118654768 0.70710678118654724
    0.70710678118654746 0.70710678118654768
    0.70710678118654757 0.70710678118654746
    0.70710678118654746 0.70710678118654757
    0.70710678118654757 0.70710678118654746
    0.70710678118654746 0.70710678118654757
    0.70710678118654757 0.70710678118654746
    0.70710678118654746 0.70710678118654757
    0.70710678118654757 0.70710678118654746
    0.70710678118654746 0.70710678118654757
    0.70710678118654757 0.70710678118654746
    0.70710678118654746 0.70710678118654757
    0.70710678118654757 0.70710678118654746
    0.70710678118654746 0.70710678118654757
    0.70710678118654757 0.70710678118654746
    0.70710678118654746 0.70710678118654757
    0.70710678118654757 0.70710678118654746
    0.70710678118654746 0.70710678118654757

    As you can see, you run the risk that y and z enter a cycle. At this point,
    their difference does not get any smaller.


    I just know enough about numerical analysis to stay away from these
    problems. Dealing with floats and doubles is *hard* and best left to
    experts who can cope with a total break down of Calculus based intuition
    and reasoning.


    Best

    Kai-Uwe Bux
     
    Kai-Uwe Bux, Dec 3, 2005
    #14
    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. tony

    how to solve file.copy problem

    tony, Aug 23, 2004, in forum: ASP .Net
    Replies:
    3
    Views:
    3,844
    Raghavendra T V
    Aug 26, 2004
  2. Replies:
    0
    Views:
    1,266
  3. Rick Osborn
    Replies:
    10
    Views:
    3,982
    Jon A. Cruz
    Feb 8, 2004
  4. hector
    Replies:
    5
    Views:
    428
    CBFalconer
    Dec 5, 2006
  5. Marc Girod

    how to solve when no root privilege?

    Marc Girod, Sep 25, 2010, in forum: Perl Misc
    Replies:
    4
    Views:
    112
Loading...

Share This Page