1108.sagews

f = (1+x)^(1/2)
x0 = 0
xn = .1
def h(n):
    return (xn-x0)/n
    
def xset(n):
    return seq([x0 + i*h(n) for i in range(0,n+1,1)])
    
def trap1(f,h):
    var('x')
    var('y')
    approx = (h(1)/2)*(f(x0)+f(xn))
    actual = integral(f(x),x,x0,xn)
    error = n(actual - approx)
    fder = derivative(f(x),x,2)
    epsilon = solve(error == (h(1)^3/12)*fder(x),x)
    print 'approx: ', approx
    print 'error: ', error
    print 'mu: ', epsilon
    
trap1(f,h)

approx: 0.102440442408508 error: 1.93795829371e-05 mu: [ x == 1/158400*(-I*sqrt(3)*1702823^(2/3) - 800*198^(2/3) - 1702823^(2/3))*198^(1/3), x == 1/158400*(I*sqrt(3)*1702823^(2/3) - 800*198^(2/3) - 1702823^(2/3))*198^(1/3), sqrt(x + 1) == 1/3960*(-1)^(1/3)*198^(2/3)*1702823^(1/3) ]

x1=seq([0+i*.5 for i in range(0,10,1)])
x1[2]

1.00000000000000

def simp1(f,h):
    var('x')
    m = 2
    xi = xset(m)
    approx = (h(m)/3)*(f(xi[0])+4*f(xi[1])+f(xi[2]))
    actual = integral(f(x),x,x0,xn)
    error = n(actual - approx)
    fder = derivative(f(x),x,4)
    epsilon = solve(error == (h(m)^5/90)*fder(x),x)
    print 'approx: ', approx
    print 'error: ', error
    print 'mu: ', epsilon
    
simp1(f,h)

approx: 0.102459819242567 error: 2.74887818064e-09 mu: [ x == 1/307199999*(363784763^(2/7)*e^(4/7*I*pi) - 307199999^(2/7))*307199999^(5/7), x == -1/307199999*(307199999^(2/7)*e^(6/7*I*pi) - 363784763^(2/7))*307199999^(5/7)*e^(-6/7*I*pi), x == -1/307199999*(307199999^(2/7)*e^(2/7*I*pi) - 363784763^(2/7))*307199999^(5/7)*e^(-2/7*I*pi), x == 1/307199999*(363784763^(2/7)*e^(2/7*I*pi) - 307199999^(2/7))*307199999^(5/7), x == 1/307199999*(363784763^(2/7)*e^(6/7*I*pi) - 307199999^(2/7))*307199999^(5/7), x == -1/307199999*(307199999^(2/7)*e^(4/7*I*pi) - 363784763^(2/7))*307199999^(5/7)*e^(-4/7*I*pi), sqrt(x + 1) == 1/307199999*(-1)^(1/7)*307199999^(6/7)*363784763^(1/7) ]

def simp2(f,h):
    var('x')
    m = 3
    xi = xset(m)
    approx = (3*h(m)/8)*(f(xi[0])+3*f(xi[1])+3*f(xi[2])+f(xi[3]))
    actual = integral(f(x),x,x0,xn)
    error = n(actual - approx)
    fder = derivative(f(x),x,4)
    epsilon = solve(error == (3*h(m)^5/80)*fder(x),x)
    print 'approx: ', approx
    print 'error: ', error
    print 'mu: ', epsilon
    
simp2(f,h)

approx: 0.102459820769143 error: 1.22230135735e-09 mu: [ x == 1/10800000*(409064423^(2/7)*e^(4/7*I*pi) - 4*2700000^(2/7))*2700000^(5/7), x == -1/10800000*(4*2700000^(2/7)*e^(6/7*I*pi) - 409064423^(2/7))*2700000^(5/7)*e^(-6/7*I*pi), x == -1/10800000*(4*2700000^(2/7)*e^(2/7*I*pi) - 409064423^(2/7))*2700000^(5/7)*e^(-2/7*I*pi), x == 1/10800000*(409064423^(2/7)*e^(2/7*I*pi) - 4*2700000^(2/7))*2700000^(5/7), x == 1/10800000*(409064423^(2/7)*e^(6/7*I*pi) - 4*2700000^(2/7))*2700000^(5/7), x == -1/10800000*(4*2700000^(2/7)*e^(4/7*I*pi) - 409064423^(2/7))*2700000^(5/7)*e^(-4/7*I*pi), sqrt(x + 1) == 1/5400000*(-1)^(1/7)*2700000^(6/7)*409064423^(1/7) ]

def simp3(f,h):
    var('x')
    m = 4
    xi = xset(m)
    approx = (2*h(m)/45)*(7*f(xi[0])+32*f(xi[1])+12*f(xi[2])+32*f(xi[3])+7*f(xi[4]))
    actual = integral(f(x),x,x0,xn)
    error = n(actual - approx)
    fder = derivative(f(x),x,6)
    epsilon = solve(error == (8*h(m)^7/945)*fder(x),x)
    print 'approx: ', approx
    print 'error: ', error
    print 'mu: ', epsilon
    
simp3(f,h)

approx: 0.102459821990859 error: 5.85295700795e-13 mu: [ x == 1/1310719999999*(1708538092185^(2/11)*e^(4/11*I*pi) - 1310719999999^(2/11))*1310719999999^(9/11), x == 1/1310719999999*(1708538092185^(2/11)*e^(8/11*I*pi) - 1310719999999^(2/11))*1310719999999^(9/11), x == -1/1310719999999*(1310719999999^(2/11)*e^(10/11*I*pi) - 1708538092185^(2/11))*1310719999999^(9/11)*e^(-10/11*I*pi), x == -1/1310719999999*(1310719999999^(2/11)*e^(6/11*I*pi) - 1708538092185^(2/11))*1310719999999^(9/11)*e^(-6/11*I*pi), x == -1/1310719999999*(1310719999999^(2/11)*e^(2/11*I*pi) - 1708538092185^(2/11))*1310719999999^(9/11)*e^(-2/11*I*pi), x == 1/1310719999999*(1708538092185^(2/11)*e^(2/11*I*pi) - 1310719999999^(2/11))*1310719999999^(9/11), x == 1/1310719999999*(1708538092185^(2/11)*e^(6/11*I*pi) - 1310719999999^(2/11))*1310719999999^(9/11), x == 1/1310719999999*(1708538092185^(2/11)*e^(10/11*I*pi) - 1310719999999^(2/11))*1310719999999^(9/11), x == -1/1310719999999*(1310719999999^(2/11)*e^(8/11*I*pi) - 1708538092185^(2/11))*1310719999999^(9/11)*e^(-8/11*I*pi), x == -1/1310719999999*(1310719999999^(2/11)*e^(4/11*I*pi) - 1708538092185^(2/11))*1310719999999^(9/11)*e^(-4/11*I*pi), sqrt(x + 1) == 1/1310719999999*(-1)^(1/11)*1310719999999^(10/11)*1708538092185^(1/11) ]

def ech(n):
    return ((xn-x0)/(n+2))
    
def exset(n):
    return seq([x0 + i*ech(n) for i in range(0,n+2,1)])

def mid1(f,h):
    var('x')
    m = 0
    xi = exset(m)
    approx = (2*ech(m))*f(xi[0])
    actual = integral(f(x),x,x0,xn)
    error = n(actual - approx)
    print 'approx: ', approx
    print 'error: ', error
    
mid1(f,h)

approx: 0.100000000000000 error: 0.00245982199144

def mid2(f,h):
    var('x')
    m = 1
    xi = exset(m)
    approx = (3*ech(m)/2)*(f(xi[0])+f(xi[1]))
    actual = integral(f(x),x,x0,xn)
    error = n(actual - approx)
    print 'approx: ', approx
    print 'error: ', error
    
mid2(f,h)

approx: 0.100826502273256 error: 0.00163331971819

def mid3(f,h):
    var('x')
    m = 2
    xi = exset(m)
    approx = (4*ech(m)/3)*(2*f(xi[0])-f(xi[1])+2*f(xi[2]))
    actual = integral(f(x),x,x0,xn)
    error = n(actual - approx)
    print 'approx: ', approx
    print 'error: ', error
    
mid3(f,h)

approx: 0.101232243887536 error: 0.00122757810391

def mid4(f,h):
    var('x')
    m = 3
    xi = exset(m)
    approx = (5*ech(m)/24)*(11*f(xi[0])+f(xi[1])+f(xi[2])+11*f(xi[3]))
    actual = integral(f(x),x,x0,xn)
    error = n(actual - approx)
    print 'approx: ', approx
    print 'error: ', error
    
mid4(f,h)

approx: 0.101478948131835 error: 0.000980873859609 0.000980873859609

%latex
Problem 4.3.21a
\begin{equation}
f(x) = \left(1+x\right)^{1/2}
\end{equation}
Are the accuracies of the approximations consistent with the error formulas?\\
n=1, Trapezoidal Rule\\
approx:  0.102440442408508\\
error:  1.93795829371e-05\\
\\
n=2, Simpson's Rule\\
approx:  0.102459819242567\\
error:  2.74887818064e-09\\
\\
n=3, Simpson's Three-eigths rule\\
approx:  0.102459820769143\\
error:  1.22230135735e-09\\
\\
n=4, Simpson's Rule\\
approx:  0.102459821990859\\
error:  5.85295700795e-13\\
\\
n=0, Midpoint Rule\\
approx:  0.100000000000000\\
error:  0.00245982199144\\
\\
n=1, Midpoint Rule\\
approx:  0.100826502273256\\
error:  0.00163331971819\\
\\
n=2, Midpoint Rule\\
approx:  0.101232243887536\\
error:  0.00122757810391\\
\\
n=3, Midpoint Rule\\
approx:  0.101478948131835\\
error:  0.000980873859609