MAT 142 Day 4

1341 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)) 
       

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