Views : 542

Do this first!

Please make a copy of this file called Assignment_5_turnin.sagews. That file (or a corresponding Jupyter notebook with extension .ipynb) will be graded by the grader.

Assignment 5

Please do the following problems.

5.9: Theory Exercises

These problems may be written by hand and turned in, or they may be typeset in this Jupyter notebook or in a Sage worksheet.

1. Problem 5.9.3

n=3
f(x)=e^(-x^2)
a=0
b=1
h=(b-a)/n
Trap=h/2*(f(a)+2*(f(a+h)+f(a+2*h))+f(b))
print 'The answer for this integral calculated with the Composite Trapezoidal Rule is', Trap.n(digits=5)
fp(x)=diff(f(x),x)
fpp(x)=diff(fp(x),x)
p1=plot(f(x),(x,0,1), color="green",legend_label="F(x)")
p2=plot(fp(x),(x,0,1), color="blue",legend_label="First Derivative")
p3=plot(fpp(x),(x,0,1), color="red",legend_label="Second Derivative")
show(p1+p2+p3)
#max of fpp(x) is at x=0 where y=-2

error=(b-a)*h^2/12*fpp(0)
print 'The error of approximation is', error.n(digits=5)

The answer for this integral calculated with the Composite Trapezoidal Rule is 0.73999

The error of approximation is -0.018519

2. Problem 5.9.4

n=2
f(x)=e^(-x^2)
a=0
b=1
h=(b-a)/n
Simpson=h/3*(f(a)+4*f(a+h)+f(b))
print 'The answer for this integral calculated with the Composite Simpson Rule is', Simpson.n()
fp(x)=diff(f(x),x)
fpp(x)=diff(fp(x),x)
fppp(x)=diff(fpp(x),x)
fpppp(x)=diff(fppp(x),x)
fppppp(x)=diff(fpppp(x),x)
p1=plot(f(x),(x,0,1), color="green",legend_label="F(x)")
p2=plot(fp(x),(x,0,1), color="blue",legend_label="First Derivative")
p3=plot(fpp(x),(x,0,1), color="red",legend_label="Second Derivative")
p4=plot(fpp(x),(x,0,1), color="orange",legend_label="Third Derivative")
p5=plot(fppp(x),(x,0,1), color="purple",legend_label="Fourth Derivative")
p6=plot(fpppp(x),(x,0,1), color="yellow",legend_label="Fifth Derivative")
show(p1+p2+p3+p4+p5+p6)
#max of fpppp(x) is at x=0.525 where y=4

error=(b-a)*h^4/180*fpppp(0.525)
print 'The error of approximation is', error.n(digits=4)

The answer for this integral calculated with the Composite Simpson Rule is 0.747180428909510

The error of approximation is -3.820e-6

3. Problem 5.9.16

n=4
f(x)=arctan(x)
a=0
b=1
h=(b-a)/n
Trap=h/2*(f(a)+2*(f(a+h)+f(a+2*h)+f(a+3*h))+f(b))
print 'The answer for this integral calculated with the Composite Trapezoidal Rule is', Trap.n(digits=5)

fp(x)=diff(f(x),x)
fpp(x)=diff(fp(x),x)
fppp(x)=diff(fpp(x),x)
p1=plot(f(x),(x,0,1), color="green",legend_label="F(x)")
p2=plot(fp(x),(x,0,1), color="blue",legend_label="First Derivative")
p3=plot(fpp(x),(x,0,1), color="red",legend_label="Second Derivative")
p4=plot(fppp(x),(x,0,1), color="yellow",legend_label="Third Derivative")
show(p1+p2+p3+p4)
#max of fpp(x) is at x=0.575 where y=-0.7

error=(b-a)*h^2/12*fpp(0.575)
print 'The error of approximation is', error.n(digits=5)

The answer for this integral calculated with the Composite Trapezoidal Rule is 0.43621

The error of approximation is -0.0033829

4. Problem 5.9.17

n=4
f(x)=arctan(x)
a=0
b=1
h=(b-a)/n
Simpson=h/3*(f(a)+4*f(a+h)+2*f(a+2*h)+4*f(a+3*h)+f(b))
print 'The answer for this integral calculated with the Composite Simpson Rule is', Simpson.n(digits=5)
fp(x)=diff(f(x),x)
fpp(x)=diff(fp(x),x)
fppp(x)=diff(fpp(x),x)
fpppp(x)=diff(fppp(x),x)
fppppp(x)=diff(fpppp(x),x)
p1=plot(f(x),(x,0,1), color="green",legend_label="F(x)")
p2=plot(fp(x),(x,0,1), color="blue",legend_label="First Derivative")
p3=plot(fpp(x),(x,0,1), color="red",legend_label="Second Derivative")
p4=plot(fpp(x),(x,0,1), color="orange",legend_label="Third Derivative")
p5=plot(fppp(x),(x,0,1), color="purple",legend_label="Fourth Derivative")
p6=plot(fpppp(x),(x,0,1), color="yellow",legend_label="Fifth Derivative")
show(p1+p2+p3+p4+p5+p6)
#max of fpppp(x) is at x=0 where y=-2

error=(b-a)*h^4/180*fpppp(0)
print 'The error of approximation is', error.n(digits=8)

The answer for this integral calculated with the Composite Simpson Rule is 0.43888

The error of approximation is 0.00000000

5. Problem 5.9.18

x0=0.2
x1=0.4
x2=0.6
x3=0.8
x4=1.0
y0=7.11
y1=10.21
y2=20.98
y3=32.91
y4=1.12
h=(1-.2)/4
trap=(h/2)*(y0+y4)+h*(y1+y2+y3)
print 'The estimate of the integral using the Composite Trapezoidal Method is', trap.n(digits=5)

The estimate of the integral using the Composite Trapezoidal Method is 13.643

6. Problem 5.9.19

x0=0.2
x1=0.4
x2=0.6
x3=0.8
x4=1.0
y0=7.11
y1=10.21
y2=20.98
y3=32.91
y4=1.12
h=(1-.2)/4
simpson=h/3*(y0+4*y1+2*y2+4*y3+y4)
print 'The estimate of the integral using the Composite Simpson Method is', simpson.n(digits=5)

The estimate of the integral using the Composite Simpson Method is 14.845

7. Problem 5.9.21

f(x)=e^(-x^2)
a=-1
b=1
h=(b-a)/n
error=0.0001
fp(x)=diff(f(x),x)
fpp(x)=diff(fp(x),x)
fppp(x)=diff(fpp(x),x)
p1=plot(f(x),(x,-1,1), color="green",legend_label="F(x)")
p2=plot(fp(x),(x,-1,1), color="blue",legend_label="First Derivative")
p3=plot(fpp(x),(x,-1,1), color="red",legend_label="Second Derivative")
p4=plot(fppp(x),(x,-1,1), color="orange",legend_label="Third Derivative")
show(p1+p2+p3+p4)
abs(fpp(0))
#0.0001=(1/6)*2*(2/n)^2
#0.0001*6/2=(2/n)^2
#sqrt(0.0001*3)=2/n
n=2/sqrt(0.0001*3)
print n
print 'Round up to 116 nodes.'

2
115.470053837925
Round up to 116 nodes.

5.10: Computer Exercises.

Please use Sage to do the following problem. Do not use the built-in Python functions for the trapezoidal or Simpson's rule.

8. Problem 5.10.1

def CompositeTrap (f,a,b,n):
    h=(b-a)/n
    sum=1/2*f(a)
    for i in [1,2,..,n-1]:
        sum+=f(a+i*h)
    sum+=1/2*f(b)
    sum=sum*h
    return sum

f(x)=exp(-x^2)
fp(x)=diff(f(x),x)
fpp(x)=diff(fp(x),x)
a=0
b=1
print 'n', '\t','trapezoidal', '\t\t', 'actual error', '\t\t', 'estimated error'
print '------------------------------------------------------------------------'
for power in [0,1,2,3,4,5]:
    n=2^power
    h=(b-a)/n
    trap=CompositeTrap(f,a,b,n).n()
    actualerror=(trap-numerical_integral(f,a,b)[0]).n()
    estimatederror=(b-a)*h^2/12*fpp(0).n(digits=5)
    print n, '\t', trap,'\t',actualerror, '\t', estimatederror

n 	trapezoidal 		actual error 		estimated error
------------------------------------------------------------------------
0.683939720585721 	-0.0628844122267058 	-0.16667
0.731370251828563 	-0.0154538809838640 	-0.041667
0.742984097800381 	-0.00384003501204577 	-0.010417
0.745865614845695 	-0.000958517966731742 	-0.0026042
0.746584596788222 	-0.000239536024205456 	-0.00065104
0.746764254652294 	-0.0000598781601327536 	-0.00016276

9. Problem 5.10.2

def CompositeSimpson(f,a,b,n):
    if int(n)%2==1:
        return "ERROR-n should be even!!!"
    h=(b-a)/n
    sum=f(a)
    for i in [1,2..,n-1]:
        if(i%2==1):
            sum+=4*f(a+i*h)
        else:
            sum+=2*f(a+i*h)
    sum+=f(b)
    sum=sum*h/3
    return sum
f(x)=exp(-x^2)
fp(x)=diff(f(x),x)
fpp(x)=diff(fp(x),x)
fppp(x)=diff(fpp(x),x)
fpppp(x)=diff(fppp(x),x)
a=0
b=1
print 'n', '\t','simpsons', '\t\t', 'actual error', '\t\t', 'estimated error'
print '------------------------------------------------------------------------'
for power in [1,2,3,4,5]:
    n=2^power
    h=(b-a)/n
    simpsons=CompositeSimpson(f,a,b,n).n()
    actualerror=(simpsons-numerical_integral(f,a,b)[0]).n()
    estimatederror=(b-a)*h^4/180*fpppp(0.525).n(digits=5)
    print n, '\t', trap,'\t',actualerror, '\t', estimatederror

n 	simpsons 		actual error 		estimated error
------------------------------------------------------------------------
0.746764254652294 	0.000356296097083320 	-3.8202e-6
0.746764254652294 	0.0000312469785601621 	-2.3876e-7
0.746764254652294 	1.98771503967521e-6 	-1.4923e-8
0.746764254652294 	1.24623303232596e-7 	-9.3267e-10
0.746764254652294 	7.79455788801897e-9 	-5.8292e-11

10. Problem 5.10.5

f(x)=exp(cos(pi*x))
fp(x)=diff(f(x),x)
fpp(x)=diff(fp(x),x)
a=-1
b=1
p1=plot(f(x),(x,-1,1), color="green",legend_label="F(x)")
p2=plot(fp(x),(x,-1,1), color="blue",legend_label="First Derivative")
p3=plot(fpp(x),(x,-1,1), color="red",legend_label="Second Derivative")
show(p1+p2+p3)
print 'n', '\t','trapezoidal', '\t\t', 'actual error', '\t\t', 'estimated error'
print '------------------------------------------------------------------------'
for power in [0,1,2,3,4,5]:
    n=2^power
    h=(b-a)/n
    trap=CompositeTrap(f,a,b,n).n()
    actualerror=(trap-numerical_integral(f,a,b)[0]).n()
    estimatederror=(b-a)*h^2/12*fpp(0).n(digits=5)
    print n, '\t', trap,'\t',actualerror, '\t', estimatederror

n 	trapezoidal 		actual error 		estimated error
------------------------------------------------------------------------
0.735758882342885 	-1.79637287316113 	-17.886
3.08616126963049 	0.554029514126471 	-4.4714
2.54308063481524 	0.0109488793112269 	-1.1178
2.53213215392898 	3.98424961023380e-7 	-0.27946
2.53213175550402 	-4.44089209850063e-16 	-0.069866
2.53213175550402 	0.000000000000000 	-0.017466

11. Problem 5.10.6

f(x)=exp(cos(pi*x))
fp(x)=diff(f(x),x)
fpp(x)=diff(fp(x),x)
fppp(x)=diff(fpp(x),x)
fpppp(x)=diff(fppp(x),x)
a=-1
b=1
print 'n', '\t','simpsons', '\t\t', 'actual error', '\t\t', 'estimated error'
print '------------------------------------------------------------------------'
for power in [1,2,3,4,5]:
    n=2^power
    h=(b-a)/n
    simpsons=CompositeSimpson(f,a,b,n).n()
    actualerror=(simpsons-numerical_integral(f,a,b)[0]).n()
    estimatederror=(b-a)*h^4/180*fpppp(0).n(digits=5)
    print n, '\t', trap,'\t',actualerror, '\t', estimatederror

n 	simpsons 		actual error 		estimated error
------------------------------------------------------------------------
2.53213175550402 	1.33749697655567 	11.768
2.53213175550402 	-0.170077998960521 	0.73552
2.53213175550402 	-0.00364909520379397 	0.045970
2.53213175550402 	-1.32808321229305e-7 	0.0028731
2.53213175550402 	-4.44089209850063e-16 	0.00017957