MAT 142 Day 4

1818 days ago by kcrisman

This worksheet will make it easy for you to calculate various approximations (and check your work on the ones you should do by hand). The first thing to do here is evaluate (by Shift-Enter or clicking "evaluate") the first cell here, which will enable you to do all the calculations you need.  Note that the Simpson is the version from class, where $n$ is the number of subdivisions for midpoint and trapezoid separately.

def left(f,a,b,n): Deltax=(b-a)/n;c=a;est=0 for i in range(n): est+=f(c) c+=Deltax return(float(est*Deltax)) def right(f,a,b,n): Deltax=(b-a)/n;c=a;est=0 for i in range(n): c+=Deltax est+=f(c) return(float(est*Deltax)) def midpoint(f,a,b,n): Deltax=(b-a)/n;c=a+Deltax/2;est=0 for i in range(n): est+=f(c) c+=Deltax return(float(est*Deltax)) def trapezoid(f,a,b,n): Deltax=(b-a)/n;c=a;est=0 for i in range(n): est+=f(c)+f(c+Deltax) c+=Deltax return(float(est*Deltax/2)) def Simpson(f,a,b,n): return((2*midpoint(f,a,b,n)+trapezoid(f,a,b,n))/3)

If you click in the cell below and evaluate it, you will get an interactive thing where you can try different functions and get the approximations (and their errors) for all the methods we have discussed.  Just change the function ($f$), endpoints ($a$ and $b$) or number of intervals ($n$) to try others out!

@interact def _(f=x^5,a=0,b=2,n=(2,range(1,100))): f(x) = f try: value=integral(f,x,a,b).n() except TypeError: value=numerical_integral(f,a,b)[0] print "Left approx is",left(f,a,b,n),"with an error of",value-left(f,a,b,n) print "Right approx is",right(f,a,b,n),"with an error of",value-right(f,a,b,n) print "Midpoint approx is",midpoint(f,a,b,n),"with an error of",value-midpoint(f,a,b,n) print "Trapezoid approx is",trapezoid(f,a,b,n),"with an error of",value-trapezoid(f,a,b,n) print "Simpson approx is",Simpson(f,a,b,n),"with an error of",value-Simpson(f,a,b,n)

Click to the left again to hide and once more to show the dynamic interactive window

This cell will enable you to try seeing what the error is for a specific method, but with many different $n$.  Just type in your function and endpoints as before, but this time put the $n$ you want to try in a list.  To try a different method, replace the line "method=trapezoid" with "method=Simpson" or "method=midpoint".

The sample I did has 2, 4, 8, 16, and 32, because I wanted to see what happened to the error as $n$ doubled.  In this case, the error appears to divide by four each time I double $n$.  (Remember, don't get into too big of $n$ - not because Sage can't do it, but because I am not a master programmer and can't avoid recursion issues.)

@interact def _(f=e^(-x^2),a=-1,b=1,method=['trapezoid','Simpson','midpoint']): f(x) = f n_list=[2,4,8,16,32,64,128] try: value=integral(f,x,a,b).n() except TypeError: value=numerical_integral(f,a,b)[0] for n in n_list: print "The error with n =",n,"and using the",eval(method).func_name,"method is",value-eval(method)(f,a,b,n)

method

Click to the left again to hide and once more to show the dynamic interactive window

The rest of this worksheet is eye candy.

Do you want to know what Simpson's rule looks like?

sage.interacts.calculus.simpson_integration()

Simpson integration

$f(x)=$
# divisions
Integration interval
slider:
keyboard:
Computations form

Simpson integration

$f(x)=$
# divisions
Integration interval
slider:
keyboard:
Computations form

Or the trapezoid rule?

sage.interacts.calculus.trapezoid_integration()

Trapezoid integration

$f(x)=$
# divisions
Integration interval
slider:
keyboard:
Computations form

Trapezoid integration

$f(x)=$
# divisions
Integration interval
slider:
keyboard:
Computations form

The following example is another way to do a lot of different types.  It is due to Nick Alexander and Marshall Hampton.

@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'], 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) 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]) 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''' <div class="math"> \begin{align*} \int_{a}^{b} {f(x) \, dx} & = %s \\\ \sum_{i=1}^{%s} {f(x_i) \, \Delta x} & = %s \\\ & = %s \\\ & = %s . \end{align*} </div> ''' % (numerical_answer, number_of_subdivisions, sum_html, num_html, estimated_answer))