166SWC, 7.7 - Example 7: Simpson's (Parabolic) Rule and Error Bounds Discussion (Answers)

Recall the Mathematica commands for approximation of
b
∫
a
f(x)x
with
M
n
or
T
n
:
In[]:=
MidpointSum[f_,a_,b_,n_]:=
n
∑
i=1
f[a+(2*i-1)*(b-a)/(2*n)]*(b-a)/n
In[]:=
Trapezoid[f_,a_,b_,n_]:=
n
∑
i=1
(1/2)(f[a+(i-1)*(b-a)/n]+f[a+i*(b-a)/n])*(b-a)/n

1. Construct a command 'Simpson[f_,a_,b_,n_]' that will perform Simpson's rule and compute
S
n
.
Note that n must be even and Δx=(b-a)/n.

In[]:=
Simpson1[f_,a_,b_,n_]:=1/3*Trapezoid[f,a,b,n/2]+2/3*MidpointSum[f,a,b,n/2]
In[]:=
Simpson2[f_,a_,b_,n_]:=(b-a)/(3*n)*Sum[f[a+(2*i-2)*(b-a)/n]+​​ 4*f[a+(2*i-1)*(b-a)/n]+f[a+2*i*(b-a)/n],{i,1,n/2}]
In[]:=
Simpson3[f_,a_,b_,n_]:=(b-a)/(3*n)*
n/2
∑
i=1
(f[a+(2*i-2)*(b-a)/n]+​​ 4*f[a+(2*i-1)*(b-a)/n]+f[a+2*i*(b-a)/n])

2.Testthiscommandontheintegral​​
1
∫
0
sin(x)xwithn=6andcomparetotheresultswefoundforapproximationsforthisintegralvia
𝐿
n
,
𝑅
n
,
𝑀
n
and
𝑇
n
inExample1.

Simpson1[Sin,0,1,6]
Out[]=
2
3
1
3
Sin
1
6
+
1
3
Sin
1
2
+
1
3
Sin
5
6
+
1
3
1
6
Sin
1
3
+
1
6
Sin
1
3
+Sin
2
3
+
1
6
Sin
2
3
+Sin[1]
N[%]
Out[]=
0.4597
Simpson2[Sin,0,1,6]
Out[]=
1
18
4Sin
1
6
+2Sin
1
3
+4Sin
1
2
+2Sin
2
3
+4Sin
5
6
+Sin[1]
N[%]
Out[]=
0.4597
Simpson3[Sin,0,1,6]
Out[]=
1
18
4Sin
1
6
+2Sin
1
3
+4Sin
1
2
+2Sin
2
3
+4Sin
5
6
+Sin[1]
N[%]
Out[]=
0.4597
Comparing the result of
S
6
=0.4597
to what we got in Example 1, namely
L
10
=0.417241
,
R
10
=0.501388
,
M
10
=0.459889
, and
T
10
=0.459315
, our estimate is closer to those for
M
10
and
T
10
, with about half the number of subintervals. Here’s the “actual” value of the integral:
NIntegrate[Sin[x],{x,0,1}]
Out[]=
0.459698

3. Find an upper bound for using Simpson’s Rule to estimate
1
∫
0
sin(
2
x
)x
with even positive integer n. Test this out with n = 10. Here are some ideas that may be useful:

To define the function
sin(
2
x
) with Mathematica, enter the command:
In[]:=
f[x_]:=Sin[x^2]
Here is how to find the fourth derivative of f:
f''''[x]
Out[]=
6(-4
2
x
Cos[
2
x
]-2Sin[
2
x
])+2x(-8xCos[
2
x
]-2x(2Cos[
2
x
]-4
2
x
Sin[
2
x
]))
If this is too complicated, you can use the 'Simplify' command:
Simplify[%]
Out[]=
-48
2
x
Cos[
2
x
]+4(-3+4
4
x
)Sin[
2
x
]
Here is one way to graph f''''[x]:
Plot[f''''[x],{x,0,1},Frame->True]
Out[]=
From this graph, we see that on [0, 1],

(4)
f
(x)<=30
. Taking
K
4
=30
in Equation (8), we have
I-
S
n
<=
5
K
4
(b-a)
180
4
n
. Thus with n = 10, we have an upper bound for
|I-
S
10
|
:
(30*(1)^5)/(180*10^4)
Out[]=
1
60000
N[%]
Out[]=
0.0000166667
Simpson1[f,0,1,10]
Out[]=
2
3
1
5
Sin
1
100
+
1
5
Sin
9
100
+
1
5
Sin
1
4
+
1
5
Sin
49
100
+
1
5
Sin
81
100
+
1
3
1
10
Sin
1
25
+
1
10
Sin
1
25
+Sin
4
25
+
1
10
Sin
4
25
+Sin
9
25
+
1
10
Sin
9
25
+Sin
16
25
+
1
10
Sin
16
25
+Sin[1]
N[%]
Out[]=
0.31026
Integrate[Sin[x^2],{x,0,1}]//N
Out[]=
0.310268
%-%%
Out[]=
8.06729×
-6
10

4.Findthesmallestevenpositiveintegernsothat​​I-
S
n
<0.001for
2
∫
-1
-
2
x
e
x.ComparethistowhatwefoundinExamples5and6.

In[]:=
g[x_]:=E^(-x^2)
g''''[x]
Out[]=
-6-2
-
2
x

+4
-
2
x

2
x
-2x8
-
2
x

x-2x-2
-
2
x

+4
-
2
x

2
x

Simplify[%]
Out[]=
4
-
2
x

(3-12
2
x
+4
4
x
)
Plot[g''''[x],{x,-1,2},Frame->True]
Out[]=
Solve[g'''''[x]==0,x]
Out[]=
{x0},x-
5
2
+
5
2
,x
5
2
+
5
2
,x-
1
2
(5-
10
)
,x
1
2
(5-
10
)

N[%]
Out[]=
{{x0.},{x-2.02018},{x2.02018},{x-0.958572},{x0.958572}}
Ng''''-1,-
1
2
5-
10

,0,
1
2
5-
10

,2
Out[]=
{-7.35759,-7.41948,12.,-7.41948,1.39199}
Thus,on[-1,2],|g''''(x)|≤12,sowecantake
K
4
=12inEquation(8)'serrorestimateforSimpson'sRule.Taken≥((
K
4
(b-a)^5)/(180*0.001))^(1/4)
((12*(3)^5)/(180*0.001))^(1/4)