# 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"+1 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 "hprime" is the first\ # derivative of "h".\ \ 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; # -------------------------------------------------------------------------------- >