MA 531 HW 4 Sample Code.sagews

#5.10.1.  Create a table of values obtained using the Composite Trapezoidal Rule, with n=1,2,4,8,16 and 32, for int(e^(-x^2)dx, 0,1).  Also output the estimated errors, and the actual errors in each case.

#NOTE:  This is just skeleton code, you will need to make sure to flesh it out with f(x),a,b,etc.

#For outputting the estimated errors, I wrote a function to find them using the bound given in equation (5.10).

def EstimatedTrapError(a,b,n,M):
    h=(b-a)/n
    error=(b-a)*M*h^2/12
    return error
    
#To determine the estimated errors, using (5.10) on pg. 185, I plotted f''(x) to determine a good choice for M.

#Then once you have a suitable M, I edited the books code for producing the trapezoidal estimate table to include the estimated error

print 'trapezoidal','\t\t','error','\t\t\t','estimatedtraperror'
print '----------------------------------------------------------------'
for power in [0,1,2,3,4,5]:
    n=2^power
    trapezoidal=CompositeTrapezoidal(f,0,1,n).n()
    error=(trapezoidal-integral(f(x),x,0,1)).n()
    estimatedtraperror=EstimatedTrapError(a,b,n,M).n()
    print trapezoidal,'\t',error,'\t',estimatedtraperror
    
#5.10.2.  Create a table of values obtained using the Composite Simpson's Rule, with n=2,4,8,16 and 32, for int(e^(-x^2)dx, 0,1).  Also output the estimated errors, and the actual errors in each case.

#For 5.10.2, try to write you own function for the Estimated Simpson's error, by editing the one I wrote for 5.10.1.  Similarly, edit the table to include your estimated error in the output.