How to Solve a Root...

[Mark]

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]

Re: How to Solve a Root...

(just thought i'd share. no question involved)

[int2str@gmail.com]

Re: How to Solve a Root...


Mark wrote:[color=blue]
> 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;[/color]

I thought it was bad to directly compare floats?
[color=blue]
> else return f2(x,z);
> }
> ---------------------------
>
> will work for cube roots and others if you just tweak the formula.[/color]

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

[gottlobfrege@gmail.com]

Re: How to Solve a Root...


Mark wrote:[color=blue]
> always wondered how a computer might go about doing this...
> i'm taking calculus at university... and we went over newton's method.[/color]
[color=blue]
> will work for cube roots and others if you just tweak the formula.[/color]

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

[Mark]

Re: How to Solve a Root...

int2str@gmail.c om wrote:
[color=blue]
> I thought it was bad to directly compare floats?[/color]

i have no idea... never heard that before.
[color=blue]
> 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.[/color]

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.


gottlobfrege@gm ail.com wrote:[color=blue]
> You might want to look into the Runge-Kutta method. Quite robust. I
> coded it up once many years ago...[/color]

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 P]

Re: How to Solve a Root...

Mark wrote:[color=blue]
> int2str@gmail.c om wrote:
>
>[color=green]
>>I thought it was bad to directly compare floats?[/color]
>
>
> i have no idea... never heard that before.
>[/color]

"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).
[color=blue]
>
> gottlobfrege@gm ail.com wrote:
>[color=green]
>>You might want to look into the Runge-Kutta method. Quite robust. I
>>coded it up once many years ago...[/color]
>[/color]

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

-Mark

[Julián Albo]

Re: How to Solve a Root...

Mark wrote:
[color=blue]
> 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*[/color]

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

[Michiel.Salters@tomtom.com]

Re: How to Solve a Root...


Mark P wrote:
[color=blue]
> Runge-Kutta is a numerical method for solving differential equations.
> Quite different from root finding.[/color]

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.

[int2str@gmail.com]

Re: How to Solve a Root...


Mark wrote:[color=blue]
> int2str@gmail.c om wrote:[color=green]
> > 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.[/color]
>
> 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*[/color]

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.
[color=blue]
>
> but yes, maybe a for loop would be more efficient...but ... oh well.
>[/color]

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

[mlimber]

Re: How to Solve a Root...

Mark P wrote:[color=blue]
> Mark wrote:[color=green]
> > int2str@gmail.c om wrote:
> >
> >[color=darkred]
> >>I thought it was bad to directly compare floats?[/color]
> >
> >
> > i have no idea... never heard that before.
> >[/color]
>
> "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).[/color]

See the FAQ:

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

Cheers! --M

[Neil Cerutti]

Re: How to Solve a Root...

On 2005-12-02, mlimber <mlimber@gmail. com> wrote:[color=blue]
> Mark P wrote:[color=green]
>> Mark wrote:[color=darkred]
>> > int2str@gmail.c om wrote:
>> >
>> >
>> >>I thought it was bad to directly compare floats?
>> >
>> >
>> > i have no idea... never heard that before.
>> >[/color]
>>
>> "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).[/color]
>
> See the FAQ:
>
> http://www.parashift.com/c++-faq-lit...html#faq-29.17[/color]

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

[mlimber]

Re: How to Solve a Root...


Neil Cerutti wrote:[color=blue]
> On 2005-12-02, mlimber <mlimber@gmail. com> wrote:[color=green]
> > Mark P wrote:[color=darkred]
> >> Mark wrote:
> >> > int2str@gmail.c om 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).[/color]
> >
> > See the FAQ:
> >
> > http://www.parashift.com/c++-faq-lit...html#faq-29.17[/color]
>
> This algorithm will not be helped much by a better equality
> operation.[/color]
[snip]

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

Cheers! --M

[Greg]

Re: How to Solve a Root...

int2str@gmail.c om wrote:[color=blue]
> Mark wrote:[color=green]
> > 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;[/color]
>
> I thought it was bad to directly compare floats?[/color]

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

[Kai-Uwe Bux]

Re: How to Solve a Root...

Greg wrote:
[color=blue]
> int2str@gmail.c om wrote:[color=green]
>> Mark wrote:[color=darkred]
>> > 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;[/color]
>>
>> I thought it was bad to directly compare floats?[/color]
>
> 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.[/color]

This reasoning is based on calculus and does not apply to floating point
numbers of limitedprecisio n.
[color=blue]
> 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.[/color]

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::setprecisi on(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_li mits<double>::e psilon() ) {
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.7038421630859 375 0.71484375
0.7084483725484 4606 0.7038421630859 375
0.7065492759819 042 0.7084483725484 4606
0.7073373965913 512 0.7065492759819 042
0.7070112039747 2073 0.7073373965913 512
0.7071463614289 3661 0.7070112039747 2073
0.7070903849467 5236 0.7071463614289 3661
0.7071135724626 0598 0.7070903849467 5236
0.7071039681017 7689 0.7071135724626 0598
0.7071079463964 9821 0.7071039681017 7689
0.7071062985394 2519 0.7071079463964 9821
0.7071069811052 9846 0.7071062985394 2519
0.7071066983774 4955 0.7071069811052 9846
0.7071068154871 9212 0.7071066983774 4955
0.7071067669787 5419 0.7071068154871 9212
0.7071067870716 0801 0.7071067669787 5419
0.7071067787488 7562 0.7071067870716 0801
0.7071067821962 6433 0.7071067787488 7562
0.7071067807683 0913 0.7071067821962 6433
0.7071067813597 8755 0.7071067807683 0913
0.7071067811147 892 0.7071067813597 8755
0.7071067812162 708 0.7071067811147 892
0.7071067811742 3575 0.7071067812162 708
0.7071067811916 4727 0.7071067811742 3575
0.7071067811844 3515 0.7071067811916 4727
0.7071067811874 2254 0.7071067811844 3515
0.7071067811861 8508 0.7071067811874 2254
0.7071067811866 9767 0.7071067811861 8508
0.7071067811864 8529 0.7071067811866 9767
0.7071067811865 7333 0.7071067811864 8529
0.7071067811865 368 0.7071067811865 7333
0.7071067811865 5201 0.7071067811865 368
0.7071067811865 4569 0.7071067811865 5201
0.7071067811865 4824 0.7071067811865 4569
0.7071067811865 4724 0.7071067811865 4824
0.7071067811865 4768 0.7071067811865 4724
0.7071067811865 4746 0.7071067811865 4768
0.7071067811865 4757 0.7071067811865 4746
0.7071067811865 4746 0.7071067811865 4757
0.7071067811865 4757 0.7071067811865 4746
0.7071067811865 4746 0.7071067811865 4757
0.7071067811865 4757 0.7071067811865 4746
0.7071067811865 4746 0.7071067811865 4757
0.7071067811865 4757 0.7071067811865 4746
0.7071067811865 4746 0.7071067811865 4757
0.7071067811865 4757 0.7071067811865 4746
0.7071067811865 4746 0.7071067811865 4757
0.7071067811865 4757 0.7071067811865 4746
0.7071067811865 4746 0.7071067811865 4757
0.7071067811865 4757 0.7071067811865 4746
0.7071067811865 4746 0.7071067811865 4757
0.7071067811865 4757 0.7071067811865 4746
0.7071067811865 4746 0.7071067811865 4757

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