%!PS
%%BoundingBox: 10 40 600 600
/doos { newpath dup 3 2 roll dup 6 1 roll exch moveto
3 2 roll dup 4 1 roll exch lineto
dup 4 1 roll lineto exch lineto closepath stroke } def
%
% 10 40 600 600 doos
%
/Times-Roman findfont 0.13 scalefont setfont
%
0 20 translate 150 150 scale
0.006 setlinewidth
%
% X-axis and Y-axis:
newpath 0.2 0.27 moveto 3.75 0 rlineto
0.77 0.2 moveto 0 3.5 rlineto stroke
%
0.77 0.27 translate
%
% Start & increments:
/x -0.5 def /dx { x /x x 0.01 add def } def
%
/fun { dup 0.5 eq { pop 1000 } { dup 1 exch 2 mul sub dup mul
exch dup mul exch div } ifelse } def
newpath dx dup fun moveto
361 { dx dup fun dup 3.3 gt { moveto }{ lineto } ifelse } repeat stroke
%
0.002 setlinewidth newpath
%
/Times-Roman findfont 0.10 scalefont setfont
%
2.5 -0.08 moveto (x) show
-0.08 1.8 moveto (y) show
-0.08 -0.08 moveto (O) show
0.8 2.5 moveto (y\50x\51 = x^2/\50 1-2.x\51^2) show
%
newpath -0.5 0.25 moveto 3.5 0 rlineto stroke
-0.15 0.27 moveto (1/4) show
newpath 0.5 0 moveto 0 3.5 rlineto stroke
0.45 -0.08 moveto (1/2) show
newpath 0.25 0 moveto 0 0.25 rlineto stroke
0.20 -0.08 moveto (1/4) show
newpath 1 0 moveto 1 1 lineto 0 1 lineto stroke
1.0 -0.08 moveto (1) show
-0.07 1.0 moveto (1) show
%
% Coarsening 1/3 gives 1:
%
newpath 1 3 div 0 moveto 0 1 rlineto stroke
%
% Function y = sqrt(w) with Newton-Rhapson.
%
% Theory: y(n+1) = y(n) - f(y(n))/f'(y(n))
%
% Here: y = sqrt(w) or y^2 = w . Take f(y) = y^2 - w , then f'(y) = 2.y
% ==> y(n+1) = y(n) - (y(n)^2 - w)/(2.y(n)) = y(n) - y(n)/2 + w/y(n)/2
%
% Conclusion: y(n+1) = 1/2.( y(n) + w / y(n) )
% --------------------------------
/wortel
{ /y 0.5 def /w exch def w 0 le { /y 0 def }
{ 7 { /y y w y div add 0.5 mul def } repeat } ifelse y } def
%
/x2 3 2 div 5 wortel 2 div add def
/x3 3 2 div 5 wortel 2 div sub def
%
% Stationary points with multiplicity 2:
%
newpath x2 0 moveto x2 x2 fun lineto 0 x2 fun lineto
x3 0 lineto x3 x3 fun lineto 0 x3 fun lineto x2 0 lineto stroke
showpage