# 7.1 > with(student); [D, Diff, Doubleint, Int, Limit, Lineint, Product, Sum, Tripleint, changevar, combine, completesquare, distance, equate, extrema, integrand, intercept, intparts, isolate, leftbox, leftsum, makeproc, maximize, middlebox, middlesum, midpoint, minimize, powsubs, rightbox, rightsum, showtangent, simpson, slope, trapezoid, value] > bound71 := (f,a,b,n) -> abs(f(b)-f(a))*(b-a)/n; | f(b) - f(a) | (b - a) bound71 := (f, a, b, n) -> ----------------------- n # #3 Find n so abs(Int(x,x=1..3 )-L[n])<=.000005 > eqn:= bound71(x->x ,1,3,n)=.000005; > round(solve(eqn,n)) ; -5 eqn := 4/n = .5 10 800000 # So, the bound in 7.1 says to use 800,000 subdivisions to be within # .000005 of the correct answer, # if you are using Ln or Rn, 400,000 if you are using Tn. # Lets check the first bound in 7.2 to see if its better. > bound72 := (f,a,b,n,K1 ) -> K1*(b-a)^2/(2*n); 2 K1 (b - a) bound72 := (f, a, b, n, K1) -> 1/2 ----------- n # In order to apply this, we need to get a value for K[1] which needs to # satisfy K[1]>=diff(f(x),x) for # all x in the interval [a,b]. For the function x -> x we can take # K[1] to be 1 or any larger number. > solve(bound72(x->x ,1,3,n,1)= .000005,n); 400000. > # Yes this bound is twice as good for Ln or Rn. # # # 32 (a) Tabulate the right sums Rn approximating (a) # Int(x^2,x=1..2 for n = 2,8,32, and 128. # # Compare the error bound given in 7.1 with the actual error. # # While we are at it, we will compare with the bound given in the next # section also. # We can take K1 to be the derivative at x=2, 2*2. # # Here f is the function x->x^2, a is 1, b is 2, and n runs from 2^1 to # 2^7 by powers of two. printlevel := 1 > header := > `# intervals`,`right rule`,`errbd1`,`errbd2`,`realerr` : > f := x->x^2; > a := 1; b := 2; > m := 7; K1 := 4; > integral := evalf(int(f(x),x=a..b)) ; 2 f := x -> x a := 1 b := 2 m := 7 integral := 2.333333333 > header; for i from 1 to m do > Rn := evalf(rightsum(f(x),x=a..b,2^i)): > errorbound1 := evalf(bound71(f,a,b,2^i )): > errorbound2 := evalf(bound72(f,a,b,2^i,K1)): > actualerror := abs(integral-Rn): > print ( 2^i,Rn,errorbound1,errorbound2,actualerror ); od: # intervals, right rule, errbd1, errbd2, realerr 2, 3.125000000, 1.500000000, 1., .791666667 4, 2.718750000, .7500000000, .5000000000, .385416667 8, 2.523437500, .3750000000, .2500000000, .190104167 16, 2.427734375, .1875000000, .1250000000, .094401042 32, 2.380371094, .09375000000, .06250000000, .047037761 64, 2.356811523, .04687500000, .03125000000, .023478190 128, 2.345062256, .02343750000, .01562500000, .011728923 # Notice that errorbound1 and errorbound2 are never smaller than the # error. The # theorems guarantee this. Also notice that errorbound2 is closer to # the actual error, # at least in this case. We don't have a theorem for this. # # We can check the other problems in the same way. Just copy down the # input cells above and # change the inputs appropriately. # Or we can define a word checkerror to do the rest of the problems by # copying down the # input cell from above and turning it into a procedure by adding a # 'proc' line at the top and # the word 'end' at the bottom. > checkerror := proc(f,a,b,m,K1) > header := > `# intervals`,`right rule`,`errbd1`,`errbd2`,`realerr` : > integral := evalf(int(f(x),x=a..b)): > print(header); for i from 1 to m do > Rn := evalf(rightsum(f(x),x=a..b,2^i)): > errorbound1 := evalf(bound71(f,a,b,2^i )): > errorbound2 := evalf(bound72(f,a,b,2^i,K1)): > actualerror := abs(integral-Rn): > print ( 2^i,Rn,errorbound1,errorbound2,actualerror ) ; od: > end; Warning, `header` is implicitly declared local Warning, `integral` is implicitly declared local Warning, `i` is implicitly declared local Warning, `Rn` is implicitly declared local Warning, `errorbound1` is implicitly declared local Warning, `errorbound2` is implicitly declared local Warning, `actualerror` is implicitly declared local # checkerror for Int(sin(x),x=2..3. Here we can use K1 = 1 (why?) > checkerror(x->sin(x),2,3,7,1); # intervals, right rule, errbd1, errbd2, realerr 2, .3697960761, .3840887094, .2500000000, .2040495840 4, .4748315853, .1920443547, .1250000000, .0990140748 8, .5250871819, .09602217739, .06250000000, .0487584782 16, .5496533048, .04801108867, .03125000000, .0241923553 32, .5617961875, .02400554434, .01562500000, .0120494726 64, .5678325991, .01200277216, .007812500000, .0060130610 128, .5708420483, .006001386084, .003906250000, .0030036118 # Note that again that errorbound1 > errorbound2 > realerror. # # Let's check this over another interval, say 0 to Pi. > checkerror(x->sin(x),0,Pi,7,1); # intervals, right rule, errbd1, errbd2, realerr 2, 1.570796327, 0, 2.467401101, .429203673 4, 1.896118898, 0, 1.233700551, .103881102 8, 1.974231603, 0, .6168502753, .025768397 16, 1.993570344, 0, .3084251376, .006429656 32, 1.998393362, 0, .1542125688, .001606638 64, 1.999598389, 0, .07710628441, .000401611 128, 1.999899601, 0, .03855314220, .000100399 # Whoa! errorbound1 is obviously out to lunch. # Why can't you use the error bound in 7.1 to estimate the error in # using Ln to approximate # Int(sin(x),x=0..Pi)? # Checkerror for Int(sqrt(x),x=1..4). Here we can use K1 = 1/2 (why) > checkerror(x->sqrt(x),1,4,7,.5); # intervals, right rule, errbd1, errbd2, realerr 2, 5.371708245, 1.500000000, 1.125000000, .705041578 4, 5.030092593, .7500000000, .5625000000, .363425926 8, 4.851246679, .3750000000, .2812500000, .184580012 16, 4.759684864, .1875000000, .1406250000, .093018197 32, 4.713358600, .09375000000, .07031250000, .046691933 64, 4.690058391, .04687500000, .03515625000, .023391724 128, 4.678373972, .02343750000, .01757812500, .011707305 > # Same behaviour -- errrorbound1 > errorbound2 > realerror # # Question: # Why can't you use the error bound in 7.2 to estimate the error in # using Ln to approximate # Int(sqrt(x),x=0..1)? > # 32 (b) Do the same using Tn rather than Rn. # # We can change the definition of checkerror appropriately. First copy # down the definition. # The bound in theorem 1 is cut in half for the trapezoid rule. The # bound in theorem 2 # does not apply to the trapezoid rule, but the bound in theorem 3 does # # > bound73 := proc(f,a,b,n,K2) > evalf(K2*(b-a)^3/(12*n^2)) end; bound73 := proc(f, a, b, n, K2) evalf(1/12*K2*(b - a)^3/n^2) end > checkerror := proc(f,a,b,m,K2) > header := > `# intervals`,`trap rule`,`errbd1`,`errbd3`,`realerr` : > integral := evalf(int(f(x),x=a..b)): > print(header); for i from 1 to m do > Tn := evalf(trapezoid(f(x),x=a..b,2^i)): > errorbound1 := evalf(bound71(f,a,b,2^i ))/2: > errorbound3 := evalf(bound73(f,a,b,2^i,K2)) : > actualerror := abs(integral-Tn): > print ( 2^i,Rn,errorbound1,errorbound3,actualerror ) ; od: > end; Warning, `header` is implicitly declared local Warning, `integral` is implicitly declared local Warning, `i` is implicitly declared local Warning, `Tn` is implicitly declared local Warning, `errorbound1` is implicitly declared local Warning, `errorbound3` is implicitly declared local Warning, `actualerror` is implicitly declared local checkerror := proc(f, a, b, m, K2) local header, integral, i, Tn, errorbound1, errorbound3, actualerror; header := `# intervals`, `trap rule`, errbd1, errbd3, realerr; integral := evalf(int(f(x), x = a .. b)); print(header); for i to m do Tn := evalf(trapezoid(f(x), x = a .. b, 2^i)); errorbound1 := 1/2*evalf(bound71(f, a, b, 2^i)); errorbound3 := evalf(bound73(f, a, b, 2^i, K2)); actualerror := abs(integral - Tn); print(2^i, Rn, errorbound1, errorbound3, actualerror) od end > # Checkerror for Int(x^2,x=1..3. Here we can take K2 to be 2 (why?) > checkerror(x->x^2,1,3,7,2); # intervals, trap rule, errbd1, errbd3, realerr 2, 2.345062256, 4.000000000, .3333333333, .333333333 4, 2.345062256, 2.000000000, .08333333333, .083333333 8, 2.345062256, 1.000000000, .02083333333, .020833333 16, 2.345062256, .5000000000, .005208333333, .005208333 32, 2.345062256, .2500000000, .001302083333, .001302083 64, 2.345062256, .1250000000, .0003255208333, .000325521 128, 2.345062256, .06250000000, .00008138020833, .000081380 # We can see that the error bound from theorem 3 is much better than # that from theorem 1. # # Checkerror for Int(sin(x),x=2..3 Here we can take K2 to be 1 (why?) > checkerror(sin,2,3,7,1); # intervals, trap rule, errbd1, errbd3, realerr 2, 2.345062256, .1920443547, .02083333333, .0120052293 4, 2.345062256, .09602217735, .005208333333, .0029918974 8, 2.345062256, .04801108870, .001302083333, .0007473895 16, 2.345062256, .02400554434, .0003255208333, .0001868109 32, 2.345062256, .01200277217, .00008138020833, .0000467005 64, 2.345062256, .006001386080, .00002034505208, .0000116750 -5 -5 128, 2.345062256, .003000693042, .5086263021 10 , .29188 10 # Again theorem 3 does better. # # Checkerror for Int(sqrt(x),x=1..4) Here we can take K2 to be 1/4 # (why?) # # > checkerror(sqrt,1,4,7,1/4); # intervals, trap rule, errbd1, errbd3, realerr 2, 2.345062256, .7500000000, .1406250000, .044958422 4, 2.345062256, .3750000000, .03515625000, .011574074 8, 2.345062256, .1875000000, .008789062500, .002919988 16, 2.345062256, .09375000000, .002197265625, .000731803 32, 2.345062256, .04687500000, .0005493164063, .000183067 64, 2.345062256, .02343750000, .0001373291016, .000045774 128, 2.345062256, .01171875000, .00003433227539, .000011445 # Same behaviour.