SAS Algorithmically(1): Newton-Raphson method

This post was kindly contributed by From a Logical Point of View » SAS - go there to comment and to read the full post.

A good reference for the basic algorithms of Newton-Raphson method to calculate the square root of a number, see

http://mathforum.org/library/drmath/view/52644.html

And the SAS codes(self-explanatory):

data root;
        /*question: find the square root of 4*/
        x=4;

        /*first choose a rough approximation of sqrt(4);
        actually, you can start with any numbers*/
        y0=1;

        count=0;/*init count number*/

        do until (w<1e-8); /*set a small tolerance error*/
                count=count+1;   /*accumulate count number*/
                y=(y0+x/y0)/2;   /*Newton’s formula*/
                w=abs(y-y0); /*if close, exit;*/
                y0=y;        /* otherwise, keep the new one*/
        end;

        output;
run;

The outputs:

x y0 count w y

4 2 6 2.2204E-15 2

After 6 iterations, Newton-Raphson(also called divide-and-average) gets an approximated square root. See what happed during each iteration compared the output generated by SAS function,sqrt():

data root;
        x=4;
        y0=1;
        count=0;
        do until (w<1e-8);
                count=count+1;
                y=(y0+x/y0)/2;
                w=abs(y-y0);
                y0=y;
                if y =sqrt(x) then is_eq_sqrt="YES";
                else is_eq_sqrt="NO";
                output;
        end;
run;

Outputs:

x       y0         count        w               y     is_eq_sqrt

4    2.50000      1      1.50000    2.50000     NO
4    2.05000      2      0.45000    2.05000     NO
4    2.00061      3      0.04939    2.00061     NO
4    2.00000      4      0.00061    2.00000     NO
4    2.00000      5      0.00000    2.00000     NO
4    2.00000      6      0.00000   2.00000     YES

What’s the difference between count 5 and 6 since their y values look the same? We reset the tolerance value to 1e-3 rather than 1e-8, and get the outputs:

x       y0      count       w          y              is_eq_sqrt

4    2.50000      1      1.50000    2.50000      NO
4    2.05000      2      0.45000    2.05000      NO
4    2.00061      3      0.04939    2.00061      NO
4    2.00000      4      0.00061    2.00000      NO

The system get a faster convergence at an higher error rate, with an approximated value little away from sqrt(4).

We should have a deep understanding of how SAS stores numeric values, which deserves a full session to discuss, to unearth the mystery. Some basic references:

This post was kindly contributed by From a Logical Point of View » SAS - go there to comment and to read the full post.

Post navigation