Contact
CoCalc Logo Icon
StoreFeaturesDocsShareSupport News AboutSign UpSign In
| Download

All published worksheets from http://sagenb.org

Views: 168739
Image: ubuntu2004
def bisect_method(f, a, b, eps): try: f = f._fast_float_(f.variables()[0]) except AttributeError: pass intervals = [(a,b)] two = float(2); eps = float(eps) while True: c = (a+b)/two fa = f(a); fb = f(b); fc = f(c) if abs(fc) < eps: return c, intervals if fa*fc < 0: a, b = a, c elif fc*fb < 0: a, b = c, b else: raise ValueError, "f must have a sign change in the interval (%s,%s)"%(a,b) intervals.append((a,b)) html("<h1>Double Precision Root Finding Using Bisection</h1>") @interact def _(f = cos(x) - x, a = float(0), b = float(1), eps=(-3,(-16..-1))): eps = 10^eps print "eps = %s"%float(eps) try: time c, intervals = bisect_method(f, a, b, eps) except ValueError: print "f must have opposite sign at the endpoints of the interval" show(plot(f, a, b, color='red'), xmin=a, xmax=b) else: print "root =", c print "f(c) = %r"%f(c) print "iterations =", len(intervals) P = plot(f, a, b, color='red') h = (P.ymax() - P.ymin())/ (1.5*len(intervals)) L = sum(line([(c,h*i), (d,h*i)]) for i, (c,d) in enumerate(intervals) ) L += sum(line([(c,h*i-h/4), (c,h*i+h/4)]) for i, (c,d) in enumerate(intervals) ) L += sum(line([(d,h*i-h/4), (d,h*i+h/4)]) for i, (c,d) in enumerate(intervals) ) show(P + L, xmin=a, xmax=b)

Double Precision Root Finding Using Bisection

eps 
[removed]
[removed]
[removed]
[removed]
def newton_method(f, c, eps, maxiter=100): x = f.variables()[0] fprime = f.derivative(x) try: g = f._fast_float_(x) gprime = fprime._fast_float_(x) except AttributeError: g = f; gprime = fprime iterates = [c] for i in xrange(maxiter): fc = g(c) if abs(fc) < eps: return c, iterates c = c - fc/gprime(c) iterates.append(c) return c, iterates html("<h1>Double Precision Root Finding Using Newton's Method</h1>") @interact def _(f = x^2 - 2, c = float(0.5), eps=(-3,(-16..-1)), interval=float(0.5)): eps = 10^(eps) print "eps = %s"%float(eps) time z, iterates = newton_method(f, c, eps) print "root =", z print "f(c) = %r"%f(z) n = len(iterates) print "iterations =", n html(iterates) P = plot(f, z-interval, z+interval, rgbcolor='blue') h = P.ymax(); j = P.ymin() L = sum(point((w,(n-1-float(i))/n*h), rgbcolor=(float(i)/n,0.2,0.3), pointsize=10) + \ line([(w,h),(w,j)],rgbcolor='black',thickness=0.2) for i,w in enumerate(iterates)) show(P + L, xmin=z-interval, xmax=z+interval)

Double Precision Root Finding Using Newton's Method

eps 
interval 
[removed]
[removed]
[removed]
[removed]
@interact def _(q1=(-1,(-3,3)), q2=(-2,(-3,3)), cmap=['autumn', 'bone', 'cool', 'copper', 'gray', 'hot', 'hsv', 'jet', 'pink', 'prism', 'spring', 'summer', 'winter']): x,y = var('x,y') f = q1/sqrt((x+1)^2 + y^2) + q2/sqrt((x-1)^2+(y+0.5)^2) C = contour_plot(f, (-2,2), (-2,2), plot_points=30, contours=15, cmap=cmap) show(C, figsize=3, aspect_ratio=1) show(plot3d(f, (x,-2,2), (y,-2,2)), figsize=5, viewer='tachyon')
q1 
q2 
cmap 
[removed]
[removed]
[removed]
[removed]
html('<h2>Tangent line grapher</h2>') @interact def tangent_line(f = input_box(default=sin(x)), xbegin = slider(0,10,1/10,0), xend = slider(0,10,1/10,10), x0 = slider(0, 1, 1/100, 1/2)): prange = [xbegin, xend] x0i = xbegin + x0*(xend-xbegin) var('x') df = diff(f) tanf = f(x0i) + df(x0i)*(x-x0i) fplot = plot(f, prange[0], prange[1]) print 'Tangent line is y = ' + tanf._repr_() tanplot = plot(tanf, prange[0], prange[1], rgbcolor = (1,0,0)) fmax = f.find_maximum_on_interval(prange[0], prange[1])[0] fmin = f.find_minimum_on_interval(prange[0], prange[1])[0] show(fplot + tanplot, xmin = prange[0], xmax = prange[1], ymax = fmax, ymin = fmin)

Tangent line grapher

xbegin 
xend 
x0 
[removed]
[removed]
[removed]
[removed]
var('x') @interact def midpoint(n = slider(1,100,1,4), f = input_box(default = "x^2", type = str), start = input_box(default = "0", type = str), end = input_box(default = "1", type = str)): a = N(start) b = N(end) func = sage_eval(f, locals={'x':x}) dx = (b-a)/n midxs = [q*dx+dx/2 + a for q in range(n)] midys = [func(x_val) for x_val in midxs] rects = Graphics() for q in range(n): xm = midxs[q] ym = midys[q] rects = rects + line([[xm-dx/2,0],[xm-dx/2,ym],[xm+dx/2,ym],[xm+dx/2,0]], rgbcolor = (1,0,0)) + point((xm,ym), rgbcolor = (1,0,0)) min_y = find_minimum_on_interval(func,a,b)[0] max_y = find_maximum_on_interval(func,a,b)[0] html('<h3>Numerical integrals with the midpoint rule</h3>') html('$\int_{a}^{b}{f(x) dx} {\\approx} \sum_i{f(x_i) \Delta x}$') print "\n\nSage numerical answer: " + str(integral_numerical(func,a,b,max_points = 200)[0]) print "Midpoint estimated answer: " + str(RDF(dx*sum([midys[q] for q in range(n)]))) show(plot(func,a,b) + rects, xmin = a, xmax = b, ymin = min_y, ymax = max_y)
start 
end 
[removed]
[removed]
[removed]
[removed]
x = var('x') @interact def _(f=sin(x), g=cos(x), xrange=input_box((0,1)), yrange='auto', a=1, action=selector(['f', 'df/dx', 'int f', 'num f', 'den f', '1/f', 'finv', 'f+a', 'f-a', 'f*a', 'f/a', 'f^a', 'f(x+a)', 'f(x*a)', 'f+g', 'f-g', 'f*g', 'f/g', 'f(g)'], width=15, nrows=5, label="h = "), do_plot = ("Draw Plots", True)): try: f = SR(f); g = SR(g); a = SR(a) except TypeError, msg: print msg[-200:] print "Unable to make sense of f,g, or a as symbolic expressions." return if not (isinstance(xrange, tuple) and len(xrange) == 2): xrange = (0,1) h = 0; lbl = '' if action == 'f': h = f lbl = 'f' elif action == 'df/dx': h = f.derivative(x) lbl = '\\frac{df}{dx}' elif action == 'int f': h = f.integrate(x) lbl = '\\int f dx' elif action == 'num f': h = f.numerator() lbl = '\\text{numer(f)}' elif action == 'den f': h = f.denominator() lbl = '\\text{denom(f)}' elif action == '1/f': h = 1/f lbl = '\\frac{1}{f}' elif action == 'finv': h = solve(f == var('y'), x)[0].rhs() lbl = 'f^{-1}(y)' elif action == 'f+a': h = f+a lbl = 'f + a' elif action == 'f-a': h = f-a lbl = 'f - a' elif action == 'f*a': h = f*a lbl = 'f \\times a' elif action == 'f/a': h = f/a lbl = '\\frac{f}{a}' elif action == 'f^a': h = f^a lbl = 'f^a' elif action == 'f^a': h = f^a lbl = 'f^a' elif action == 'f(x+a)': h = f(x+a) lbl = 'f(x+a)' elif action == 'f(x*a)': h = f(x*a) lbl = 'f(xa)' elif action == 'f+g': h = f+g lbl = 'f + g' elif action == 'f-g': h = f-g lbl = 'f - g' elif action == 'f*g': h = f*g lbl = 'f \\times g' elif action == 'f/g': h = f/g lbl = '\\frac{f}{g}' elif action == 'f(g)': h = f(g) lbl = 'f(g)' html('<center><font color="red">$f = %s$</font></center>'%latex(f)) html('<center><font color="green">$g = %s$</font></center>'%latex(g)) html('<center><font color="blue"><b>$h = %s = %s$</b></font></center>'%(lbl, latex(h))) if do_plot: P = plot(f, xrange, color='red', thickness=2) + \ plot(g, xrange, color='green', thickness=2) + \ plot(h, xrange, color='blue', thickness=2) if yrange == 'auto': show(P, xmin=xrange[0], xmax=xrange[1]) else: yrange = sage_eval(yrange) show(P, xmin=xrange[0], xmax=xrange[1], ymin=yrange[0], ymax=yrange[1])
xrange 
yrange 
h =  
Draw Plots 
[removed]
[removed]
[removed]
[removed]
# ideas from 'A simple tangent line grapher' by Marshall Hampton # http://wiki.sagemath.org/interact State = Data = None # globals to allow incremental changes in interaction data @interact def newtraph(f = input_box(default=8*sin(x)*exp(-x)-1, label='f(x)'), xmin = input_box(default=0), xmax = input_box(default=4*pi), x0 = input_box(default=3, label='x0'), show_calcs = ("Show Calcs",True), step = ['Next','Prev', 'Reset'] ): global State, Data state = [f,xmin,xmax,x0,show_calcs] if (state != State) or (step == 'Reset'): # when any of the controls change Data = [ 1 ] # reset the plot State = state elif step == 'Next': N, = Data Data = [ N+1 ] elif step == 'Prev': N, = Data if N > 1: Data = [ N-1 ] N, = Data df = diff(f) theplot = plot( f, xmin, xmax ) theplot += text( '\n$x_0$', (x0,0), rgbcolor=(1,0,0), vertical_alignment="bottom" if f(x0) < 0 else "top" ) theplot += points( [(x0,0)], rgbcolor=(1,0,0) ) Trace = [] def Err( msg, Trace=Trace ): Trace.append( '<font color="red"><b>Error: %s!!</b></font>' % (msg,) ) def Disp( s, color="blue", Trace=Trace ): Trace.append( """<font color="%s">$ %s $</font>""" % (color,s,) ) Disp( """f(x) = %s""" % (latex(f),) ) Disp( """f'(x) = %s""" % (latex(df),) ) stop = False is_inf = False xi = x0 for i in range(N): fi = RR(f(xi)) fpi = RR(df(xi)) theplot += points( [(xi,fi)], rgbcolor=(1,0,0) ) theplot += line( [(xi,0),(xi,fi)], linestyle=':', rgbcolor=(1,0,0) ) # vert dotted line Disp( """i = %d""" % (i,) ) Disp( """~~~~x_{%d} = %.4g""" % (i,xi) ) Disp( """~~~~f(x_{%d}) = %.4g""" % (i,fi) ) Disp( """~~~~f'(x_{%d}) = %.4g""" % (i,fpi) ) if fpi == 0.0: Err( 'Derivative is 0 at iteration %d' % (i+1,) ) is_inf = True show_calcs = True else: xip1 = xi - fi/fpi Disp( r"""~~~~x_{%d} = %.4g - ({%.4g})/({%.4g}) = %.4g""" % (i+1,xi,fi,fpi,xip1) ) if abs(xip1) > 10*(xmax-xmin): Err( 'Derivative is too close to 0!' ) is_inf = True show_calcs = True elif not ((xmin - 0.5*(xmax-xmin)) <= xip1 <= (xmax + 0.5*(xmax-xmin))): Err( 'x value out of range; probable divergence!' ) stop = True show_calcs = True if is_inf: xl = xi - 0.05*(xmax-xmin) xr = xi + 0.05*(xmax-xmin) yl = yr = fi else: xl = min(xi,xip1) - 0.01*(xmax-xmin) xr = max(xi,xip1) + 0.01*(xmax-xmin) yl = -(xip1-xl)*fpi yr = (xr-xip1)*fpi theplot += text( '\n$x_{%d}$' % (i+1,), (xip1,0), rgbcolor=(1,0,0), vertical_alignment="bottom" if f(xip1) < 0 else "top" ) theplot += points( [(xip1,0)], rgbcolor=(1,0,0) ) theplot += line( [(xl,yl),(xr,yr)], rgbcolor=(1,0,0) ) # tangent if stop or is_inf: break epsa = 100.0*abs((xip1-xi)/xip1) nsf = 2 - log(2.0*epsa)/log(10.0) Disp( r"""~~~~~~~~\epsilon_a = \left|(%.4g - %.4g)/%.4g\right|\times100\%% = %.4g \%%""" % (xip1,xi,xip1,epsa) ) Disp( r"""~~~~~~~~num.~sig.~fig. \approx %.2g""" % (nsf,) ) xi = xip1 show( theplot, xmin=xmin, xmax=xmax ) if show_calcs: for t in Trace: html( t )
f(x) 
xmin 
xmax 
x0 
Show Calcs 
step 
[removed]
[removed]
[removed]
[removed]
var('x') x0 = 0 f = sin(x)*e^(-x) p = plot(f,-1,5, thickness=2) dot = point((x0,f(x0)),pointsize=80,rgbcolor=(1,0,0)) @interact def _(order=(1..12)): ft = f.taylor(x,x0,order) pt = plot(ft,-1, 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 = -.5, ymax = 1) taylor_series_animated.gif
order 
[removed]
[removed]
[removed]
[removed]
Traceback (most recent call last): File "<stdin>", line 1, in <module> File "/home/sage/sagenb/sage_notebook/worksheets/MathBear/6/code/16.py", line 21, in <module> exec compile(ur'taylor_series_animated.gif' + '\n', '', 'single') File "/home/sage/sage_install/sage-a/local/lib/python2.5/site-packages/SQLAlchemy-0.4.6-py2.5.egg/", line 1, in <module> NameError: name 'taylor_series_animated' is not defined
from sage.plot.plot3d.platonic import index_face_set def cuboid(v1,v2,**kwds): """ Cuboid defined by corner points v1 and v2. """ ptlist = [] for vi in (v1,v2): for vj in (v1,v2): for vk in (v1,v2): ptlist.append([vi[0],vj[1],vk[2]]) f_incs = [[0, 2, 6, 4], [0, 1, 3, 2], [0, 1, 5, 4], [1, 3, 7, 5], [2, 3, 7, 6], [4, 5, 7, 6]] if 'aspect_ratio' not in kwds: kwds['aspect_ratio'] = [1,1,1] return index_face_set(f_incs,ptlist,enclosed = True, **kwds) var('x,y') R16 = RealField(16) npi = RDF(pi) sin,cos = math.sin,math.cos html("<h1>The midpoint rule for a function of two variables</h1>") @interact def midpoint2d(func = input_box('y*sin(x)/x+sin(y)',type=str,label='function of x and y'), nx = slider(2,20,1,3,label='x subdivisions'), ny = slider(2,20,1,3,label='y subdivisions'), x_start = slider(-10,10,.1,0), x_end = slider(-10,10,.1,3*npi), y_start= slider(-10,10,.1,0), y_end= slider(-10,10,.1,3*npi)): f = sage_eval('lambda x,y: ' + func) delx = (x_end - x_start)/nx dely = (y_end - y_start)/nx xvals = [RDF(x_start + (i+1.0/2)*delx) for i in range(nx)] yvals = [RDF(y_start + (i+1.0/2)*dely) for i in range(ny)] num_approx = 0 cubs = [] darea = delx*dely for xv in xvals: for yv in yvals: num_approx += f(xv,yv)*darea cubs.append(cuboid([xv-delx/2,yv-dely/2,0],[xv+delx/2,yv+dely/2,f(xv,yv)], opacity = .5, rgbcolor = (1,0,0))) html("$$\int_{"+str(R16(y_start))+"}^{"+str(R16(y_end))+"} "+ "\int_{"+str(R16(x_start))+"}^{"+str(R16(x_end))+"} "+func+"\ dx \ dy$$") html('<p style="text-align: center;">Numerical approximation: ' + str(num_approx)+'</p>') p1 = plot3d(f,(x,x_start,x_end),(y,y_start,y_end)) show(p1+sum(cubs))

The midpoint rule for a function of two variables

function of x and y 
x subdivisions 
y subdivisions 
x_start 
x_end 
y_start 
y_end 
[removed]
[removed]
[removed]
[removed]
npi = RDF(pi) from math import cos,sin def rot(t): return matrix([[cos(t),sin(t)],[-sin(t),cos(t)]]) def pursuit(n,x0,y0,lamb,steps = 100, threshold = .01): paths = [[[x0,y0]]] for i in range(1,n): rx,ry = list(rot(2*npi*i/n)*vector([x0,y0])) paths.append([[rx,ry]]) oldpath = [x[-1] for x in paths] for q in range(steps): diffs = [[oldpath[(j+1)%n][0]-oldpath[j][0],oldpath[(j+1)%n][1]-oldpath[j][1]] for j in range(n)] npath = [[oldpath[j][0]+lamb*diffs[j][0],oldpath[j][1]+lamb*diffs[j][1]] for j in range(n)] for j in range(n): paths[j].append(npath[j]) oldpath = npath return paths html('<h3>Curves of Pursuit</h3>') @interact def curves_of_pursuit(n = slider([2..20],default = 6, label="# of points"),steps = slider([2^i for i in range(1,10)],default = 10, label="# of steps"), stepsize = slider(srange(.01,1,.01),default = .2, label="stepsize"), colorize = checkbox(default = False)): outpaths = pursuit(n,1,0,stepsize, steps = steps) mcolor = (0,0,0) outer = line([q[0] for q in outpaths]+[outpaths[0][0]], rgbcolor = mcolor) if colorize: colors = [hue(j/steps,1,1) for j in range(len(outpaths[0]))] else: colors = [(0,0,0) for j in range(len(outpaths[0]))] nested = sum([line([q[j] for q in outpaths]+[outpaths[0][j]], rgbcolor = colors[j]) for j in range(len(outpaths[0]))]) lpaths = [line(x, rgbcolor = mcolor) for x in outpaths] show(sum(lpaths)+nested, axes = False, figsize = [5,5], xmin = -1, xmax = 1, ymin = -1, ymax =1)

Curves of Pursuit

# of points 
# of steps 
stepsize 
colorize 
[removed]
[removed]
[removed]
[removed]
var('x,y') @interact def example(clr=Color('orange'), f=4*x*exp(-x^2-y^2), xrange='(-2, 2)', yrange='(-2,2)', zrot=(0,pi), xrot=(0,pi), zoom=(1,(1/2,3)), square_aspect=('Square Frame', False), tachyon=('Ray Tracer', True)): xmin, xmax = sage_eval(xrange); ymin, ymax = sage_eval(yrange) P = plot3d(f, (x, xmin, xmax), (y, ymin, ymax), color=clr) html('<h1>Plot of $f(x,y) = %s$</h1>'%latex(f)) aspect_ratio = [1,1,1] if square_aspect else [1,1,1/2] show(P.rotate((0,0,1), -zrot).rotate((1,0,0),xrot), viewer='tachyon' if tachyon else 'jmol', figsize=6, zoom=zoom, frame=False, frame_aspect_ratio=aspect_ratio)
clr 
xrange 
yrange 
zrot 
xrot 
zoom 
Square Frame 
Ray Tracer 
[removed]
[removed]
[removed]
[removed]
var('s,t') g(s) = ((0.57496*sqrt(121 - 16.0*s^2))/sqrt(10.+ s)) def P(color, rng): return parametric_plot3d((cos(t)*g(s), sin(t)*g(s), s), (s,rng[0],rng[1]), (t,0,2*pi), plot_points = [150,150], rgbcolor=color, frame = False, opacity = 1) colorlist = ['red','blue','red','blue'] @interact def _(band_number = selector(range(1,5)), current_color = Color('red')): html('<h1 align=center>Egg Painter</h1>') colorlist[band_number-1] = current_color egg = sum([P(colorlist[i],[-2.75+5.5*(i/4),-2.75+5.5*(i+1)/4]) for i in range(4)]) show(egg)
band_number 
current_color 
[removed]
[removed]
[removed]
[removed]
@interact def color_experimenter(expression=input_box('', 'Expression', str), color=Color('red')): if expression: try: plot(SR(expression), rgbcolor=color).show() except TypeError: print "There's a problem with your expression."
Expression 
color 
[removed]
[removed]
[removed]
[removed]
def error_msg(msg): print '<html><p style="font-family:Arial, sans-serif;color:#000"><span style="color:red;font-weight:bold">Error</span>: %s</p></html>' % msg @interact def interactive_2d_plotter(expression=input_box('sin(x)', 'Expression', str), x_range=range_slider(-10,10,1,(0,10), label='X Range'), square=checkbox(True, 'Square'), axes=checkbox(False, 'Show Axes')): if expression: try: expression = SR(expression) # turn string into a Sage expression except TypeError: print error_msg('This is not an expression.') return try: xmin, xmax = x_range if square or not axes: print "var('%s')\nplot(%s).show(%s%s%s)" % (expression.variables()[0], repr(expression), 'aspect_ratio=1' if square else '', ', ' if square and not axes else '', 'axes=False' if not axes else '') if square: plot(expression, xmin, xmax).show(aspect_ratio=1, axes=axes) else: plot(expression, xmin, xmax).show(axes=axes) else: print "var('%s')\nplot(%s)" % (expression.variables()[0], repr(expression)) plot(expression, xmin, xmax).show(axes=axes) except ValueError: print error_msg('This expression has more than one variable.') return except TypeError: print error_msg("This expression contains an unknown function.") return
Expression 
X Range 
Square 
Show Axes 
[removed]
[removed]
[removed]
[removed]
x,y = var('x,y') @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, 0.2)), steps=(600,(0, 1400)) ): old_f = f f = f.function(x,y) old_g = g g = g.function(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*f(xx,yy), yy+step_size*g(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"<h2>$ \frac{dx}{dt} = %s$ $ \frac{dy}{dt} = %s$</h2>"%(latex(old_f),latex(old_g))) print "Step size: %s"%step_size print "Steps: %s"%steps result.show(xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax)
xmin 
xmax 
ymin 
ymax 
start_x 
start_y 
step_size 
steps 
[removed]
[removed]
[removed]
[removed]
# Solve ODEs using sophisticated Methods like Runga-Kutta-Fehlberg # by Harald Schilly, April 2008 # (jacobian doesn't work, please fix ...) var('x y') @interact def _(fin = input_box(default=y+exp(x/10)-1/3*((x-1/2)^2+y^3)*x-x*y^3), gin=input_box(default=x^3-x+1/100*exp(y*x^2+x*y^2)-0.7*x), xmin=input_box(default=-1), xmax=input_box(default=1.8), ymin=input_box(default=-1.3), ymax=input_box(default=1.5), x_start=(-1,(-2,2)), y_start=(0,(-2,2)), error=(0.5,(0,1)), t_length=(23,(0, 100)) , num_of_points = (1500,(5,2000)), algorithm = selector([ ("rkf45" , "runga-kutta-felhberg (4,5)"), ("rk2" , "embedded runga-kutta (2,3)"), ("rk4" , "4th order classical runga-kutta"), ("rk8pd" , 'runga-kutta prince-dormand (8,9)'), ("rk2imp" , "implicit 2nd order runga-kutta at gaussian points"), ("rk4imp" , "implicit 4th order runga-kutta at gaussian points"), ("bsimp" , "implicit burlisch-stoer (requires jacobian)"), ("gear1" , "M=1 implicit gear"), ("gear2" , "M=2 implicit gear") ])): f(x,y)=fin g(x,y)=gin ff = f._fast_float_(*f.args()) gg = g._fast_float_(*g.args()) #solve path = [] err = error xerr = 0 for yerr in [-err, 0, +err]: T=ode_solver() T.algorithm=algorithm T.function = lambda t, yp: [ff(yp[0],yp[1]), gg(yp[0],yp[1])] T.jacobian = lambda t, yp: [[diff(fun,dval)(yp[0],yp[1]) for dval in [x,y]] for fun in [f,g]] T.ode_solve(y_0=[x_start + xerr, y_start + yerr],t_span=[0,t_length],num_points=num_of_points) path.append(line([p[1] for p in T.solution])) #plot vector_field = plot_vector_field( (f,g), (x,xmin,xmax), (y,ymin,ymax) ) starting_point = point([x_start, y_start], pointsize=50) show(vector_field + starting_point + sum(path), aspect_ratio=1, figsize=[8,9])
fin 
gin 
xmin 
xmax 
ymin 
ymax 
x_start 
y_start 
error 
t_length 
num_of_points 
algorithm 
[removed]
[removed]
[removed]
[removed]
def cobweb(a_function, start, mask = 0, iterations = 20, xmin = 0, xmax = 1): ''' Returns a graphics object of a plot of the function and a cobweb trajectory starting from the value start. INPUT: a_function: a function of one variable start: the starting value of the iteration mask: (optional) the number of initial iterates to ignore iterations: (optional) the number of iterations to draw, following the masked iterations xmin: (optional) the lower end of the plotted interval xmax: (optional) the upper end of the plotted interval EXAMPLES: sage: f = lambda x: 3.9*x*(1-x) sage: show(cobweb(f,.01,iterations=200), xmin = 0, xmax = 1, ymin=0) ''' basic_plot = plot(a_function, xmin = xmin, xmax = xmax) id_plot = plot(lambda x: x, xmin = xmin, xmax = xmax) iter_list = [] current = start for i in range(mask): current = a_function(current) for i in range(iterations): iter_list.append([current,a_function(current)]) current = a_function(current) iter_list.append([current,current]) cobweb = line(iter_list, rgbcolor = (1,0,0)) return basic_plot + id_plot + cobweb var('x') @interact def cobwebber(f_text = input_box(default = "3.8*x*(1-x)",label = "function", type=str), start_val = slider(0,1,.01,.5,label = 'start value'), its = slider([2^i for i in range(0,12)],default = 16, label="iterations")): def f(x): return eval(f_text) show(cobweb(f, start_val, iterations = its))
function 
start value 
iterations 
[removed]
[removed]
[removed]
[removed]