# Math 210 Sec. 201
# Solutions to Assignment 4
#
# # 1.(a) Consider the integral
#
# 1
# /
# | 1/2
# | x dx
# |
# /
# 0
#
# # with ME(n), TE(n) and SE(n) the errors for Midpoint, Trapezoid and Simpson's Rules.
# # Plot n^(3/2) ME(n), n^(3/2) TE(n), and n^(3/2) SE(n) (for even n) on the same graph.
# # What can you conclude?
#
> f:= x -> x^(1/2): a:= 0: b:= 1: J:= int(f(x), x=a..b);\
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): \
ME:= n -> J - midrule(n): TE:= n -> J - traprule(n): SE:= n -> J - simpson(n):
J := 2/3
> plot({[seq([nn, ME(nn)*nn^(3/2)], nn = 1 .. 20)], \
[seq([nn, TE(nn)*nn^(3/2)], nn = 1 .. 20)],\
[seq([2*nn, SE(2*nn)*(2*nn)^(3/2)], nn= 1 .. 10)]});
# All three appear to approach nonzero constants as n -> infinity. Conclude that all three
# errors are of order n^(-3/2), rather than n^(-2) or n^(-4) as would happen for a smooth
# integrand. The fact that x^(1/2) is not differentiable at x=0 is the problem.
#
# # (b) Approximately what n would be needed to make the Simpson's Rule error be less
# # than 10^(-9)?
#
#
# According to the graph, S(n) is approximately .08 n^(-3/2). This would be 10^(-9)
# at the following value of n:
> (.08/10^(-9))^(2/3);
185663.5533
# So the minimum n would be approximately 186000.
#
# # (c) Transform
#
# infinity
# /
# | 1
# | --------------- dx
# | 3 1/2
# / x + (x + 1)
# 0
#
# # into a proper integral with a smooth integrand, and use Simpson's Rule to find an
# # approximation with error less than 10^(-6).
#
# The integral is improper only at infinity (the square root causes no problem for x > 0, nor
# is the denominator ever 0). The change of variables t = 1/(x+1) maps 0 < x < infinity to
# 0 < t < 1.
> with(student, changevar):\
J:= Int(1/(x+(x^3+1)^(1/2)), x=0..infinity);\
infinity
/
| 1
J := | --------------- dx
| 3 1/2
/ x + (x + 1)
0
> J2:= changevar(t=1/(x+1), J, t);
1
/
| 1
J2 := | ------------------------------ dt
| / / 3 \1/2\
/ |1 - t |(1 - t) | | 2
0 |----- + |-------- + 1| | t
| t | 3 | |
\ \ t / /
# Some simplification must be done, to see whether we have a singularity at t=0.
> simplify(J2);
1
/
| 1
| --------------------------------- dt
| / / 2\1/2 \
/ | |1 - 3 t + 3 t | |
0 t |1 - t + |--------------| t|
| | 3 | |
\ \ t / /
# Maple isn't simplifying this as much as we'd want, because it's reluctant to take the t^3
# out of the square root. Since t is positive, this would be all right. The "symbolic" option
# to "simplify" can be used here.
> simplify(", symbolic);
1
/
| 1
| ----------------------------- dt
| 2 2 3 1/2
/ t - t + (t - 3 t + 3 t )
0
# Yes, there is a problem. The integrand has a 1/sqrt(t) singularity at 0. A further
# transformation may smooth this out: t = s^2.
> J3:= changevar(t=s^2, ", s);
1
/
| s
J3 := | 2 ------------------------------- ds
| 2 4 2 4 6 1/2
/ s - s + (s - 3 s + 3 s )
0
> J4:= simplify(", symbolic);
1
/
| 1
J4 := - 2 | ------------------------------- ds
| 3 2 4 1/2
/ - s + s - (1 - 3 s + 3 s )
0
> f:= s -> 2/(s - s^3 + (1 - 3*s^2 + 3*s^4)^(1/2));
2
f := s -> -----------------------------
3 2 4 1/2
s - s + (1 - 3 s + 3 s )
> a:= 0: b:= 1:
> S1:= evalf(simpson(10)); S2:= evalf(simpson(20));
S1 := 2.013966056
S2 := 2.013929309
# Richardson extrapolation would indicate that the error in S2 is approximately
> Eest:= abs(S2-S1)/15;
-5
Eest := .2449800000*10
# a bit larger than we want. But n=40 should do it.
> S3:= evalf(simpson(40));
S3 := 2.013926945
> Eest:= abs(S3-S2)/15;
-6
Eest := .1576000000*10
# and we might as well use the improved approximation
> R:= (16*S3 - S2)/15;
R := 2.013926788
# Just for comparison, here's Maple's value for the integral.
> evalf(J);
2.013926789
#
# # 2.(a) Catalan's constant is defined as
#
#
# infinity
# ----- i
# \ (-1)
# C = ) ----------
# / 2
# ----- (2 i + 1)
# i = 0
#
#
# # By finding bounds on the tail of this series, determine the value of C with an error less
# # than 10^(-8). Compare to Maple's "Catalan". Hint: consider the terms in pairs:
#
#
#
# infinity
# -----
# \
# C = ) (a + a )
# / 2j 2j+1
# -----
# j = 0
#
#
# The pair for j contains the terms for i=2*j and i=2*j+1, so we have
> Sum(1/(4*j+1)^2 - 1/(4*j+3)^2, j=0 .. infinity);
infinity
-----
\ / 1 1 \
) |---------- - ----------|
/ | 2 2|
----- \(4 j + 1) (4 j + 3) /
j = 0
> f:= unapply(normal(1/(4*j+1)^2 - 1/(4*j+3)^2),j);
2 j + 1
f := j -> 8 ---------------------
2 2
(4 j + 1) (4 j + 3)
# This is a decreasing function of j, so the Integral Test estimates should apply. Upper
# and lower bounds for the "tail", the sum from N+1 to infinity, are
> upperbd:= int(f(x), x=N .. infinity);
1
upperbd := ---------------------
2 (4 N + 1) (4 N + 3)
> lowerbd:= int(f(x), x=N+1 .. infinity);
1
lowerbd := ---------------------
2 (4 N + 5) (4 N + 7)
# The difference between these is
> difference:=normal(upperbd - lowerbd);
N + 1
difference := 16 ---------------------------------------
(4 N + 1) (4 N + 3) (4 N + 5) (4 N + 7)
# Our approximation for the tail will be the average of the upper and lower bounds, which
# will have error at most half of this. For what N will the estimated error be 10^(-8)?
> fsolve(difference/2 = 10^(-8), N);
-.9999999888
# Silly of me, I didn't tell Maple that N was positive.
> fsolve(difference/2 = 10^(-8), N, N = 1 .. infinity);
145.2023119
> MyCatalan:= evalf(Sum(f(j), j = 0 .. 146) + subs(N=146, upperbd+lowerbd)/2);
MyCatalan := .9159655942
> evalf(Catalan);
.9159655942
>
# # (b) Euler's constant gamma is defined by
#
#
# / n \
# |----- |
# | \ |
# gamma = limit | ) 1/i - ln(n)|
# n -> infinity | / |
# |----- |
# \i = 1 /
#
#
# # By finding bounds on the tail of the convergent series
#
#
# infinity
# -----
# \
# ) (1/i + ln(i) - ln(i + 1))
# /
# -----
# i = 1
#
#
# # determine the value of gamma with an error less than 10^(-8). Compare to Maple's
# # "gamma".
#
# Note that the partial sum to i=N is
> p:= N -> Sum(1/i, i=1..N) - ln(N+1);
/ N \
|----- |
| \ |
p := N -> | ) 1/i| - ln(N + 1)
| / |
|----- |
\i = 1 /
> f:= x -> 1/x + ln(x) - ln(x+1);
f := x -> 1/x + ln(x) - ln(x + 1)
# This is positive and decreasing, with limit 0.
> normal(diff(f(x),x));
1
- ----------
2
x (x + 1)
# So the Integral Test bounds apply. The upper and lower bounds for the tail are
> upperbd:= int(f(x), x=N .. infinity);
upperbd := - ln(N) + ln(N + 1) N - ln(N) N + ln(N + 1) - 1
> lowerbd:= int(f(x), x=N+1 .. infinity);
lowerbd := - 2 ln(N + 1) + ln(N + 2) N - ln(N + 1) N + 2 ln(N + 2) - 1
> difference:=upperbd - lowerbd;
difference :=
- ln(N) + 2 ln(N + 1) N - ln(N) N + 3 ln(N + 1) - ln(N + 2) N - 2 ln(N + 2)
> fsolve(difference/2 = 10^(-8), N, N = 1 .. infinity);
4999.122055
# A bit large for my taste! Let's try the technique from Lesson 22, based on the Trapezoid
# Rule. We required (D@@2)(f)(t) to be decreasing. To test this, look at (D@@3)(f)(t).
> normal((D@@3)(f)(t));
2
6 t + 8 t + 3
- 2 --------------
4 3
t (t + 1)
> Approx:= (f,N) -> Sum(f(n),n=1..N) \
+ Int(f(t),t=N .. infinity) - f(N)/2 \
- (D(f)(N+1)+D(f)(N-1))/24;
/ N \ infinity
|----- | /
| \ | |
Approx := (f,N) -> | ) f(n)| + | f(t) dt - 1/2 f(N)
| / | |
|----- | /
\n = 1 / N
- 1/24 D(f)(N + 1) - 1/24 D(f)(N - 1)
> ErrEst:= (f,N) -> (D(f)(N+1)-D(f)(N-1))/24;
ErrEst := (f,N) -> 1/24 D(f)(N + 1) - 1/24 D(f)(N - 1)
> ErrEst(f,N);
1 1 1 1 1 1
- ----------- + ---------- - ---------- + ----------- - ---------- + ----
2 24 (N + 1) 24 (N + 2) 2 24 (N - 1) 24 N
24 (N + 1) 24 (N - 1)
> fsolve(" = 10^(-8), N = 1 .. infinity);
70.39108520
> evalf(Approx(f, 71));
.577215665
> evalf(ErrEst(f,71));
-8
.9662771958*10
> evalf(gamma);
.5772156649
>
#
# # 3.(a) Let
#
#
# 2 Pi
# /
# 1 |
# f(x) = ---- | exp(x cos(t)) dt
# 2 Pi |
# /
# 0
#
#
# # Maple can't evaluate this integral, but in fact it turns out to be a Maple function:
# # "BesselI(0,x)". Evaluate f and "BesselI(0,...)" at several points to convince yourself that
# # they are the same.
#
> f:= x -> 1/(2*Pi)*Int(exp(x*cos(t)), t=0 .. 2*Pi);
2 Pi
/
|
| exp(x cos(t)) dt
|
/
0
f := x -> 1/2 ----------------------
Pi
> evalf([seq(f(kk) = BesselI(0,kk), kk = -1 .. 3)]);
[1.266065878 = 1.266065878, .9999999995 = 1., 1.266065878 = 1.266065878,
2.279585302 = 2.279585302, 4.880792584 = 4.880792586]
#
# # (b) Find the Taylor series of f(x) at x=0 by starting with the Taylor series of exp. Hint:
# 2 Pi
# /
# | n
# | cos (t) dt = 0 if n is odd
# |
# /
# 0 (1 - n)
# 2 Pi n!
# = -------------- if n is even
# 2
# ((n/2)!)
#
> exp(x*cos(t)) = Sum((x*cos(t))^n/n!, n = 0 .. infinity);
infinity
----- n
\ (x cos(t))
exp(x cos(t)) = ) -----------
/ n!
-----
n = 0
> J := 1/(2*Pi)*int(op(2,"), t = 0 .. 2*Pi);
2 Pi infinity
/ ----- n
| \ (x cos(t))
| ) ----------- dt
| / n!
/ -----
0 n = 0
J := 1/2 -----------------------------
Pi
# It's not so simple to get Maple to interchange the integral and the sum, so I'll just do it
# by hand.
> J:= 1/(2*Pi)*sum(Int((x*cos(t))^n/n!, t=0..2*Pi), n=0 .. infinity);
infinity 2 Pi
----- / n
\ | (x cos(t))
) | ----------- dt
/ | n!
----- /
n = 0 0
J := 1/2 -----------------------------
Pi
# Maple also won't do the integral of cos(t)^n in general (although it does do it for any
# given value of n.
> int(cos(t)^n, t=0..2*Pi);
2 Pi
/
| n
| cos(t) dt
|
/
0
> seq(int(cos(t)^nn, t=0 .. 2*Pi) , nn = 0 .. 6);
2 Pi, 0, Pi, 0, 3/4 Pi, 0, 5/8 Pi
> seq(2^(1-2*kk)*Pi*(2*kk)!/(kk!)^2, kk = 0 .. 3);
2 Pi, Pi, 3/4 Pi, 5/8 Pi
> J:= expand(1/(2*Pi)*Sum(x^(2*k)/(2*k)! * 2^(1-2*k)*Pi*(2*k)!/(k!)^2, k = 0 .. infinity));
infinity
----- k 2
\ (x )
J := ) -----------
/ 2 k 2
----- (k!) (2 )
k = 0
#
# # (c) Check your answer for (b) by using "taylor" on f(x) and BesselI(0,x).
> expand(1/(2*Pi)*sum(x^(2*k)/(2*k)! * 2^(1-2*k)*Pi*(2*k)!/(k!)^2, k = 0 .. 5));
2 4 6 8 10
1 + 1/4 x + 1/64 x + 1/2304 x + 1/147456 x + 1/14745600 x
> taylor(f(x), x=0, 11);
2 4 6 8 10 11
1 + 1/4 x + 1/64 x + 1/2304 x + 1/147456 x + 1/14745600 x + O(x )
> taylor(BesselI(0,x), x=0, 11);
2 4 6 8 10 11
1 + 1/4 x + 1/64 x + 1/2304 x + 1/147456 x + 1/14745600 x + O(x )
>
# # (d) Study the convergence of the Taylor series to BesselI(0,x), using the methods we
# # used for exp in Lesson 23. Are the results similar?
> P:= (n,t) -> subs(x=t, convert(taylor(BesselI(0,x), x=0, n+1), polynom));
P := (n,t) -> subs(x = t, convert(taylor(BesselI(0, x), x = 0, n + 1), polynom))
# Of course we only need to consider P(n,t) for even n, because the Taylor series has no odd
# terms.
> plot({ BesselI(0,x), seq(P(2*nn, x), nn=1 .. 5) }, x=-5 .. 5, y=0 .. 20, colour = black);
> plot({ seq(BesselI(0,x) - P(2*nn, x), nn = 1 .. 10) }, x = 0 .. 10, y = 0 .. 20, colour=black);
# A similar effect: the curves seem to march off to the right at regular intervals. The
# horizontal distance between the curves is approximately 2/E this time.
> plot({ seq(BesselI(0,x+2*nn/E) - P(2*nn, x+2*nn/E), nn = 1 .. 20) },\
x=-1 .. 3, y = 0 .. 3, color=black);
>
# # (e) The function f(x) satisfies the differential equation x f''(x) + f'(x) - x f(x) = 0. What
# # is x S[10]''(x) + S[10]'(x) - x S[10](x), where S[10] is the Taylor polynomial of degree
# # 10?
> P(10,x);
2 4 6 8 10
1 + 1/4 x + 1/64 x + 1/2304 x + 1/147456 x + 1/14745600 x
> x*diff(P(10,x),x$2) + diff(P(10,x),x) - x*P(10,x);
2 4 6 8 3
x (1/2 + 3/16 x + 5/384 x + 7/18432 x + 1/163840 x ) + 1/2 x + 1/16 x
5 7 9
+ 1/384 x + 1/18432 x + 1/1474560 x
2 4 6 8 10
- x (1 + 1/4 x + 1/64 x + 1/2304 x + 1/147456 x + 1/14745600 x )
> expand(");
11
- 1/14745600 x
>
#
# # 4.(a) Find the Fourier series for the function f(x) defined on [-Pi, Pi] by f(x) = 1 for x >
# # 0 and -1 for x < 0.
#
# The function is odd, so there will be only sin terms, and the coefficient can be found by
# integrating on [0,Pi].
> b:= unapply( 2/Pi*int(sin(n*t), t=0 .. Pi),n);
cos(Pi n)
- --------- + 1/n
n
b := n -> 2 -----------------
Pi
# Maple doesn't know the general formula for cos(Pi n) when n is an integer. Note that
# b(n) is 0 when n is even, and 4/(n*Pi) when n is odd.
> seq(b(nn), nn = 1 .. 6);
4 4 4
----, 0, ----, 0, ----, 0
Pi 3 Pi 5 Pi
> S:= Sum(4/((2*k+1)*Pi) * sin((2*k+1)*x), k = 1 .. infinity);
infinity
-----
\ sin((2 k + 1) x)
S := ) 4 ----------------
/ (2 k + 1) Pi
-----
k = 1
# # (b) Produce an animation of the partial sums, so that the Gibbs phenomenon can be
# # seen. Near what points does this occur?
> PartialSum:= n -> Sum(4/((2*k+1) *Pi) * sin((2*k+1)*x), k = 0 .. n-1);
n - 1
-----
\ sin((2 k + 1) x)
PartialSum := n -> ) 4 ----------------
/ (2 k + 1) Pi
-----
k = 0
> value(PartialSum(3));
sin(x) sin(3 x) sin(5 x)
4 ------ + 4/3 -------- + 4/5 --------
Pi Pi Pi
> frame:= n -> plot(PartialSum(n), x = 0 .. 2*Pi):
> with(plots,display): display([seq(frame(nn), nn = 1 .. 8)], insequence=true);
# The Gibbs phenomenon occurs near the integer multiples of Pi, which are the points
# where
# f is discontinuous.
#
# # (c) Show that the ratio of the size of the "overshoot" to the size of the jump is the same
# # as it was for the example we studied in Lesson 26.
> diff(PartialSum(n), x);
n - 1
-----
\ cos((2 k + 1) x)
) 4 ----------------
/ Pi
-----
k = 0
> value(");
sin(x n) cos(x n)
4 -----------------
sin(x) Pi
# So this will be 0 when sin(n x) = 0 or cos(n x) = 0. The highest peak seems to be the first
# one to the right of x=0, which will be at x = Pi/(2 n) where cos(n x) = 0.
> S:= subs(x=Pi/(2*n), PartialSum(n));
n - 1 (2 k + 1) Pi
----- sin(1/2 ------------)
\ n
S := ) 4 ---------------------
/ (2 k + 1) Pi
-----
k = 0
# By analogy with the case done in Lesson 26, we might expect this to be a Riemann sum
# for an integral. If I take t = (k + 1/2)/n, we are adding 2/(Pi n) sin(Pi t)/t for n equally
# spaced values of t from 1/(2 n) to 1 - 1/(2 n). The Midpoint Rule approximation for
# int(sin(Pi t)/t, t=0 .. 1) with n intervals is as follows:
> g:= t -> sin(Pi * t)/t;
sin(Pi t)
g := t -> ---------
t
> Mrule:= 1/n * Sum(g((k+1/2)/n), k = 0 .. n-1);
n - 1 Pi (k + 1/2)
----- sin(------------) n
\ n
) -------------------
/ k + 1/2
-----
k = 0
Mrule := -------------------------
n
> simplify(S/Mrule);
2
----
Pi
# So as n -> infinity, S should approach a limit of
> 2/Pi * int(g(t), t=0..1);
Si(Pi)
2 ------
Pi
# Since f(x) = 1 for x > 0 and -1 for x < 0, the "overshoot" is
> " - 1;
Si(Pi)
2 ------ - 1
Pi
# and the ratio of overshoot to height of jump is
> evalf("/2);\
.0894898721
# which is the same as we found for the example in Lesson 26.
#
