Graphs of Elementary Numerical Integrals

%hide # by Nate Iverson (based on the work of Nick Alexander and Marshall Hampton) var('x') @interact def midpoint(f = input_box(default = sin(x^2) + 2, type = SR,label="f(x)="), interval=input_box(default=[-1,3], label="Interval="), number_of_subdivisions = slider(1, 250, 1, default=4, label="Number of regions="), endpoint_rule = selector(['Midpoint', 'Left', 'Right', 'Upper', 'Lower','Trapezoid','Simpson'], nrows=1, label="Estimation rule"), shade_region=("Shade regions",True)): if shade_region: shadecolor=(.85,.85,1) else: shadecolor=(1,1,1) 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) elif endpoint_rule == 'Trapezoid' or endpoint_rule == 'Simpson': xs.append(q*dx + a) if endpoint_rule == 'Trapezoid' or endpoint_rule == 'Simpson': xs.append(b) estimated_answer=0 ys = [ func(x) for x in xs ] rects = Graphics() poly = Graphics() for q in range(number_of_subdivisions): xm = q*dx + dx/2 + a x = xs[q] y = ys[q] if endpoint_rule == 'Trapezoid': rects += line([[xm-dx/2,0],[xm-dx/2,y],[xm+dx/2,ys[q+1]],[xm+dx/2,0]], rgbcolor = (0,0,0)) rects += point((x, y), rgbcolor = (0,0,0)) poly += polygon([[xm-dx/2,0],[xm-dx/2,y],[xm+dx/2,ys[q+1]],[xm+dx/2,0]], rgbcolor = shadecolor) elif endpoint_rule == 'Simpson': x0 = xm-dx/2 x1 = xm x2 = xm+dx/2 y0 =func(x0) y1 =func(x1) y2 =func(x2) var('z') P=y0*(z-x1)*(z-x2)/((x0-x1)*(x0-x2))+y1*(z-x2)*(z-x0)/((x1-x2)*(x1-x0))+y2*(z-x1)*(z-x0)/((x2-x1)*(x2-x0)) rects += plot(P,xmin=x0,xmax=x2,rgbcolor=(0,0,0), fill="axis", fillcolor=shadecolor) rects += line([[x0,y0],[x0,0]],rgbcolor=(0,0,0)) rects += line([[x2,y2],[x2,0]],rgbcolor=(0,0,0)) rects += point((x0,y0),rgbcolor = (0,0,0)) rects += point((x1,y1),rgbcolor = (0,0,0)) rects += point((x2,y2),rgbcolor = (0,0,0)) estimated_answer +=integral_numerical(P,x0,x2,max_points = 200)[0] else: rects += line([[xm-dx/2,0],[xm-dx/2,y],[xm+dx/2,y],[xm+dx/2,0]], rgbcolor = (0,0,0)) rects += point((x, y), rgbcolor = (0,0,0)) poly +=polygon([[xm-dx/2,0],[xm-dx/2,y],[xm+dx/2,y],[xm+dx/2,0]], rgbcolor = shadecolor) # 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('<h3>Numerical integrals with the midpoint rule</h3>') graph=plot(func,a,b) + rects + poly graph.show(xmin = a, xmax = b) def cap(x): # print only a few digits of precision if abs(x) < 1e-4: return 0 return RealField(20)(x) if endpoint_rule == 'Trapezoid': 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 ])) elif endpoint_rule == 'Simpson': 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 ])) else: 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] if endpoint_rule == 'Trapezoid': estimated_answer = dx/2 * (ys[0]+ys[number_of_subdivisions]+2*sum([ ys[q] for q in range(1,number_of_subdivisions)])) html(r''' <div class="math"> \begin{align*} \int_{a}^{b} {f(x) \, dx} & = %s \\\ T_{%s}(f) & = %s \\\ Error & = %s \end{align*} </div> ''' % (numerical_answer, number_of_subdivisions, estimated_answer,estimated_answer-numerical_answer)) elif endpoint_rule == 'Simpson': html(r''' <div class="math"> \begin{align*} \int_{a}^{b} {f(x) \, dx} & = %s \\\ S_{%s}(f) & = %s \\\ Error & = %s \end{align*} </div> ''' % (numerical_answer, 2*number_of_subdivisions+1, estimated_answer,estimated_answer-numerical_answer)) else: 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 . \\\ Error & = %s \end{align*} </div> ''' % (numerical_answer, number_of_subdivisions, sum_html, num_html, estimated_answer,estimated_answer-numerical_answer)) #