# Math 210 Sec. 201
# Solutions to Assignment 3
#
# 1.
> r:= cos(theta)^2:
> L:= Int(sqrt(r^2 + diff(r,theta)^2), theta = -Pi/2 .. Pi/2);
1/2 Pi
/
| 4 2 2 1/2
L := | (cos(theta) + 4 cos(theta) sin(theta) ) dtheta
|
/
- 1/2 Pi
> L:= int(cos(theta)*sqrt(cos(theta)^2+4*sin(theta)^2), theta = -Pi/2 .. Pi/2);
1/2 Pi
/
| 2 2 1/2
L := | cos(theta) (cos(theta) + 4 sin(theta) ) dtheta
|
/
- 1/2 Pi
> with(student):
> changevar(x=sin(theta),L,x);
1/2 1/2
2 + 1/3 3 arcsinh(3 )
> convert(",ln);
1/2 1/2
2 + 1/3 3 ln(2 + 3 )
> evalf(");
2.760345997
# 2. (a)
> f:= x -> x^11:
# The following definitions are copied from Lesson 18.
> h:= n ->(b-a)/n:
> X:= (i,n) -> a + i*(b-a)/n:
> midrule:= n ->
> h(n)*sum(f((X(i,n)+X(i+1,n))/2), i = 0 .. n-1);
> traprule:= n ->
> h(n) * sum((f(X(i,n))+f(X(i+1,n)))/2, i=0..n-1);
> simpson:= n ->
> h(n) * sum(1/3*f(X(2*i,n))+4/3*f(X(2*i+1,n))+1/3*f(X(2*i+2,n)), i=0..n/2-1);
> J:= int(f(x), x=a..b);
/n - 1 \
|----- |
| \ |
midrule := n -> h(n) | ) f(1/2 X(i, n) + 1/2 X(i + 1, n))|
| / |
|----- |
\i = 0 /
/n - 1 \
|----- |
| \ |
traprule := n -> h(n) | ) (1/2 f(X(i, n)) + 1/2 f(X(i + 1, n)))|
| / |
|----- |
\i = 0 /
simpson := n -> h(n)
/1/2 n - 1 \
| ----- |
| \ |
| ) (1/3 f(X(2 i, n)) + 4/3 f(X(2 i + 1, n)) + 1/3 f(X(2 i + 2, n)))|
| / |
| ----- |
\ i = 0 /
b
/
|
J := | f(x) dx
|
/
a
# Now define ME, TE and SE as functions of n.
> ME:= unapply(J - midrule(n), n);
10 8 3
12 12 a b a b
ME := n -> 1/12 b - 1/12 a - (b - a) (7665/2048 ----- + 1705/96 -----
9 5
n n
3 8 8 3 10
a b 9 2 a b a b
- 48895/1024 ----- + 1/12 b a n + 63875/2048 ----- + 7665/2048 -----
7 9 9
n n n
11 4 7 9 2
a a b 7 4 a b
+ 77/64 --- - 38325/1024 ----- + 1/12 a b n - 1705/96 -----
3 9 5
n n n
5 6 7 4 5 6
a b a b 6 5 a b
+ 341/192 ----- - 1705/192 ----- + 1/12 a b n + 17885/1024 -----
5 5 9
n n n
3 8 10 7 4 2 9
a b a b a b a b
+ 63875/2048 ----- + 1705/192 ----- + 23749/512 ----- + 29337/1024 -----
9 5 7 7
n n n n
8 3 6 5
a b 10 10 a b
- 77/64 ----- + 1/12 a b n + 1/12 a b n + 341/192 -----
3 5
n n
7 4 2 9 10 4 7
a b a b a b a b
- 38325/1024 ----- - 89425/6144 ----- - 231/64 ----- + 23749/512 -----
9 9 3 7
n n n n
8 3 9 2 10
9 2 a b b a a b
+ 1/12 a b n - 48895/1024 ----- - 1705/96 ----- + 11/24 -----
7 5 n
n n
10 10 6 5 2 9
a b a b a b a b
+ 11/24 ----- + 1705/192 ----- - 9779/512 ----- + 231/64 -----
n 5 7 3
n n n
6 5 8 3 3 8
a b b a a b 3 8
+ 17885/1024 ----- + 1705/96 ----- - 77/64 ----- + 1/12 a b n
9 5 3
n n n
4 7 10 9 2
a b b a 4 7 a b
- 1705/192 ----- - 9779/1024 ----- + 1/12 a b n - 89425/6144 -----
5 7 9
n n n
10 9 2 5 6
b a a b a b 8 3
- 231/64 ----- + 29337/1024 ----- - 9779/512 ----- + 1/12 a b n
3 7 7
n n n
9 2 10 11
a b 5 6 a b a 11
+ 231/64 ----- + 1/12 a b n - 9779/1024 ----- - 11/24 --- + 1/12 n b
3 7 n
n n
11 11 11 11
a a a b
- 341/192 --- + 1397/1024 --- - 2555/6144 --- - 2555/6144 ---
5 7 9 9
n n n n
11 11 11 11
b b b b 11
+ 1397/1024 --- - 341/192 --- + 77/64 --- - 11/24 --- + 1/12 a n)/n
7 5 3 n
n n n
> TE:= unapply(J - traprule(n), n);
10 8 3
12 12 a b a b
TE := n -> 1/12 b - 1/12 a - (b - a) (- 15/4 ----- - 55/3 -----
9 5
n n
3 8 8 3 10 11
a b 9 2 a b a b a
+ 385/8 ----- + 1/12 b a n - 125/4 ----- - 15/4 ----- - 11/8 ---
7 9 9 3
n n n n
4 7 9 2 5 6 7 4
a b 7 4 a b a b a b
+ 75/2 ----- + 1/12 a b n + 55/3 ----- - 11/6 ----- + 55/6 -----
9 5 5 5
n n n n
5 6 3 8 10 7 4
6 5 a b a b a b a b
+ 1/12 a b n - 35/2 ----- - 125/4 ----- - 55/6 ----- - 187/4 -----
9 9 5 7
n n n n
2 9 8 3 6 5
a b a b 10 10 a b
- 231/8 ----- + 11/8 ----- + 1/12 a b n + 1/12 a b n - 11/6 -----
7 3 5
n n n
7 4 2 9 10 4 7
a b a b a b a b 9 2
+ 75/2 ----- + 175/12 ----- + 33/8 ----- - 187/4 ----- + 1/12 a b n
9 9 3 7
n n n n
8 3 9 2 10 10 10
a b b a a b a b a b
+ 385/8 ----- + 55/3 ----- - 11/12 ----- - 11/12 ----- - 55/6 -----
7 5 n n 5
n n n
6 5 2 9 6 5 8 3 3 8
a b a b a b b a a b
+ 77/4 ----- - 33/8 ----- - 35/2 ----- - 55/3 ----- + 11/8 -----
7 3 9 5 3
n n n n n
4 7 10 9 2
3 8 a b b a 4 7 a b
+ 1/12 a b n + 55/6 ----- + 77/8 ----- + 1/12 a b n + 175/12 -----
5 7 9
n n n
10 9 2 5 6 9 2
b a a b a b 8 3 a b
+ 33/8 ----- - 231/8 ----- + 77/4 ----- + 1/12 a b n - 33/8 -----
3 7 7 3
n n n n
10 11 11 11
5 6 a b a 11 a a
+ 1/12 a b n + 77/8 ----- + 11/12 --- + 1/12 n b + 11/6 --- - 11/8 ---
7 n 5 7
n n n
11 11 11 11 11 11
a b b b b b
+ 5/12 --- + 5/12 --- - 11/8 --- + 11/6 --- - 11/8 --- + 11/12 ---
9 9 7 5 3 n
n n n n n
11
+ 1/12 a n)/n
> SE:= unapply(J - simpson(n), n);
10 8 3
12 12 a b a b
SE := n -> 1/12 b - 1/12 a - (b - a) (1275 ----- + 1100/3 -----
9 5
n n
3 8 8 3 10 11
a b 9 2 a b a b a
- 8085/2 ----- + 1/12 b a n + 10625 ----- + 1275 ----- + 11/2 ---
7 9 9 3
n n n n
4 7 9 2 5 6 7 4
a b 7 4 a b a b a b
- 12750 ----- + 1/12 a b n - 1100/3 ----- + 110/3 ----- - 550/3 -----
9 5 5 5
n n n n
5 6 3 8 10 7 4
6 5 a b a b a b a b
+ 1/12 a b n + 5950 ----- + 10625 ----- + 550/3 ----- + 3927 -----
9 9 5 7
n n n n
2 9 8 3 6 5
a b a b 10 10 a b
+ 4851/2 ----- - 11/2 ----- + 1/12 a b n + 1/12 a b n + 110/3 -----
7 3 5
n n n
7 4 2 9 10 4 7
a b a b a b a b 9 2
- 12750 ----- - 14875/3 ----- - 33/2 ----- + 3927 ----- + 1/12 a b n
9 9 3 7
n n n n
8 3 9 2 10 6 5 2 9
a b b a a b a b a b
- 8085/2 ----- - 1100/3 ----- + 550/3 ----- - 1617 ----- + 33/2 -----
7 5 5 7 3
n n n n n
6 5 8 3 3 8 4 7
a b b a a b 3 8 a b
+ 5950 ----- + 1100/3 ----- - 11/2 ----- + 1/12 a b n - 550/3 -----
9 5 3 5
n n n n
10 9 2 10 9 2
b a 4 7 a b b a a b
- 1617/2 ----- + 1/12 a b n - 14875/3 ----- - 33/2 ----- + 4851/2 -----
7 9 3 7
n n n n
5 6 9 2 10
a b 8 3 a b 5 6 a b
- 1617 ----- + 1/12 a b n + 33/2 ----- + 1/12 a b n - 1617/2 -----
7 3 7
n n n
11 11 11 11 11
11 a a a b b
+ 1/12 n b - 110/3 --- + 231/2 --- - 425/3 --- - 425/3 --- + 231/2 ---
5 7 9 9 7
n n n n n
11 11
b b 11
- 110/3 --- + 11/2 --- + 1/12 a n)/n
5 3
n n
# b) To plot these, we need to take particular values of a and b,
# e.g. a=0, b=1.
> a:= 0; b:= 1;
> plot({n^2 * ME(n), n^2* TE(n)}, n = 1 .. 20);
a := 0
b := 1
# The Midpoint error is on the top, and the Trapezoid error is on the bottom. We see that the Midpoint
# error has about half the absolute value of the Trapezoid error, and the opposite sign.
# The Midpoint error is approximately 0.45/n^2, and the Trapezoid error approximately -0.9/n^2.
> plot(n^4*SE(n), n = 2 .. 40);
# The Simpson's Rule error is approximately -5.4/n^4 for large n.
#
# (c) For this we need a and b to be symbolic again.
> a:= 'a': b:= 'b':
> Mcoeffs:= [seq(coeff(expand(ME(n)),n,-2*jj), jj=1..5)];
Mcoeffs := [
11 11 11 10 2 11 10 2 11 12 11 11 11 12
---- b a + ---- b a - ---- a b + ---- b - ---- b a - ---- a ,
12 24 24 24 12 24
77 11 77 8 4 77 9 3 77 11 231 10 2 77 9 3
- ---- b a + ---- a b - ---- a b + ---- b a - --- b a + ---- b a
16 64 16 16 32 16
231 10 2 77 12 77 12 77 8 4
+ --- a b + ---- a - ---- b - ---- b a ,
32 64 64 64
1705 8 4 1705 10 2 1705 9 3 341 7 5 1705 8 4 1705 9 3
---- b a + ---- b a - ---- b a + --- a b - ---- a b + ---- a b
64 64 48 32 64 48
1705 10 2 341 11 341 5 7 341 12 341 12 341 11
- ---- a b - --- b a - --- a b + --- b - --- a + --- b a ,
64 32 32 192 192 32
1397 11 96393 8 4 9779 9 3 1397 11 9779 9 3 96393 8 4
- ---- b a + ----- a b - ---- a b + ---- b a + ---- b a - ----- b a
128 1024 128 128 128 1024
9779 10 2 9779 10 2 4191 5 7 4191 7 5 1397 12 1397 12
+ ---- a b - ---- b a + ---- a b - ---- a b - ---- b + ---- a
256 256 64 64 1024 1024
,
12775 11 12775 11 140525 8 4 28105 5 7 140525 9 3
----- b a - ----- b a + ------ b a - ----- a b - ------ b a
3072 3072 2048 512 3072
28105 10 2 28105 7 5 140525 8 4 140525 9 3 28105 10 2
+ ----- b a + ----- a b - ------ a b + ------ a b - ----- a b
1536 512 2048 3072 1536
2555 12 2555 12
+ ---- b - ---- a ]
6144 6144
> shouldbe := [ seq((b-a)^(2*jj)*((D@@(2*jj-1))(f)(b) - (D@@(2*jj-1))(f)(a)), jj = 1 .. 5) ];
2 10 10 4 8 8
shouldbe := [(b - a) (11 b - 11 a ), (b - a) (990 b - 990 a ),
6 6 6 8 4 4
(b - a) (55440 b - 55440 a ), (b - a) (1663200 b - 1663200 a ),
10 2 2
(b - a) (19958400 b - 19958400 a )]
> Mcs:= seq(normal(Mcoeffs[kk]/shouldbe[kk]), kk = 1 .. 5);
31 127 73
Mcs := 1/24, -7/5760, ------, - ---------, ----------
967680 154828800 3503554560
> Tcoeffs:= [seq(coeff(expand(TE(n)),n,-2*jj), jj=1..5)];
Tcoeffs := [
11 11 12 11 10 2 11 12 11 10 2 11
11/6 b a - ---- b + ---- a b + ---- a - ---- b a - 11/6 b a ,
12 12 12 12
12 9 3 10 2 8 4 9 3 11
11/8 b + 11/2 a b - 33/4 a b + 11/8 b a - 11/2 b a - 11/2 b a
10 2 8 4 11 12
+ 33/4 b a - 11/8 a b + 11/2 b a - 11/8 a ,
9 3 10 2 8 4 12 8 4 9 3
- 110/3 a b + 55/2 a b + 55/2 a b + 11/6 a - 55/2 b a + 110/3 b a
7 5 10 2 12 11 5 7 11
- 11 a b - 55/2 b a - 11/6 b - 11 b a + 11 a b + 11 b a,
12 9 3 5 7 12 9 3 10 2 11
11/8 b + 77 a b - 66 a b - 11/8 a - 77 b a - 77/2 a b - 11 b a
8 4 8 4 11 7 5 10 2
+ 759/8 b a - 759/8 a b + 11 b a + 66 a b + 77/2 b a ,
8 4 9 3 11 12 12 9 3
- 275/4 b a + 275/6 b a + 25/6 b a - 5/12 b + 5/12 a - 275/6 a b
7 5 10 2 8 4 10 2 11
- 55 a b - 55/3 b a + 275/4 a b + 55/3 a b - 25/6 b a
5 7
+ 55 a b ]
> Tcs:= seq(normal(Tcoeffs[kk]/shouldbe[kk]), kk = 1 .. 5);
Tcs := -1/12, 1/720, -1/30240, 1/1209600, -1/47900160
> Scoeffs:= [seq(coeff(expand(SE(n)),n,-2*jj), jj=1..5)];
Scoeffs := [0,
8 4 11 10 2 8 4 9 3 9 3
- 11/2 b a - 22 b a + 33 a b + 11/2 a b - 22 a b + 22 b a
11 10 2 12 12
+ 22 b a - 33 b a - 11/2 b + 11/2 a ,
10 2 11 9 3 5 7 11 8 4
550 b a - 220 b a + 2200/3 a b - 220 a b + 220 b a + 550 b a
10 2 8 4 7 5 9 3 12 12
- 550 a b - 550 a b + 220 a b - 2200/3 b a + 110/3 b - 110/3 a
,
11 10 2 9 3 8 4 10 2 11
- 924 b a + 3234 a b - 6468 a b - 15939/2 b a - 3234 b a + 924 b a
8 4 7 5 9 3 5 7 12
+ 15939/2 a b - 5544 a b + 6468 b a + 5544 a b - 231/2 b
12
+ 231/2 a ,
8 4 10 2 9 3 5 7 11
23375 b a + 18700/3 b a + 46750/3 a b - 18700 a b - 4250/3 b a
10 2 11 9 3 7 5 8 4
- 18700/3 a b + 4250/3 b a - 46750/3 b a + 18700 a b - 23375 a b
12 12
+ 425/3 b - 425/3 a
]
> Scs:= seq(normal(Scoeffs[kk]/shouldbe[kk]), kk = 1 .. 5);
17
Scs := 0, -1/180, 1/1512, -1/14400, -------
2395008
# (d) Let's try it for x^13.
> f:= x -> x^13:
# We must repeat the definitions of J, ME, TE, SE etc. The coefficients should go up to jj=6 this time.
> J:= int(f(x), x=a..b):
> ME:= unapply(J - midrule(n), n):
> TE:= unapply(J - traprule(n), n):
> SE:= unapply(J - simpson(n), n):
> Mcoeffs:= [seq(coeff(expand(ME(n)),n,-2*jj), jj=1..6)]:
> Tcoeffs:= [seq(coeff(expand(TE(n)),n,-2*jj), jj=1..6)]:
> Scoeffs:= [seq(coeff(expand(SE(n)),n,-2*jj), jj=1..6)]:
> shouldbe := [ seq((b-a)^(2*jj)*((D@@(2*jj-1))(f)(b) - (D@@(2*jj-1))(f)(a)), jj = 1 .. 6) ]:
> Mcs:= seq(normal(Mcoeffs[kk]/shouldbe[kk]), kk = 1 .. 6);
> Tcs:= seq(normal(Tcoeffs[kk]/shouldbe[kk]), kk = 1 .. 6);
> Scs:= seq(normal(Scoeffs[kk]/shouldbe[kk]), kk = 1 .. 6);
31 127 73 1414477
Mcs := 1/24, -7/5760, ------, - ---------, ----------, - ----------------
967680 154828800 3503554560 2678117105664000
691
Tcs := -1/12, 1/720, -1/30240, 1/1209600, -1/47900160, -------------
1307674368000
17 21421
Scs := 0, -1/180, 1/1512, -1/14400, -------, - -----------
2395008 29719872000
# (e)
> a:= 0: b:= 1: f:= exp:
> J:= int(f(x), x=a..b):
> TE:= unapply(J - traprule(n), n);
exp(1) (1 + exp(1/n)) 1 + exp(1/n)
1/2 --------------------- - 1/2 ------------
exp(1/n) - 1 exp(1/n) - 1
TE := n -> exp(1) - 1 - --------------------------------------------
n
> formula:= n -> sum( Tcs[j] *(b-a)^(2*j)*((D@@(2*j-1))(f)(b) - (D@@(2*j-1))(f)(a))/ n^(2*j), j = 1 .. 5);
5
----- (2 j) (2 j - 1) (2 j - 1)
\ Tcs[j] (b - a) (D (f)(b) - D (f)(a))
formula := n -> ) ---------------------------------------------------------
/ (2 j)
----- n
j = 1
> Digits:= 30:
> seq(evalf(nn^12 * (TE(nn)-formula(nn))), nn = 1 .. 10);
Digits := 30
-9 -9 -9
.88554587380646207944*10 , .90226022060161953*10 , .905424963218222*10 ,
-9 -9 -9
.9065378753292*10 , .907053921534*10 , .90733448936*10 ,
-9 -9 -9 -9
.9075037465*10 , .9076136333*10 , .907688984*10 , .90774289*10
# It certainly appears to be approaching a limit, perhaps somewhere around .908 * 10^(-9).
#
# (f) Since all derivatives of f(x) = 1/(2 + cos(x)) are the same at 0 and at 2 Pi, the formula would say that
# the error is O(n^(-m)) for every m, i.e. it decreases faster than any power of m.
> n:= 20: a:= 0: b:= 2*Pi: f:= x -> 1/(2+cos(x)):
> J:= int(f(x), x=a .. b);
1/2
J := 2/3 3 Pi
> evalf(J - traprule(n));
-10
-.2640573654902402695*10
> evalf(J - midrule(n));
-10
.2640573654883181638*10
> evalf(J - simpson(n));
-5
.461370711086269911080003*10
# The Trapezoid and Midpoint Rule errors are almost exactly equal in magnitude, and both very small.
# The Simpson's Rule error is much larger. This is rather surprising at first sight. The point is that the
# terms that are the main contributions to the errors for most functions happen to be 0 for this f. The error
# that there is, goes rapidly to 0 as
# n grows. There is no reason to expect the Midpoint Rule error to be approximately -1/2 times the
# Trapezoid Rule error, as was the case for those terms that happened to be 0 in this case. So the Simpson's
# Rule error for n = 20, which is a 2/3 : 1/3 mixture of the Midpoint and Trapezoid Rule errors for n=10, is
# of the same order of magnitude as those errors. But these are much larger than the Midpoint and
# Trapezoid Rule errors for n=20, because the errors converge to 0 so rapidly as n -> infinity.
#
# 3.(a) I'll keep most of the definitions from question 2, changing what is needed. Other definitions come
# from Lesson 19. I'll use Digits = 20, because some of the errors are very small.
> a:= 1: b:= 3: n:= 'n': Digits:= 20:
> f:= x -> 1/(x^2 - 2*x + 5):
> J:= int(f(x), x=a .. b); Je:=evalf(J);
J := 1/8 Pi
Je := .39269908169872415481
> for jj from 1 to 5 do R[jj,1]:= evalf(traprule(2^jj)) od:
> for kk from 2 to 5 do
> for jj from kk to 5 do
> R[jj,kk]:= (4^(kk-1)*R[jj,kk-1]-R[jj-1,kk-1])/(4^(kk-1)-1) od od;
> [seq(R[jj,1],jj=1..5)]; [seq(R[jj,2],jj=2..5)];\
> [seq(R[jj,3],jj=3..5)]; [seq(R[jj,4],jj=4..5)];
> R[5,5];
[.38750000000000000000, .39139705882352941176, .39237356181138612617,
.39261770150517361183, .39267873664687180415]
[.39269607843137254900, .39269906280733836430, .39269908140310277370,
.39269908169410453493]
[.39269926176573608532, .39269908264282040100, .39269908171350465235]
[.39269907979959951713, .39269908169875360872]
.39269908170620127181
# (b) For int(x^11, x=0 .. 1):
> a:= 0: b:= 1: f:= x -> x^11:
> J:= int(f(x), x=a..b); Je:= evalf(J);
J := 1/12
Je := .083333333333333333333
> for jj from 1 to 5 do R[jj,1]:= evalf(traprule(2^jj)) od:
> for kk from 2 to 5 do
> for jj from kk to 5 do
> R[jj,kk]:= (4^(kk-1)*R[jj,kk-1]-R[jj-1,kk-1])/(4^(kk-1)-1) od od;
> seq([ Je - R[nn,nn], R[nn,nn-1]-R[nn,nn]], nn=2..5);
[-.014159838358561197919, .038187742233276367188],
[-.000346307953198750807, .000863345650335152949],
-5 -5
[-.1179364820321403*10 , .5392634193412958*10 ],
-9 -8
[-.388051072754*10 , .4605378004879*10 ]
# In each of these cases the actual error in R[n,n] is less in absolute value than R[n,n-1] - R[n,n].
> f:= x -> x^41: J:= int(f(x), x=a..b); Je:= evalf(J);
J := 1/42
Je := .023809523809523809524
> for jj from 1 to 5 do R[jj,1]:= evalf(traprule(2^jj)) od:
> for kk from 2 to 5 do
> for jj from kk to 5 do
> R[jj,kk]:= (4^(kk-1)*R[jj,kk-1]-R[jj-1,kk-1])/(4^(kk-1)-1) od od;
> seq([ Je - R[nn,nn], R[nn,nn-1]-R[nn,nn]], nn=2..5);
[-.059526323670175724408, .041666038130131959938],
[-.015824956803071526056, .002731335429194012396],
[-.001879668742061791717, .000217895125953277100],
-5
[-.000062025551328686623, .7100168713801191*10 ]
# In each of these cases the actual error in R[n,n] is greater in absolute value than R[n,n-1] - R[n,n].
# (c) No, it isn't, as we see above for k = 41.
# (d) R[n,n] is exactly correct for integrating powers up to x^(2*n-1). So for x^41 we would need
# R[21,21]. This would require calculating the Trapezoid Rule with up
# to 2^21 intervals.
> 2^21;
2097152
# That's not very practical.
