MAT 141, Day 38

2120 days ago by kcrisman

We can evaluate Riemann sums numerically.  Remember, to use this worksheet, you must:

  • Log in to Sage
  • Come back to this worksheet
  • Click "Edit a Copy" at the upper left
  • Evaluate cells by clicking in them and then clicking "evaluate" in the lower left of the cell!

Here, we define all of our stuff from 4.3 # 3b, estimating $\int_0^\pi \sin(x)\; dx$ with $n=8$.  Evaluate this cell!

n=8 f(x)=sin(x) a = 0 b = pi Deltax = (b-a)/n 
       

Now we will write the Riemann sum with left-hand endpoints.  Don't worry about the syntax.  The point is that you can just replace the stuff above with your own functions, and just use this.

My_Sum = sum([Deltax*f(y) for y in srange(a,b,step=Deltax)]) My_Sum.n() 
       
1.97423160194555
1.97423160194555

Here is $\int_1^3 \sqrt{x+1}\, dx$ with $n=20$.  Here, I've kept in an intermediate step so that you see what we're summing.

n = 20 Deltax = (3-1)/n f(x)=sqrt(x+1) My_Sum = sum([Deltax*f(y) for y in srange(1,3,step=Deltax)]) print My_Sum My_Sum.n() 
       
2/5*sqrt(1/5) + 3/10*sqrt(3/10) + 3/10*sqrt(2/5) + 1/5*sqrt(3/5) +
1/10*sqrt(2) + 1/10*sqrt(21/10) + 1/10*sqrt(11/5) + 1/10*sqrt(23/10) +
1/10*sqrt(5/2) + 1/10*sqrt(13/5) + 1/10*sqrt(14/5) + 1/10*sqrt(29/10) +
1/10*sqrt(3) + 1/10*sqrt(31/10) + 1/10*sqrt(33/10) + 1/10*sqrt(17/5) +
1/10*sqrt(7/2) + 1/10*sqrt(37/10) + 1/10*sqrt(19/5) + 1/10*sqrt(39/10)
3.41833964137064
2/5*sqrt(1/5) + 3/10*sqrt(3/10) + 3/10*sqrt(2/5) + 1/5*sqrt(3/5) + 1/10*sqrt(2) + 1/10*sqrt(21/10) + 1/10*sqrt(11/5) + 1/10*sqrt(23/10) + 1/10*sqrt(5/2) + 1/10*sqrt(13/5) + 1/10*sqrt(14/5) + 1/10*sqrt(29/10) + 1/10*sqrt(3) + 1/10*sqrt(31/10) + 1/10*sqrt(33/10) + 1/10*sqrt(17/5) + 1/10*sqrt(7/2) + 1/10*sqrt(37/10) + 1/10*sqrt(19/5) + 1/10*sqrt(39/10)
3.41833964137064

These give left-hand sums.  For a more general Riemann sum, one might have to get a little more hands-on.  Here's a midpoint sum for the same one.

n = 20 Deltax = (3-1)/n f(x)=sqrt(x+1) My_Sum = sum([Deltax*f(y+Deltax/2) for y in srange(1,3,step=Deltax)]) My_Sum.n() 
       
3.44775839078764
3.44775839078764

And here's a random more general sum.  We won't actually do these this semester, though they prove quite useful in technical work.

f(x)=x^2 My_Sum = sum([Deltax*f(y) for y in [0,1,5,6] for Deltax in [1,2,2,1]]) My_Sum.n() 
       
372.000000000000
372.000000000000

This interactive Sagelet from Nick Alexander and Marshall Hampton is very nice to see the difference between different Riemann sums and partitions.  You can use this to help with your homework on any sum with eight or more boxes.  Just click where it says "hide" and evaluate that big cell.

%hide var('x') @interact def midpoint(f = input_box(default = sin(x^2) + 2, type = SR), interval=range_slider(-10, 10, .5, default=(0, 4), label="Interval"), number_of_subdivisions = slider(1, 150, 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('<h3>Numerical integrals with the midpoint rule</h3>') 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*} \text{The most accurate answer is }\int_{a}^{b} {f(x) \, dx} & \approx %s \\\ \text{The approximation you chose is }\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