#***************************************************************************** # Copyright (C) 2008 Dean Moore # # Distributed under the terms of the GNU General Public License (GPL) # http://www.gnu.org/licenses/ #***************************************************************************** ################################################################################# # # # What this program does, needed input, boilerplate info & disclaimers: # # # # This program was written for the mathematical package SAGE (by implication, # # using the Python programming language), and animates a sliding derivative # # line to a given one-variable single-valued function, defined as *f* (below). # # # # A picture paints a thousand words, and it's best to look at some sample # # output. In short, it's an animated illustration of a standard Calc I # # concept. # # # # This really acts like two programs, depending on value of parameter *byArcs*: # # # # When I first wrote this I simply divided the interval in an even "Riemannian" # # mesh. Trouble is, where function was "flat" the "derivative" line traveled # # slow, and function was steep, the line traveled fast. # # # # Ignored code for a long time, but decided to fix this, dividing curve into # # intervals based on arc length, not an even mesh of x-values. # # # # Was loathe to completely throw away old code, however: # # # # Setting *byArcs* to *false* & it behaves as before, dividing x-axis into an # # even mesh -- this runs faster, avoiding CPU-intensive calculations of arc # # length. Set *byArcs* true and it divides curve into even (up to certain # # unavoidable error limits due to numerical calculations) arcs, each arc # # (nearly) equal every other arc. # # # # Other info: # # # # Once a few parameters are defined the program is designed to churn without # # further input, doing all needed calculations. # # # # Caveat: don't use unbounded functions, e.g., f(x) = 1/x on (0,1]. No # # functions with unbounded first derivatives, e.g., f(x) = x^(1/3) about x = 0. # # # # Doubtless there are other ways to break the code of which I have not # # thought. If the user finds it amusing, fine with me. # # # # See code for user-input material. # # # # The code assumes *f* "nice" with a "nice" derivative. Where there are a few # # checks, it isn't designed to be idiot-proof. # # # ################################################################################# ################################################################################# # # # Why did I write it? A challenge to learn more, to do it in SAGE. I looked # # at < http://en.wikipedia.org/wiki/Derivative > & thought, "Gee, maybe some # # more lively images. It's 2008." # # # # Written by Dean Moore, original in February, 2008, "arc length fix" in # # November 2008. # # # # Egomania: stepdaughter Katie tells me her teacher used a version of it in her # # HS calculus class, illustrating derivatives. Glad it's coming to good use, # # obvious statements. # # # ################################################################################# # User-defined material: ################################################################################# # Do we evenly space points on x-axis, or "follow the arc," evenly spacing by # # arc length? By arc length looks nice, but going straight along x-axis is # # faster in terms of computer time. # # # # So set *byArcs* to *true* and we divide curve into segments of (nearly) equal # # arc length, set it false & it's a faster computation, but final picture isn't # # as nice. # # # ################################################################################# byArcs = true ################################################################################# # Following parameters apply regardless of value of *byArcs* # ################################################################################# f(x) = x*sin(x**2) + 1 # The function. a = -1 # Bottom x-value b = 3 font_size = 10 function_string = text("$f(x)=x* sin(x^2) + 1$", (2.7,5.5), fontsize=font_size, rgbcolor=(0,0,0)) deriv_string = text("$f\ ^'(x)=sin(x^2) + 2x^2* cos(x^2)$", (2.7,4.8), fontsize=font_size, rgbcolor=(0,0,0)) x_axis = text("x-axis", (3.2,0.4), fontsize=font_size, rgbcolor=(0,0,0)) y_axis = text("y-axis", (-0.5, -1.5), fontsize=font_size, rgbcolor=(0,0,0)) zero_one = text("(0, 0)", (-0.5, 1.4), fontsize=font_size, rgbcolor=(0,0,0)) zero_two = text("(1.3552, 2.3076)", (1.3, 2.7), fontsize=font_size, rgbcolor=(0,0,0)) zero_three = text("(2.1945,-1.1828)", (2.2, -1.5), fontsize=font_size, rgbcolor=(0,0,0)) zero_four = text("(2.8137, 3.8081)", (2.8, 4.1), fontsize=font_size, rgbcolor=(0,0,0)) givenPoints = [0] thickness_of_curve = 0.75 # Thickness of function's curve. color_of_curve = (0, 0, 1) # By name. color_of_pos_deriv = (0,1,0) color_of_neg_deriv = (1,0,0) color_of_zero_deriv = (0,0,0) length_of_deriv_line = 3 # Length of "sliding" derivative line. thickness_deriv_line = 1.5 # Thickness of derivative line. cutoff_slope_for_zero_slope_line = 0.05 # Probability of zero-slope line is near # zero. Make a cut-off for coloring a line # the "zero slope" color. sliding_point_size = 20 # By name; how big is the tracing point. color_of_sliding_point = (0, 0, 0) # A point traces the curve; this is its color. figure_size = 4 # Regulates size of final picture. The bigger, the # more bytes / image, and the slower it runs. delay_between_frames = 0 # Time between "frames" of final "movie." ################################################################################# # Following parameters make sense when *byArcs* = *true* # ################################################################################# arc_length_tolerance =0.1 # Computing lengths tends to involve numeric # integrals. Accuracy of calcuation. arc_tolerance =0.1 # Ditto, rarely get exact arc lengths. Accuracy. distance_between_points_on_arc =0.1 # Distance between points on the arc ################################################################################# # Following makes sense if *byArcs* = *false: # ################################################################################# number_of_plot_points = 100 # Number of points at which the function & derivative # are sampled: the more, the smoother runs the final # "movie" -- and the slower the program runs, the # slower the "movie" & the more bytes to the image. ################################################################################# # End user-defined material. # ################################################################################# # The rest of the program churns without user input. The disinterested user # may safely ignore the rest, unless s/he wishes to modify the code -- or I # made some error I didn't catch. ################################################################################# # Several short utility routines follow. # ################################################################################# ################################################################################# # # # Function *color_of_derivative_line*: routine to color derivative line, # # different if up, down or flat # # # ################################################################################# def color_of_derivative_line(bottom_x, bottom_y, top_x, top_y): slope = (top_y - bottom_y) / (top_x - bottom_x) # Actually getting a perfect zero-slope line is near-zero probability. # Make a cut-off. if (slope > cutoff_slope_for_zero_slope_line): return color_of_pos_deriv if (slope < -cutoff_slope_for_zero_slope_line): return color_of_neg_deriv return color_of_zero_deriv # End short routine *color_of_derivative_line*. ################################################################################# # # # Can't make SAGE's &%*$(!! sign() facility work right # # # ################################################################################# def sign(x): if (x > 0): return 1 if (x < 0): return -1 return 0 ################################################################################# # # # Function fillInZeros: # # # # Always fill in values where derivative is zero, wherever they are. # # # ################################################################################# def fillInZeros(f, f_prime, a, b, xVals, givenPoints, function_values): from copy import copy tempxVals = copy(xVals) i = var('i') numRoots = var('numRoots') numRoots = 0 i = 0 if (i < len(xVals)): # Dumb as it seems, make the check. dfx = f_prime(xVals[i]) # Derivative of f at x, slight abuse of notation. while (i < len(xVals) - 1): i = i + 1 sgn = sign(dfx) if (i < len(xVals)): dfx = f_prime(xVals[i]) # Derivative of f at x, slight abuse of notation. while ( ((sign(dfx) == sgn) or (dfx == 0)) and (i < len(xVals)) ): i = i + 1 if (i < len(xVals)): dfx = f_prime(xVals[i]) if ( ((i givenPoints[i])): tempxVals.insert(j + 1, givenPoints[i]) j = j + 1 function_values[givenPoints[i]] = f(givenPoints[i]) return tempxVals ################################################################################# # # # Function *arc_length_a_to_b*: Arc length of function f from *a* to *b* # # # # Uses standard arc length formula from calculus # # # ################################################################################# def arc_length_a_to_b(f_prime, a, b, tolerance): g(x) = sqrt(1 + (f_prime**2)) return g.nintegral(x, a, b, tolerance)[0] # Returns a tuple, want first element # End routine ################################################################################# # # # Function *fillXandYValues*: Fill in x and y-values values: # # # # Fills in x & y-values needed for calcuation. # # # ################################################################################# def fillXandYValues(byArcs, f, a, b, number_of_plot_points, xVals, givenPoints, f_prime, function_values, arc_length_tolerance = 0.1, distance_between_points_on_arc = 0.1): if (byArcs == false): # Not using arc length -- the easy case. for i in srange(a, b + (b - a)/number_of_plot_points, (b - a)/number_of_plot_points): function_values[i] = f(i) xVals.append(i) else: # Going by arc length. More involved. xVals.append(a) function_values[a] = f(a) bottom = a top = b while (bottom < b): top = b # Each pass we make the most liberal estimate & pare down. Bisection isn't the most efficient # for root finding, but always works and ten digits of accuracy probably doesn't here matter: # What's actual arc length? length = arc_length_a_to_b(f_prime, bottom, top, arc_length_tolerance) # Was our estimate good? while (abs(length - distance_between_points_on_arc) > arc_tolerance): top += (sign(distance_between_points_on_arc - length))*(top - bottom)/2 # Isn't that hard to verify. length = arc_length_a_to_b(f_prime, bottom, top, arc_length_tolerance) # By the time we get here estimate is good if (top <= b): # Be safe. function_values[top] = f(top) xVals.append(top) bottom = top if (top != b): # Append right endpoint if not already there. function_values[b] = f(b) xVals.append(b) return fillInZeros(f, f_prime, a, b, xVals, givenPoints, function_values) # End routine *fillXandYValues* ################################################################################# # Top matter, no nice way to classify. Stuff that needs done: ################################################################################# ################################################################################# # Do a few checks, but code isn't designed to be idiot-proof # ################################################################################# cutoff_slope_for_zero_slope_line = abs(cutoff_slope_for_zero_slope_line) delay_between_frames = abs(delay_between_frames) distance_between_points_on_arc = abs(distance_between_points_on_arc) figure_size = abs(figure_size) sliding_point_size = abs(sliding_point_size) thickness_of_curve = abs(thickness_of_curve) length_of_deriv_line = abs(length_of_deriv_line) thickness_deriv_line = abs(thickness_deriv_line) number_of_plot_points = abs(number_of_plot_points) thickness_of_curve = abs(thickness_of_curve) arc_length_tolerance = abs(arc_length_tolerance) arc_tolerance = abs(arc_tolerance) ################################################################################# # Done with checks # ################################################################################# ################################################################################# # Define a few variable; do a few initializations # ################################################################################# v = [] # This will hold graphic image of function's curve. xVals = [] # Holds independent variable. length_of_deriv_line /= 2 # We draw half below, half # above. User doesn't need to know this. ################################################################################# # We use each function value three times; place them # # in a dictionary for efficiency. # ################################################################################# function_values = {} # Blank dictionary to hold function values. ################################################################################# # Function's derivative # ################################################################################# f_prime = f.derivative() # Derivative of function *f* (above). def deriv_line(x0,x): # A function of two variables: return f_prime(x0)*x + ( f(x0) - f_prime(x0)*x0 ) # Draws a line through (x0, f(x0)), # slope (df/dx)(x0). The derivative # and function are evaluated at # "fixed" *x0*; the "real" independent # variable is *x*. Parens are only # for clarity, y = mx + b form. ################################################################################# # Fill in x and y-values needed for the compuation. # # # # Not trying to be fancy using overloaded function, but some parameters make no # # sense for computation by even "mesh" of x-values. # # # ################################################################################# if (byArcs): xVals = fillXandYValues(byArcs, f, a, b, number_of_plot_points, xVals, givenPoints, f_prime, function_values, arc_length_tolerance, distance_between_points_on_arc) else: xVals = fillXandYValues(byArcs, f, a, b, number_of_plot_points, xVals, givenPoints) ################################################################################# # End top matter. # ################################################################################# # The next code is more involved. We need to make a sliding derivative # line of constant length. While experimenting I tried from, say, one # below & one above on the x-axis, and it looked awful. # Define four dictionaries that will hold coordinates of line segments: bottom_x_points = {} bottom_y_points = {} top_x_points = {} top_y_points = {} # Next loop populates the dictionaries with coordinates of line segments of # constant length *length_of_deriv_line* (well, its original value, twice # its current value): for i in range(0, len(xVals)): temp = f_prime(xVals[i]) theta = RDF(arctan(temp)) # Angle of derivative line. # Use trig & basic calc # to verify next block: c = RDF(cos(theta)) # Both of these are s = RDF(sin(theta)) # used twice, never # outside this block. # Note sign patterns in next: bottom_x_points[xVals[i]] = xVals[i] - length_of_deriv_line*c # Bottom x-value of deriv line bottom_y_points[xVals[i]] = function_values[xVals[i]] - length_of_deriv_line*s # Bottom y-value of deriv line. top_x_points[xVals[i]] = xVals[i] + length_of_deriv_line*c # Top x-value of deriv line top_y_points[xVals[i]] = function_values[xVals[i]] + length_of_deriv_line*s # Top y-value of deriv line. # Done computing line segments that comprise the sliding derivative line. # We find max & min x-values for the graph, so the display looks nice: max_x_val = max(max(top_x_points.values()), max(bottom_x_points.values())) # Not perfect, but if enough points, min_x_val = min(min(top_x_points.values()), min(bottom_x_points.values())) # a decent approximation. # We find max & min y-values for graph: max_y_val = max(max(top_y_points.values()), max(bottom_y_points.values())) # See above documentation of min_y_val = min(min(top_y_points.values()), min(bottom_y_points.values())) # *max_x_val* & *min_x_val*. # Graph the sliding derivative line: sliding_deriv_line = animate(line([(bottom_x_points[xVals[i]], bottom_y_points[xVals[i]]), (top_x_points[xVals[i]], top_y_points[xVals[i]] )], thickness = thickness_deriv_line, rgbcolor = color_of_derivative_line(bottom_x_points[xVals[i]], bottom_y_points[xVals[i]], top_x_points[xVals[i]], top_y_points[xVals[i]])) \ + function_string + deriv_string + x_axis + y_axis + zero_one + zero_two + zero_three + zero_four for i in range(0, len(xVals))) # Graph the point that slides along the graph. sliding_points = animate([point((xVals[i], function_values[xVals[i]]), rgbcolor = color_of_sliding_point, pointsize = sliding_point_size) for i in range(0, len(xVals)) ], xmin = min_x_val, ymin = min_y_val, xmax = max_x_val, ymax = max_y_val, figsize = [figure_size,figure_size] ) # Of note on next: when I first did the graph, it was "wiggly." This # led to a note to sage-support, which led to a ticket, which someone # fixed (thanks!), but the fix won't come out until the next version # of SAGE. So for this code I use a work-around. # Animate the graph: graph = (plot(f, [a, b], thickness = thickness_of_curve, rgbcolor = color_of_curve)) for i in xVals: v.append(graph) # End of i-loop. curve = animate(v) # Must be done. # Now show the final "movie": (curve + sliding_points + sliding_deriv_line).show() # Done with program.