Differences between revisions 1 and 36 (spanning 35 versions)
 ⇤ ← Revision 1 as of 2008-05-07 13:19:17 → Size: 11104 Editor: schilly Comment: ← Revision 36 as of 2011-01-12 02:15:49 → ⇥ Size: 45669 Editor: MarshallHampton Comment: Deletions are marked like this. Additions are marked like this. Line 2: Line 2: goto [:interact:interact main page] goto [[interact|interact main page]]<>== Root Finding Using Bisection ==by William Stein{{{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("

Double Precision Root Finding Using Bisection

")@interactdef _(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)}}}{{attachment:bisect.png}}== Newton's Method ==Note that there is a more complicated Newton's method below.by William Stein{{{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, iterateshtml("

Double Precision Root Finding Using Newton's Method

")@interactdef _(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)}}}{{attachment:newton.png}} Line 7: Line 93: {{{@interactdef _(q1=(-1,(-3,3)), q2=(-2,(-3,3)),       cmap=['autumn', 'bone', 'cool', 'copper', 'gray', 'hot', 'hsv', {{{@interactdef _(q1=(-1,(-3,3)), q2=(-2,(-3,3)),      cmap=['autumn', 'bone', 'cool', 'copper', 'gray', 'hot', 'hsv', Line 16: Line 103: show(plot3d(f, (x,-2,2), (y,-2,2)), figsize=5, viewer='tachyon')     }}}attachment:mountains.png show(plot3d(f, (x,-2,2), (y,-2,2)), figsize=5, viewer='tachyon')}}}{{attachment:mountains.png}} Line 22: Line 109: Line 38: Line 126: attachment:tangents.png {{attachment:tangents.png}}== Numerical integrals with the midpoint rule ==by Marshall Hampton{{{var('x')@interactdef 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('

Numerical integrals with the midpoint rule

')    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)}}}{{attachment:num_int.png}}== Numerical integrals with various rules ==by Nick Alexander (based on the work of Marshall Hampton){{{# by Nick Alexander (based on the work of Marshall Hampton)var('x')@interactdef 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_maximum_on_interval(func, q*dx + a, q*dx + dx + a)[1]            xs.append(x)        elif endpoint_rule == 'Lower':            x = find_minimum_on_interval(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_minimum_on_interval(func,a,b)[0])    max_y = max(0, find_maximum_on_interval(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))}}}{{attachment:num_int2.png}}== Some polar parametric curves ==by Marshall Hampton.This is not very general, but could be modified to show other families of polar curves.{{{@interactdef para(n1 = slider(1,5,1,default = 2), n2 = slider(1,5,1,default = 3), a1 = slider(1,10,1/10,6/5), a2 = slider(1,10,1/10,6), b = slider(0,2,1/50,0)):    var('t')    html('$r=' + latex(b+sin(a1*t)^n1 + cos(a2*t)^n2)+'$')    p = parametric_plot((cos(t)*(b+sin(a1*t)^n1 + cos(a2*t)^n2), sin(t)*(b+sin(a1*t)^n1 + cos(a2*t)^n2)), (t,0, 20*pi), plot_points = 1024, rgbcolor = (0,0,0))    show(p, figsize = [5,5], xmin = -2-b, xmax = 2+b, ymin = -2-b, ymax = 2+b, axes = False)}}}{{attachment:polarcurves1.png}} Line 41: Line 242: Enter symbolic functions $f$, $g$, and $a$, a range, then click the appropriatebutton to compute and plot some combination of $f$, $g$, and $a$ along with $f$ and $g$.This is inspired by the Matlab funtool GUI. Enter symbolic functions $f$, $g$, and $a$, a range, then click the appropriate button to compute and plot some combination of $f$, $g$, and $a$ along with $f$ and $g$. This is inspired by the Matlab funtool GUI. Line 47: Line 246: Line 51: Line 249: '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)'], '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)'], Line 55: Line 253: Line 61: Line 258: return return Line 125: Line 322: Line 140: Line 335: attachment:funtool.png {{attachment:funtool.png}} Line 144: Line 338: Line 147: Line 340: This allows user to display the Newton-Raphson procedure one step at a time.It uses the heuristic that, if any of the values of the controls change,then the procedure should be re-started, else it should be continued. This allows user to display the Newton-Raphson procedure one step at a time. It uses the heuristic that, if any of the values of the controls change, then the procedure should be re-started, else it should be continued. Line 155: Line 346: State = Data = None # globals to allow incremental additions to graphics State = Data = None # globals to allow incremental changes in interaction data Line 163: Line 354: step = ['Next','Reset'] ): step = ['Next','Prev', 'Reset'] ): Line 165: Line 356: prange = [xmin,xmax] Line 168: Line 358: X = [RR(x0)] # restart the plot        df = diff(f)        Fplot = plot(f, prange[0], prange[1])        Data = [X, df, Fplot] Data = [ 1 ] # reset the plot Line 173: Line 360: X, df, Fplot = Data    i = len(X) - 1 # compute and append the next x value    xi = X[i]    fi = RR(f(xi))    fpi = RR(df(xi))    xip1 = xi - fi/fpi    X.append(xip1)    msg = xip1s = None # now check x value for reasonableness 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( 'Error: %s!!' % (msg,) )    def Disp( s, color="blue", Trace=Trace ):        Trace.append( """$%s$""" % (color,s,) )    Disp( """f(x) = %s""" % (latex(f),) )    Disp( """f'(x) = %s""" % (latex(df),) )    stop = False Line 184: Line 386: if abs(xip1) > 10E6*(xmax-xmin):        is_inf = True        show_calcs = True        msg = 'Derivative is 0!'        xip1s = latex(xip1.sign()*infinity)        X.pop()    elif not ((xmin - 0.5*(xmax-xmin)) <= xip1 <= (xmax + 0.5*(xmax-xmin))):        show_calcs = True        msg = 'x value out of range; probable divergence!'    if xip1s is None:        xip1s = '%.4g' % (xip1,)    def Disp( s, color="blue" ):        if show_calcs:            html( """$%s$""" % (color,s,) )    Disp( """f(x) = %s""" % (latex(f),) +          """~~~~f'(x) = %s""" % (latex(df),) )    Disp( """i = %d""" % (i,) +          """~~~~x_{%d} = %.4g""" % (i,xi) +          """~~~~f(x_{%d}) = %.4g""" % (i,fi) +          """~~~~f'(x_{%d}) = %.4g""" % (i,fpi) )    if msg:        html( """%s""" % (msg,) )        c = "red"    else:        c = "blue"    Disp( r"""x_{%d} = %.4g - ({%.4g})/({%.4g}) = %s""" % (i+1,xi,fi,fpi,xip1s), color=c )    Fplot += line( [(xi,0),(xi,fi)], linestyle=':', rgbcolor=(1,0,0) ) # vert dotted line    Fplot += points( [(xi,0),(xi,fi)], rgbcolor=(1,0,0) )    labi = text( '\nx%d\n' % (i,), (xi,0), rgbcolor=(1,0,0),                 vertical_alignment="bottom" if fi < 0 else "top" )    if is_inf:        xl = xi - 0.05*(xmax-xmin)        xr = xi + 0.05*(xmax-xmin)        yl = yr = fi    else:        xl = min(xi,xip1) - 0.02*(xmax-xmin)        xr = max(xi,xip1) + 0.02*(xmax-xmin)        yl = -(xip1-xl)*fpi        yr = (xr-xip1)*fpi        Fplot += points( [(xip1,0)], rgbcolor=(0,0,1) ) # new x value        labi += text( '\nx%d\n' % (i+1,), (xip1,0), rgbcolor=(1,0,0),                 vertical_alignment="bottom" if fi < 0 else "top" )    Fplot += line( [(xl,yl),(xr,yr)], rgbcolor=(1,0,0) ) # tangent    show( Fplot+labi, xmin = prange[0], xmax = prange[1] )    Data = [X, df, Fplot]}}}attachment:newtraph.png 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 )}}}{{attachment:newtraph.png}} Line 241: Line 447: {{{ {{{ Line 244: Line 455: @interactdef trans(x=input_box(u^2-v^2, label="x=",type=SR), \    y=input_box(u*v, label="y=",type=SR), \    t_val=slider(0,10,0.2,6, label="Length of curves"), \    u_percent=slider(0,1,0.05,label="u", default=.7),    v_percent=slider(0,1,0.05,label="v", default=.7),    u_range=input_box(range(-5,5,1), label="u lines"),    v_range=input_box(range(-5,5,1), label="v lines")):    thickness=4    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 g1=sum([parametric_plot((SR(u.subs(u=i))._fast_float_('v'),v.subs(u=i)._fast_float_('v')), t_min,t_max, rgbcolor=(1,0,0)) for i in u_range])    g2=sum([parametric_plot((u.subs(v=i)._fast_float_('u'),SR(v.subs(v=i))._fast_float_('u')), t_min,t_max, rgbcolor=(0,0,1)) for i in v_range])    vline_straight=parametric_plot((SR(u.subs(v=v_val))._fast_float_('u'),SR(v.subs(v=v_val))._fast_float_('u')), t_min,t_max, rgbcolor=(0,0,1), linestyle='-',thickness=thickness)    uline_straight=parametric_plot((SR(u.subs(u=u_val))._fast_float_('v'),SR(v.subs(u=u_val))._fast_float_('v')), t_min,t_max,rgbcolor=(1,0,0), linestyle='-',thickness=thickness) from functools import partial@interactdef trans(x=input_box(u2-v2, label="x=",type=SR), \         y=input_box(u*v+cos(u*v), label="y=",type=SR), \         t_val=slider(0,10,0.2,6, label="Length of curves"), \      u_percent=slider(0,1,0.05,label="u", default=.7),         v_percent=slider(0,1,0.05,label="v", default=.7),         u_range=input_box(range(-5,5,1), label="u lines"),         v_range=input_box(range(-5,5,1), label="v lines")):     thickness=4     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  g1=sum([parametric_plot((i,v), t_min,t_max, rgbcolor=(1,0,0)) for i in u_range])     g2=sum([parametric_plot((u,i), t_min,t_max, rgbcolor=(0,0,1)) for i in v_range])     vline_straight=parametric_plot((u,v_val), t_min,t_max, rgbcolor=(0,0,1), linestyle='-',thickness=thickness)     uline_straight=parametric_plot((u_val, v), t_min,t_max,rgbcolor=(1,0,0), linestyle='-',thickness=thickness) Line 264: Line 475: g3=sum([parametric_plot((x.subs(u=i)._fast_float_('v'),y.subs(u=i)._fast_float_('v')), t_min,t_max, rgbcolor=(1,0,0)) for i in u_range])    g4=sum([parametric_plot((x.subs(v=i)._fast_float_('u'),y.subs(v=i)._fast_float_('u')), t_min,t_max, rgbcolor=(0,0,1)) for i in v_range])    vline=parametric_plot((SR(x.subs(v=v_val))._fast_float_('u'),SR(y.subs(v=v_val))._fast_float_('u')), t_min,t_max, rgbcolor=(0,0,1), linestyle='-',thickness=thickness)    uline=parametric_plot((SR(x.subs(u=u_val))._fast_float_('v'),SR(y.subs(u=u_val))._fast_float_('v')), t_min,t_max,rgbcolor=(1,0,0), linestyle='-',thickness=thickness)    (g3+g4+vline+uline).save("xy_coord.png", aspect_ratio=1, figsize=[5,5], axes_labels=['$x$','$y$'])    print jsmath("x=%s, \: y=%s"%(latex(x), latex(y)))    print "
"}}}attachment:coordinate-transform-1.pngattachment:coordinate-transform-2.png xuv = fast_float(x,'u','v')     yuv = fast_float(y,'u','v')     xvu = fast_float(x,'v','u')     yvu = fast_float(y,'v','u')     g3=sum([parametric_plot((partial(xuv,i),partial(yuv,i)), t_min,t_max, rgbcolor=(1,0,0)) for i in u_range])     g4=sum([parametric_plot((partial(xvu,i),partial(yvu,i)), t_min,t_max, rgbcolor=(0,0,1)) for i in v_range])     vline=parametric_plot((partial(xvu,v_val),partial(yvu,v_val)), t_min,t_max, rgbcolor=(0,0,1), linestyle='-',thickness=thickness)     uline=parametric_plot((partial(xuv,u_val),partial(yuv,u_val)), t_min,t_max,rgbcolor=(1,0,0), linestyle='-',thickness=thickness)     (g3+g4+vline+uline).save("xy_coord.png", aspect_ratio=1, figsize=[5,5], axes_labels=['$x$','$y$'])     print jsmath("x=%s, \: y=%s"%(latex(x), latex(y)))     print "
"}}}{{attachment:coordinate-transform-1.png}} {{attachment:coordinate-transform-2.png}} Line 282: Line 493: Line 296: Line 508: attachment:taylor_series_animated.gif {{attachment:taylor_series_animated.gif}}== Illustration of the precise definition of a limit ==by John PerryI'll break tradition and put the image first. Apologies if this is Not A Good Thing.{{attachment:snapshot_epsilon_delta.png}}{{{html("

Limits: ε-δ

")html("This allows you to estimate which values of δ guarantee that f is within ε units of a limit.")html("
• Modify the value of f to choose a function.
• ")html("
• Modify the value of a to change the x-value where the limit is being estimated.
• ")html("
• Modify the value of L to change your guess of the limit.
• ")html("
• Modify the values of δ and ε to modify the rectangle.
")html("If the blue curve passes through the pink boxes, your values for δ and/or ε are probably wrong.")@interactdef delta_epsilon(f = input_box(default=(x^2-x)/(x-1)), a=input_box(default=1), L = input_box(default=1), delta=input_box(label="δ",default=0.1), epsilon=input_box(label="ε",default=0.1), xm=input_box(label="xmin",default=-1), xM=input_box(label="xmax",default=4)):    f_left_plot = plot(f,xm,a-delta/3,thickness=2)    f_right_plot = plot(f,a+delta/3,xM,thickness=2)    epsilon_line_1 = line([(xm,L-epsilon),(xM,L-epsilon)], rgbcolor=(0.5,0.5,0.5),linestyle='--')    epsilon_line_2 = line([(xm,L+epsilon),(xM,L+epsilon)], rgbcolor=(0.5,0.5,0.5),linestyle='--')    ym = min(f_right_plot.ymin(),f_left_plot.ymin())    yM = max(f_right_plot.ymax(),f_left_plot.ymax())    bad_region_1 = polygon([(a-delta,L+epsilon),(a-delta,yM),(a+delta,yM),(a+delta,L+epsilon)], rgbcolor=(1,0.6,0.6))    bad_region_2 = polygon([(a-delta,L-epsilon),(a-delta,ym),(a+delta,ym),(a+delta,L-epsilon)], rgbcolor=(1,0.6,0.6))    aL_point = point((a,L),rgbcolor=(1,0,0),pointsize=20)    delta_line_1 = line([(a-delta,ym),(a-delta,yM)],rgbcolor=(0.5,0.5,0.5),linestyle='--')    delta_line_2 = line([(a+delta,ym),(a+delta,yM)],rgbcolor=(0.5,0.5,0.5),linestyle='--')    (f_left_plot +f_right_plot +epsilon_line_1 +epsilon_line_2 +delta_line_1 +delta_line_2 +aL_point +bad_region_1 +bad_region_2).show(xmin=xm,xmax=xM)}}}== A graphical illustration of sin(x)/x -> 1 as x-> 0 ==by Wai Yan Pong{{{x=var('x')@interactdef _(x = slider(-7/10,7/10,1/20,1/2)):    html('

A graphical illustration of $\lim_{x -> 0} \sin(x)/x =1$

')    html('Below is the unit circle, so the length of the red line is |sin(x)|')    html('and the length of the blue line is |tan(x)| where x is the length of the arc.')     html('From the picture, we see that |sin(x)| $\le$ |x| $\le$ |tan(x)|.')    html('It follows easily from this that cos(x) $\le$ sin(x)/x $\le$ 1 when x is near 0.')    html('As $\lim_{x ->0} \cos(x) =1$, we conclude that $\lim_{x -> 0} \sin(x)/x =1$.')    if not (x == 0):        pretty_print("sin(x)/x = "+str(sin(float(x))/float(x)))    elif x == 0:        pretty_print("The limit of sin(x)/x as x tends to 0 is 1.")    C=circle((0,0),1, rgbcolor='black')    mvp = (cos(x),sin(x));tpt = (1, tan(x))    p1 = point(mvp, pointsize=30, rgbcolor='red'); p2 = point((1,0), pointsize=30, rgbcolor='red')    line1 = line([(0,0),tpt], rgbcolor='black'); line2 = line([(cos(x),0),mvp], rgbcolor='red')     line3 = line([(0,0),(1,0)], rgbcolor='black'); line4 = line([(1,0),tpt], rgbcolor='blue')    result = C+p1+p2+line1+line2+line3+line4    result.show(aspect_ratio=1, figsize=[3,3], axes=False)}}}{{attachment:sinelimit.png}}== Quadric Surface Plotter ==by Marshall Hampton. This is pretty simple, so I encourage people to spruce it up. In particular, it isn't set up to show all possible types of quadrics.{{{var('x,y,z')quadrics = {'Ellipsoid':x^2+y^2+z^2-1,'Elliptic paraboloid':x^2+y^2-z,'Hyperbolic paraboloid':x^2-y^2-z, '1-Sheeted Hyperboloid':x^2+y^2-z^2-1,'2-Sheeted Hyperboloid':x^2-y^2-z^2-1, 'Cone':x^2+y^2-z^2}@interactdef quads(q = selector(quadrics.keys()), a = slider(0,5,1/2,default = 1)):    f = quadrics[q].subs({x:x*a^(1/2)})    if a==0 or q=='Cone': html('
$'+latex(f)+' \$'+ '(degenerate)
')    else: html('
$'+latex(f)+'$
')    p = implicit_plot3d(f,(x,-2,2),(y,-2,2),(z,-2,2), plot_points = 75)    show(p)}}}{{attachment:quadrics.png}}== The midpoint rule for numerically integrating a function of two variables ==by Marshall Hampton{{{from sage.plot.plot3d.platonic import index_face_setdef 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("

The midpoint rule for a function of two variables

")@interactdef 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)/ny    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('

Numerical approximation: ' + str(num_approx)+'

')    p1 = plot3d(f,(x,x_start,x_end),(y,y_start,y_end))    show(p1+sum(cubs))}}}{{attachment:numint2d.png}}== Gaussian (Legendre) quadrature ==by Jason GroutThe output shows the points evaluated using Gaussian quadrature (using a weight of 1, so using Legendre polynomials). The vertical bars are shaded to represent the relative weights of the points (darker = more weight). The error in the trapezoid, Simpson, and quadrature methods is both printed out and compared through a bar graph. The "Real" error is the error returned from scipy on the definite integral.{{{from scipy.special.orthogonal import p_rootsfrom scipy.integrate import quad, trapz, simpsfrom sage.ext.fast_eval import fast_floatfrom numpy import linspaceshow_weight_graph=False# 'Hermite': {'w': e**(-x**2), 'xmin': -numpy.inf, 'xmax': numpy.inf, 'func': h_roots},# 'Laguerre': {'w': e**(-x), 'xmin': 0, 'xmax': numpy.inf, 'func': l_roots},methods = {'Legendre': {'w': 1, 'xmin': -1, 'xmax': 1, 'func': p_roots},                'Chebyshev': {'w': 1/sqrt(1-x**2), 'xmin': -1, 'xmax': 1, 'func': t_roots},                'Chebyshev2': {'w': sqrt(1-x**2), 'xmin': -1, 'xmax': 1, 'func': u_roots},                'Trapezoid': {'w': 1, 'xmin': -1, 'xmax': 1, 'func': lambda n: (linspace(-1r,1,n), numpy.array([1.0r]+[2.0r]*(n-2)+[1.0r])*1.0r/n)},                'Simpson': {'w': 1, 'xmin': -1, 'xmax': 1, 'func': lambda n: (linspace(-1r,1,n), numpy.array([1.0r]+[4.0r,2.0r]*int((n-3.0r)/2.0r)+[4.0r,1.0r])*2.0r/(3.0r*n))}}var("x")def box(center, height, area,**kwds):    width2 = 1.0*area/height/2.0    return polygon([(center-width2,0),(center+width2,0),(center+width2,height),(center-width2,height)],**kwds)        @interactdef weights(n=slider(1,30,1,default=10),f=input_box(default=3*x+cos(10*x)),show_method=["Legendre", "Chebyshev", "Chebyshev2", "Trapezoid","Simpson"]):    ff = fast_float(f,'x')    method = methods[show_method]    xcoords,w = (method['func'])(int(n))    xmin = method['xmin']    xmax = method['xmax']    plot_min = max(xmin, -10)    plot_max = min(xmax, 10)    scaled_func = f*method['w']    scaled_ff = fast_float(scaled_func)    coords = zip(xcoords,w)    max_weight = max(w)    coords_scaled = zip(xcoords,w/max_weight)    f_graph = plot(scaled_func,plot_min,plot_max)    boxes = sum(box(x,ff(x),w*ff(x),rgbcolor=(0.5,0.5,0.5),alpha=0.3) for x,w in coords)    stems = sum(line([(x,0),(x,scaled_ff(x))],rgbcolor=(1-y,1-y,1-y),thickness=2,markersize=6,alpha=y) for x,y in coords_scaled)    points = sum([point([(x,0),(x,scaled_ff(x))],rgbcolor='black',pointsize=30) for x,_ in coords])    graph = stems+points+f_graph+boxes    if show_weight_graph:        graph += line([(x,y) for x,y in coords_scaled], rgbcolor='green',alpha=0.4)        show(graph,xmin=plot_min,xmax=plot_max)    approximation = sum([w*ff(x) for x,w in coords])    integral,integral_error = scipy.integrate.quad(scaled_ff, xmin, xmax)    x_val = linspace(min(xcoords), max(xcoords),n)    y_val = map(scaled_ff,x_val)    trapezoid = integral-trapz(y_val, x_val)    simpson = integral-simps(y_val, x_val)    html("$$\sum_{i=1}^{i=%s}w_i\left(%s\\right)= %s\\approx %s =\int_{-1}^{1}%s \,dx$$"%(n,latex(f.subs(x="x_i")), approximation, integral, latex(scaled_func)))    error_data = [trapezoid, simpson, integral-approximation,integral_error]    print "Trapezoid: %s, Simpson: %s, \nMethod: %s, Real: %s"%tuple(error_data)    show(bar_chart(error_data,width=1),ymin=min(error_data), ymax=max(error_data))}}}{{attachment:quadrature1.png}}{{attachment:quadrature2.png}}== Vector Calculus, 2-D Motion ==By Rob BeezerA fast_float() version is available in a [[http://buzzard.ups.edu/sage/motion-2d.sws|worksheet]]{{{# 2-D motion and vector calculus# Copyright 2009, Robert A. Beezer# Creative Commons BY-SA 3.0 US## 2009/02/15 Built on Sage 3.3.rc0# 2009/02/17 Improvements from Jason Grout## variable parameter is t# later at a particular value named t0#var('t')## parameter range#start=0stop=2*pi## position vector definition# edit here for new example# example is wide ellipse# adjust x, extents in final show()#position=vector( (4*cos(t), sin(t)) )## graphic of the motion itself#path = parametric_plot( position(t).list(), (t, start, stop), color = "black" )## derivatives of motion, lengths, unit vectors, etc#velocity = derivative( position(t) )acceleration = derivative(velocity(t))speed = velocity.norm()speed_deriv = derivative(speed)tangent = (1/speed)*velocitydT = derivative(tangent(t))normal = (1/dT.norm())*dT## interact section# slider for parameter, 24 settings# checkboxes for various vector displays# computations at one value of parameter, t0#@interactdef _(t0 = slider(float(start), float(stop), float((stop-start)/24), float(start) , label = "Parameter"),      pos_check = ("Position", True),       vel_check = ("Velocity", False),      tan_check = ("Unit Tangent", False),      nor_check = ("Unit Normal", False),      acc_check = ("Acceleration", False),      tancomp_check = ("Tangential Component", False),      norcomp_check = ("Normal Component", False)       ):    #    # location of interest    #    pos_tzero = position(t0)    #    # various scalar quantities at point    #    speed_component = speed(t0)    tangent_component = speed_deriv(t0)    normal_component = sqrt( acceleration(t0).norm()^2 - tangent_component^2 )    curvature = normal_component/speed_component^2    #    # various vectors, mostly as arrows from the point    #    pos = arrow((0,0), pos_tzero, rgbcolor=(0,0,0))    tan = arrow(pos_tzero, pos_tzero + tangent(t0), rgbcolor=(0,1,0) )    vel = arrow(pos_tzero, pos_tzero + velocity(t0), rgbcolor=(0,0.5,0))    nor = arrow(pos_tzero, pos_tzero + normal(t0), rgbcolor=(0.5,0,0))    acc = arrow(pos_tzero, pos_tzero + acceleration(t0), rgbcolor=(1,0,1))    tancomp = arrow(pos_tzero, pos_tzero + tangent_component*tangent(t0), rgbcolor=(1,0,1) )    norcomp = arrow(pos_tzero, pos_tzero + normal_component*normal(t0), rgbcolor=(1,0,1))    #    # accumulate the graphic based on checkboxes    #    picture = path    if pos_check:        picture = picture + pos    if vel_check:        picture = picture + vel    if tan_check:        picture = picture+ tan    if nor_check:        picture = picture + nor    if acc_check:        picture = picture + acc    if tancomp_check:        picture = picture + tancomp    if norcomp_check:        picture = picture + norcomp    #    # print textual info    #    print "Position vector defined as r(t)=", position(t)    print "Speed is ", N(speed(t0))    print "Curvature is ", N(curvature)    #    # show accumulated graphical info    # adjust x-,y- extents to get best plot    #    show(picture, xmin=-4,xmax=4, ymin=-1.5,ymax=1.5,aspect_ratio=1)}}}{{attachment:motion2d.png}}== Vector Calculus, 3-D Motion ==by Rob BeezerAvailable as a [[http://buzzard.ups.edu/sage/motion-d.sws|worksheet]]{{{# 3-D motion and vector calculus# Copyright 2009, Robert A. Beezer# Creative Commons BY-SA 3.0 US### 2009/02/15 Built on Sage 3.3.rc0# 2009/02/17 Improvements from Jason Grout## variable parameter is t# later at a particular value named t0# # un-comment double hash (##) to get# time-consuming torsion computation#var('t')## parameter range#start=-4*pistop=8*pi## position vector definition# edit here for new example# example is wide ellipse# adjust figsize in final show() to get accurate aspect ratio#a=1/(8*pi)c=(3/2)*aposition=vector( (exp(a*t)*cos(t), exp(a*t)*sin(t), exp(c*t)) )## graphic of the motion itself#path = parametric_plot3d( position(t).list(), (t, start, stop), color = "black" )## derivatives of motion, lengths, unit vectors, etc#velocity = derivative( position(t) )acceleration = derivative(velocity(t))speed = velocity.norm()speed_deriv = derivative(speed)tangent = (1/speed)*velocitydT = derivative(tangent(t))normal = (1/dT.norm())*dTbinormal = tangent.cross_product(normal)## dB = derivative(binormal(t))## interact section# slider for parameter, 24 settings# checkboxes for various vector displays# computations at one value of parameter, t0#@interactdef _(t0 = slider(float(start), float(stop), float((stop-start)/24), float(start) , label = "Parameter"),      pos_check = ("Position", True),       vel_check = ("Velocity", False),      tan_check = ("Unit Tangent", False),      nor_check = ("Unit Normal", False),      bin_check = ("Unit Binormal", False),      acc_check = ("Acceleration", False),      tancomp_check = ("Tangential Component", False),      norcomp_check = ("Normal Component", False)       ):    #    # location of interest    #    pos_tzero = position(t0)    #    # various scalar quantities at point    #    speed_component = speed(t0)    tangent_component = speed_deriv(t0)    normal_component = sqrt( acceleration(t0).norm()^2 - tangent_component^2 )    curvature = normal_component/speed_component^2    ## torsion = (1/speed_component)*(dB(t0)).dot_product(normal(t0))    #    # various vectors, mostly as arrows from the point    #    pos = arrow3d((0,0,0), pos_tzero, rgbcolor=(0,0,0))    tan = arrow3d(pos_tzero, pos_tzero + tangent(t0), rgbcolor=(0,1,0) )    vel = arrow3d(pos_tzero, pos_tzero + velocity(t0), rgbcolor=(0,0.5,0))    nor = arrow3d(pos_tzero, pos_tzero + normal(t0), rgbcolor=(0.5,0,0))    bin = arrow3d(pos_tzero, pos_tzero + binormal(t0), rgbcolor=(0,0,0.5))    acc = arrow3d(pos_tzero, pos_tzero + acceleration(t0), rgbcolor=(1,0,1))    tancomp = arrow3d(pos_tzero, pos_tzero + tangent_component*tangent(t0), rgbcolor=(1,0,1) )    norcomp = arrow3d(pos_tzero, pos_tzero + normal_component*normal(t0), rgbcolor=(1,0,1))    #    # accumulate the graphic based on checkboxes    #    picture = path    if pos_check:        picture = picture + pos    if vel_check:        picture = picture + vel    if tan_check:        picture = picture+ tan    if nor_check:        picture = picture + nor    if bin_check:        picture = picture + bin    if acc_check:        picture = picture + acc    if tancomp_check:        picture = picture + tancomp    if norcomp_check:        picture = picture + norcomp    #    # print textual info    #    print "Position vector: r(t)=", position(t)    print "Speed is ", N(speed(t0))    print "Curvature is ", N(curvature)    ## print "Torsion is ", N(torsion)    print    print "Right-click on graphic to zoom to 400%"    print "Drag graphic to rotate"    #    # show accumulated graphical info    #    show(picture, aspect_ratio=[1,1,1])}}}{{attachment:motion3d.png}}== Directional Derivatives ==This interact displays graphically a tangent line to a function, illustrating a directional derivative (the slope of the tangent line).{{{var('x,y,t,z')f(x,y)=sin(x)*cos(y)pif = float(pi)line_thickness=3surface_color='blue'plane_color='purple'line_color='red'tangent_color='green'gradient_color='orange'@interactdef 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)}}}{{attachment:directional derivative.png}}== 3D graph with points and curves ==By Robert MarikThis sagelet is handy when showing local, constrained and absolute maxima and minima in two variables. Available as a [[http://user.mendelu.cz/marik/sage/3Dgraph_with_points.sws|worksheet]]{{{%hide%autox,y, t, u, v =var('x y t u v')INI_func='x^2-2*x+y^2-2*y'INI_box='-1,3.2,-1,3.2'INI_points='(1,1,\'green\'),(3/2,3/2),(0,1),(1,0),(0,0,\'black\'),(3,0,\'black\'),(0,3,\'black\')'INI_curves='(t,0,0,3,\'red\'),(0,t,0,3,\'green\'),(t,3-t,0,3)'@interactdef _(func=input_box(INI_func,label="f(x,y)=",type=str),\  bounds=input_box(INI_box,label="xmin,xmax,ymin,ymax",type=str),\  st_points=input_box(INI_points,\  label="points
(comma separated pairs, optionally with color)", type=str),\  bnd_curves=input_box(INI_curves,label="curves on boundary
(x(t),y(t),tmin,tmax,'opt_color')", type=str),\ show_planes=("Show zero planes", False), show_axes=("Show axes", True),  show_table=("Show table", True)): f=sage_eval('lambda x,y: ' + func) html(r'Function $f(x,y)=%s$ '%latex(f(x,y))) xmin,xmax,ymin,ymax=sage_eval('('+bounds+')') A=plot3d(f(x,y),(x,xmin,xmax),(y,ymin,ymax),opacity=0.5) if not(bool(st_points=='')):     st_p=sage_eval('('+st_points+',)')     html(r'')     for current in range(len(st_p)):         point_color='red'         if bool(len(st_p[current])==3):              point_color=st_p[current][2]         x0=st_p[current][0]         y0=st_p[current][1]         z0=f(x0,y0)         if show_table:              html(r'
show(plot)

## Parametric surface example

by Marshall Hampton

var('u,v')
npi = RDF(pi)
@interact
def viewer(mesh = checkbox(default = False, label = 'Show u,v meshlines'), uc = slider(-2,2,1/10,0, label = '<span style="color:red">Constant u value</span>'), vc = slider(-2,2,1/10,0, label = '<span style="color:green">Constant v value</span>'), functions = input_box([u,v^2,u^2+v])):
f1(u,v) = functions[0]
f2(u,v) = functions[1]
f3(u,v) = functions[2]
surface_plot = parametric_plot3d([f1,f2,f3], (u,-2,2), (v,-2,2), mesh = mesh, opacity = .8)
constant_u = line3d([[f1(uc,q), f2(uc,q), f3(uc,q)] for q in srange(-2,2,.01)], rgbcolor = (1,0,0), thickness = 3)
constant_v = line3d([[f1(q,vc), f2(q,vc), f3(q,vc)] for q in srange(-2,2,.01)], rgbcolor = (0,1,0), thickness = 3)
show(surface_plot + constant_u + constant_v, frame = False)

interact/calculus (last edited 2020-08-11 14:10:09 by kcrisman)