# This worksheet implements the newton package used in Chapter 14 of
# Introduction to Scientific Programming by Joseph L. Zachary.
# 
# Create an empty package.
# 
> newton := table();
# 
# Produces an animation showing how Newton's method converges to a
# solution of the function h, beginning from an initial guess g.  There
# are nframes frames created, and the portion of the x-axis to display
# is governed by domain.
# 
> newton[animateNewton] := proc (h, guess, nframes, domain)
> 
>   local i,                # Loop counter
>         x,                # Symbolic constant
>         plts,             # List of plots that makes up animation
>         g,                # Keeps track of guesses
>         hprime,           # First derivative of h
>         thePlot,          # The plot of h
>         lowRange,         # Range of y-axis
>         highRange,
>         doNewton,         # Local procedures
>         plotTangent;
# 
# Does one iteration of Newton's method by improving the guess g for the
# root of h, where the first derivative of h is hprime.
# 
>   doNewton := (h, hprime, g) -> g - h(g)/hprime(g);
# 
# Produces a plot that contains a line that is tangent to h at point
# tangent.  The extent of the x-axis is specified as domain, and the
# extent of the y-axis is specified as range.  The first derivative of h
# is hprime.
# 
>   plotTangent := proc (h, hprime, tangent, domain, range)
>     local x;
>     plot(hprime(tangent)*x + h(tangent) - hprime(tangent)*tangent,
>          x=domain,
>          range,
>          labels=[``,``]);
>     end;
# 
# Determine the minimum and maximum extent of the range of h.
# 
>   lowRange := minimize(h(x), {x}, {x=evalf(domain)});
>   highRange := maximize(h(x), {x}, {x=evalf(domain)});
# 
# Create a plot corresponding to each repetition of Newton's method.
# 
>   g := guess;
>   hprime := D(h);
>   plts := x$0;
>   thePlot := plot(h(x), x=evalf(domain), lowRange..highRange,
> labels=[``,``]);
> 
>   for i from 1 to nframes do
>     plts := plts, plots[display]({thePlot,
>                                   plotTangent(h, hprime, g,
> evalf(domain), 
>                                               lowRange..highRange)});
>     g := doNewton(h, hprime, g);
>  od;
# 
# Display the plots in sequence to form an animation.
# 
>   plots[display]([thePlot,plts], insequence=true);
> 
>  end;
>