Sage Interact

Taylor approximation

{{{id=7| # adapted from interact by Harald Schilly @interact def Taylor_approx(f=input_box(default=sin(x)*exp(-x),label="function"),x0=input_box(default=0,label="center"),order=(1..12)): p = plot(f,x0-1,x0+5, thickness=2) dot = point((x0,f(x0)),pointsize=80,rgbcolor=(1,0,0)) ft = f.taylor(x,x0,order) pt = plot(ft,x0-1, x0+5, color='green', thickness=2) html('$f(x)\;=\;%s$'%latex(f)) html('$\hat{f}(x;%s)\;=\;%s+\mathcal{O}(x^{%s})$'%(x0,latex(ft),order+1)) show(dot + p + pt, ymin = f.find_local_minimum(x0-1,x0+3)[0]-1, ymax = f.find_local_maximum(x0-1,x0+3)[0]+1) /// }}}

Animated

{{{id=108| # animation similar to the above interact f=sin(x)*exp(-x) x0=0 p = plot(f,-1,5, thickness=2) dot = point((x0,f(x0)),pointsize=80,rgbcolor=(1,0,0)) html('$f(x)\;=\;%s$'%latex(f)) #show([order,taylor(f,x,x0,order)] for order in range(1,7,1)) show(animate([dot+p+plot(taylor(f,x,x0,order),-1, 5,ymin=-.5,ymax=1, color='green', thickness=2) for order in range(1,7,1)]),delay=100) /// }}} {{{id=113| /// }}}

An example by Sunil

{{{id=2| # by Sunil Chetty @interact def compare_secant_tangent(f=input_box(default = x^2+1, label="function"), x0=input_box(default=0), delta_x=slider(0,1,1/100,default=1/2, label="Change in x")): fplot=plot(f,(x,x0-1,x0+1), thickness=2) Pplot=list_plot([(x0,f(x0))], color='red', size=75) Qplot=list_plot([(x0+delta_x,f(x0+delta_x))], color='red', size=50) Df=(f.diff(x))(x0) Tangent_line=plot(Df*(x-x0)+f(x0),(x,x0-1,x0+1), linestyle='--') if (delta_x==0): Secant=Df else: Secant=(f(x0+delta_x)-f(x0))/(delta_x) Secant_line=plot(Secant*(x-x0)+f(x0),(x,x0-1,x0+1), color='red') show(fplot+Pplot+Qplot+Tangent_line+Secant_line) /// }}}

Animated

{{{id=115| #starts the same f=x^2+1 x0=0 fplot=plot(f,(x,x0-1,x0+1), thickness=2) Pplot=list_plot([(x0,f(x0))], color='red', size=75) Df=(f.diff(x))(x0) Tangent_line=plot(Df*(x-x0)+f(x0),(x,x0-1,x0+1), linestyle='--') #new part show(animate([fplot+Pplot+Tangent_line+list_plot([(x0+(1/k),f(x0+(1/k)))], color='red', size=50)+plot(((f(x0+(1/k))-f(x0))*k)*(x-x0)+f(x0),(x,x0-1,x0+1),ymin=0.5,ymax=2,color='red') for k in range(1,20,1)]),delay=20) /// }}} {{{id=103| /// }}}

An example by Bret

{{{id=1| # by Bret Benesh # smaller graphs? side by side? labelled? leftend=-10; rightend=10; @interact def delta_epsilon (f=input_box(default=x*x, label="function f"), a=slider(leftend,rightend,1/10,0), L=input_box(default=0, label="suspected limit L"), epsilon=slider(0,1,1/100,1/2), delta=slider(0,1,1/100,1/2)): fplot=plot(f,leftend,rightend,color='purple') fdplot=plot(f,a-delta,a+delta,color='red') Lemplot=plot(L-epsilon,leftend,rightend,linestyle='--',color='green') Lepplot=plot(L+epsilon,leftend,rightend,linestyle='--',color='green') show(fplot+fdplot+Lemplot+Lepplot+text("Local picture",(0,1.1)),ymin=L-2*epsilon,ymax=L+2*epsilon,xmin=a-2*delta,xmax=a+2*delta, figsize=6) show(fplot+fdplot+Lemplot+Lepplot+text("Global picture",(0,120)),ymin=-10,ymax=100,figsize=3) ///
function f 
suspected limit L 
epsilon 
delta 
}}}

Numerical Integration with various methods

{{{id=3| # by Nick Alexander (based on the work of Marshall Hampton) var('x') @interact def midpoint(f = input_box(default = sin(x^2) + 2, type = SR), interval=range_slider(0, 10, 1, default=(0, 4), label="Interval"), number_of_subdivisions = slider(1, 20, 1, default=4, label="Number of boxes"), endpoint_rule = selector(['Midpoint', 'Left', 'Right', 'Upper', 'Lower'], nrows=1, label="Endpoint rule")): a, b = map(QQ, interval) t = sage.calculus.calculus.var('t') func = fast_callable(f(x=t), RDF, vars=[t]) dx = ZZ(b-a)/ZZ(number_of_subdivisions) xs = [] ys = [] for q in range(number_of_subdivisions): if endpoint_rule == 'Left': xs.append(q*dx + a) elif endpoint_rule == 'Midpoint': xs.append(q*dx + a + dx/2) elif endpoint_rule == 'Right': xs.append(q*dx + a + dx) elif endpoint_rule == 'Upper': x = find_local_maximum(func, q*dx + a, q*dx + dx + a)[1] xs.append(x) elif endpoint_rule == 'Lower': x = find_local_minimum(func, q*dx + a, q*dx + dx + a)[1] xs.append(x) ys = [ func(x) for x in xs ] rects = Graphics() for q in range(number_of_subdivisions): xm = q*dx + dx/2 + a x = xs[q] y = ys[q] rects += line([[xm-dx/2,0],[xm-dx/2,y],[xm+dx/2,y],[xm+dx/2,0]], rgbcolor = (1,0,0)) rects += point((x, y), rgbcolor = (1,0,0)) min_y = min(0, find_local_minimum(func,a,b)[0]) max_y = max(0, find_local_maximum(func,a,b)[0]) # html('

Numerical integrals with the midpoint rule

') show(plot(func,a,b) + rects, xmin = a, xmax = b, ymin = min_y, ymax = max_y) def cap(x): # print only a few digits of precision if x < 1e-4: return 0 return RealField(20)(x) sum_html = "%s \cdot \\left[ %s \\right]" % (dx, ' + '.join([ "f(%s)" % cap(i) for i in xs ])) num_html = "%s \cdot \\left[ %s \\right]" % (dx, ' + '.join([ str(cap(i)) for i in ys ])) numerical_answer = integral_numerical(func,a,b,max_points = 200)[0] estimated_answer = dx * sum([ ys[q] for q in range(number_of_subdivisions)]) html(r'''
\begin{align*} \int_{a}^{b} {f(x) \, dx} & = %s \\\ \sum_{i=1}^{%s} {f(x_i) \, \Delta x} & = %s \\\ & = %s \\\ & = %s . \end{align*}
''' % (numerical_answer, number_of_subdivisions, sum_html, num_html, estimated_answer)) /// }}} {{{id=104| /// }}}

Directional Derivative

{{{id=15| # author unknown var('x,y,t,z') f(x,y)=sin(x)*cos(y) pif = float(pi) line_thickness=3 surface_color='blue' plane_color='purple' line_color='red' tangent_color='green' gradient_color='orange' @interact def myfun(location=input_grid(1, 2, default=[0,0], label = "Location (x,y)", width=2), angle=slider(0, 2*pif, label = "Angle"), show_surface=("Show surface", True)): location3d = vector(location[0]+[0]) location = location3d[0:2] direction3d = vector(RDF, [cos(angle), sin(angle), 0]) direction=direction3d[0:2] cos_angle = math.cos(angle) sin_angle = math.sin(angle) df = f.gradient() direction_vector=line3d([location3d, location3d+direction3d], arrow_head=True, rgbcolor=line_color, thickness=line_thickness) curve_point = (location+t*direction).list() curve = parametric_plot(curve_point+[f(*curve_point)], (t,-3,3),color=line_color,thickness=line_thickness) plane = parametric_plot((cos_angle*x+location[0],sin_angle*x+location[1],t), (x, -3,3), (t,-3,3),opacity=0.8, color=plane_color) pt = point3d(location3d.list(),color='green', size=10) tangent_line = parametric_plot((location[0]+t*cos_angle, location[1]+t*sin_angle, f(*location)+t*df(*location)*(direction)), (t, -3,3), thickness=line_thickness, color=tangent_color) picture3d = direction_vector+curve+plane+pt+tangent_line picture2d = contour_plot(f(x,y), (x,-3,3),(y,-3,3), plot_points=100) picture2d += arrow(location.list(), (location+direction).list()) picture2d += point(location.list(),rgbcolor='green',pointsize=40) if show_surface: picture3d += plot3d(f, (x,-3,3),(y,-3,3),opacity=0.7) dff = df(location[0], location[1]) dff3d = vector(RDF,dff.list()+[0]) picture3d += line3d([location3d, location3d+dff3d], arrow_head=True, rgbcolor=gradient_color, thickness=line_thickness) picture2d += arrow(location.list(), (location+dff).list(), rgbcolor=gradient_color, width=line_thickness) show(picture3d,aspect=[1,1,1], axes=True) show(picture2d, aspect_ratio=1) /// }}} {{{id=105| /// }}}

Visualizing linear transformations

{{{id=11| # by Jason Grout # red vector $v$ is the input, determined by theta and r # blue vector $w$ is the image of $v$ via $A$ @interact def linear_transformation(A=matrix([[1,-1],[-1,1/2]]),theta=slider(0, 2*pi, .1), r=slider(0.1, 2, .1, default=1)): v=vector([r*cos(theta), r*sin(theta)]) w = A*v circles = sum([circle((0,0), radius=i, color='black') for i in [1..2]]) html("$%s %s=%s$"%tuple(map(latex, [A, v.column().n(4), w.column().n(4)]))) show(v.plot(color='red')+w.plot(color='blue')+circles,aspect_ratio=1) /// }}} {{{id=106| /// }}}

Differential Equations: Vector fields and Euler's method

{{{id=16| # by Mike Hansen (tested and updated by William Stein, and later by Dan Drake) x,y = var('x,y') from sage.ext.fast_eval import fast_float @interact def _(f = input_box(default=y), g=input_box(default=-x*y+x^3-x), xmin=input_box(default=-1), xmax=input_box(default=1), ymin=input_box(default=-1), ymax=input_box(default=1), start_x=input_box(default=0.5), start_y=input_box(default=0.5), step_size=(0.01,(0.001, 1.5)), steps=(600,(0, 1400)) ): ff = fast_float(f, 'x', 'y') gg = fast_float(g, 'x', 'y') steps = int(steps) points = [ (start_x, start_y) ] for i in range(steps): xx, yy = points[-1] try: points.append( (xx + step_size * ff(xx,yy), yy + step_size * gg(xx,yy)) ) except (ValueError, ArithmeticError, TypeError): break starting_point = point(points[0], pointsize=50) solution = line(points) vector_field = plot_vector_field( (f,g), (x,xmin,xmax), (y,ymin,ymax) ) result = vector_field + starting_point + solution html(r"$\displaystyle\frac{dx}{dt} = %s$ $ \displaystyle\frac{dy}{dt} = %s$" % (latex(f),latex(g))) result.show(xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax) /// }}} {{{id=107| /// }}}

Graph theory

{{{id=10| # by William Stein @interact def _(graph=['CycleGraph', 'CubeGraph', 'RandomGNP'], n=selector([1..10],nrows=1), p=selector([10,20,..,100],nrows=1)): print graph if graph == 'CycleGraph': print "n = %s (number of vertices)"%n G = graphs.CycleGraph(n) elif graph == 'CubeGraph': if n > 8: print "n reduced to 8" n = 8 print "n = %s (dimension)"%n G = graphs.CubeGraph(n) elif graph == 'RandomGNP': print "n = %s (number of vertices)"%n print "p = %s%% (edge probability)"%p G = graphs.RandomGNP(n, p/100.0) print G.automorphism_group() show(plot(G)) /// }}} {{{id=9| /// }}}

Adding point on an elliptic curve

{{{id=91| # by David Møller Hansen def point_txt(P,name,rgbcolor): if (P.xy()[1]) < 0: r = text(name,[float(P.xy()[0]),float(P.xy()[1])-1],rgbcolor=rgbcolor) elif P.xy()[1] == 0: r = text(name,[float(P.xy()[0]),float(P.xy()[1])+1],rgbcolor=rgbcolor) else: r = text(name,[float(P.xy()[0]),float(P.xy()[1])+1],rgbcolor=rgbcolor) return r E = EllipticCurve('37a') list_of_points = E.integral_points() html("Graphical addition of two points $P$ and $Q$ on the curve $ E: %s $"%latex(E)) def line_from_curve_points(E,P,Q,style='-',rgb=(1,0,0),length=25): """ P,Q two points on an elliptic curve. Output is a graphic representation of the straight line intersecting with P,Q. """ # The function tangent to P=Q on E if P == Q: if P[2]==0: return line([(1,-length),(1,length)],linestyle=style,rgbcolor=rgb) else: # Compute slope of the curve E in P l=-(3*P[0]^2 + 2*E.a2()*P[0] + E.a4() - E.a1()*P[1])/((-2)*P[1] - E.a1()*P[0] - E.a3()) f(x) = l * (x - P[0]) + P[1] return plot(f(x),-length,length,linestyle=style,rgbcolor=rgb) # Trivial case of P != R where P=O or R=O then we get the vertical line from the other point elif P[2] == 0: return line([(Q[0],-length),(Q[0],length)],linestyle=style,rgbcolor=rgb) elif Q[2] == 0: return line([(P[0],-length),(P[0],length)],linestyle=style,rgbcolor=rgb) # Non trivial case where P != R else: # Case where x_1 = x_2 return vertical line evaluated in Q if P[0] == Q[0]: return line([(P[0],-length),(P[0],length)],linestyle=style,rgbcolor=rgb) #Case where x_1 != x_2 return line trough P,R evaluated in Q" l=(Q[1]-P[1])/(Q[0]-P[0]) f(x) = l * (x - P[0]) + P[1] return plot(f(x),-length,length,linestyle=style,rgbcolor=rgb) @interact def _(P=selector(list_of_points,label='Point P'),Q=selector(list_of_points,label='Point Q'), marked_points = checkbox(default=True,label = 'Points'), Lines = selector([0..2],nrows=1), Axes=True): curve = E.plot(rgbcolor = (0,0,1),xmin=-5,xmax=5,plot_points=300) R = P + Q Rneg = -R l1 = line_from_curve_points(E,P,Q) l2 = line_from_curve_points(E,R,Rneg,style='--') p1 = plot(P,rgbcolor=(1,0,0),pointsize=40) p2 = plot(Q,rgbcolor=(1,0,0),pointsize=40) p3 = plot(R,rgbcolor=(1,0,0),pointsize=40) p4 = plot(Rneg,rgbcolor=(1,0,0),pointsize=40) textp1 = point_txt(P,"$P$",rgbcolor=(0,0,0)) textp2 = point_txt(Q,"$Q$",rgbcolor=(0,0,0)) textp3 = point_txt(R,"$P+Q$",rgbcolor=(0,0,0)) if Lines==0: g=curve elif Lines ==1: g=curve+l1 elif Lines == 2: g=curve+l1+l2 if marked_points: g=g+p1+p2+p3+p4 if P != Q: g=g+textp1+textp2+textp3 else: g=g+textp1+textp3 g.axes_range(xmin=-5,xmax=5,ymin=-13,ymax=13) show(g,axes = Axes) /// }}} {{{id=101| /// }}}

Coordinate Transformation

{{{id=114| # by Jason Grout var('u v') # polar coordinates #(x,y)=(u*cos(v),u*sin(v)); (u_range,v_range)=([0..6],[0..2*pi,step=pi/12]) # weird example (x,y)=(u^2-v^2,u*v+cos(u*v)); (u_range,v_range)=([-5..5],[-5..5]) thickness=4 square_length=.05 from sage.ext.fast_eval import fast_float from functools import partial @interact def trans(x=input_box(x, label="x",type=SR), y=input_box(y, label="y",type=SR), u_percent=slider(0,1,0.05,label="u", default=.7), v_percent=slider(0,1,0.05,label="v", default=.7), t_val=slider(0,10,0.2,6, label="Length"), u_range=input_box(u_range, label="u lines"), v_range=input_box(v_range, label="v lines")): x(u,v)=x y(u,v)=y u_val = min(u_range)+(max(u_range)-min(u_range))*u_percent v_val = min(v_range)+(max(v_range)-min(v_range))*v_percent t_min = -t_val t_max = t_val uvplot=sum([parametric_plot((i,v), (v,t_min,t_max), color='red',axes_labels=['u','v'],figsize=[5,5]) for i in u_range]) uvplot+=sum([parametric_plot((u,i), (u,t_min,t_max), color='blue',axes_labels=['u','v']) for i in v_range]) uvplot+=parametric_plot((u,v_val), (u,t_min,t_max), rgbcolor=(0,0,1), linestyle='-',thickness=thickness) uvplot+=parametric_plot((u_val, v), (v,t_min,t_max),rgbcolor=(1,0,0), linestyle='-',thickness=thickness) pt=vector([u_val,v_val]) du=vector([(t_max-t_min)*square_length,0]) dv=vector([0,(t_max-t_min)*square_length]) uvplot+=polygon([pt,pt+dv,pt+du+dv,pt+du],color='purple',alpha=0.7) uvplot+=line([pt,pt+dv,pt+du+dv,pt+du],color='green') T(u,v)=(x,y) xuv = fast_float(x,'u','v') yuv = fast_float(y,'u','v') xvu = fast_float(x,'v','u') yvu = fast_float(y,'v','u') xyplot=sum([parametric_plot((partial(xuv,i),partial(yuv,i)), (v,t_min,t_max), color='red', axes_labels=['x','y'],figsize=[5,5]) for i in u_range]) xyplot+=sum([parametric_plot((partial(xvu,i),partial(yvu,i)), (u,t_min,t_max), color='blue') for i in v_range]) xyplot+=parametric_plot((partial(xuv,u_val),partial(yuv,u_val)),(v,t_min,t_max),color='red', linestyle='-',thickness=thickness) xyplot+=parametric_plot((partial(xvu,v_val),partial(yvu,v_val)), (u,t_min,t_max), color='blue', linestyle='-',thickness=thickness) jacobian=abs(T.diff().det()).simplify_full() t_vals=[0..1,step=t_val*.01] vertices=[(x(*c),y(*c)) for c in [pt+t*dv for t in t_vals]] vertices+=[(x(*c),y(*c)) for c in [pt+dv+t*du for t in t_vals]] vertices+=[(x(*c),y(*c)) for c in [pt+(1-t)*dv+du for t in t_vals]] vertices+=[(x(*c),y(*c)) for c in [pt+(1-t)*du for t in t_vals]] xyplot+=polygon(vertices,color='purple',alpha=0.7) xyplot+=line(vertices,color='green') html("$T(u,v)=%s$"%(latex(T(u,v)))) html("Jacobian: $%s$"%latex(jacobian(u,v))) html("A very small region in $xy$ plane is approximately %0.4g times the size of the corresponding region in the $uv$ plane"%jacobian(u_val,v_val).n()) html.table([[uvplot,xyplot]]) /// }}}

Differential Equations: One-variable Euler's method (needs to be fixed)

{{{id=18| var('x y') @interact def euler_method(y_exact_in = input_box('-cos(x)+1.0', type = str, label = 'Exact solution = '), y_prime_in = input_box('sin(x)', type = str, label = "y' = "), start = input_box(0.0, label = 'x starting value: '), stop = input_box(6.0, label = 'x stopping value: '), startval = input_box(0.0, label = 'y starting value: '), nsteps = slider([2^m for m in range(0,10)], default = 10, label = 'Number of steps: '), show_steps = slider([2^m for m in range(0,10)], default = 8, label = 'Number of steps shown in table: ')): y_exact = lambda x: eval(y_exact_in) y_prime = lambda x,y: eval(y_prime_in) stepsize = float((stop-start)/nsteps) steps_shown = min(nsteps,show_steps) sol = [startval] xvals = [start] for step in range(nsteps): sol.append(sol[-1] + stepsize*y_prime(xvals[-1],sol[-1])) xvals.append(xvals[-1] + stepsize) sol_max = max(sol + [find_maximum_on_interval(y_exact,start,stop)[0]]) sol_min = min(sol + [find_minimum_on_interval(y_exact,start,stop)[0]]) show(plot(y_exact(x),start,stop,rgbcolor=(1,0,0))+line([[xvals[index],sol[index]] for index in range(len(sol))]),xmin=start,xmax = stop, ymax = sol_max, ymin = sol_min) if nsteps < steps_shown: table_range = range(len(sol)) else: table_range = range(0,floor(steps_shown/2)) + range(len(sol)-floor(steps_shown/2),len(sol)) html(tab_list([[i,xvals[i],sol[i]] for i in table_range], headers = ['step','x','y'])) /// }}}