# Solutions to Assignment 5
#
# 1. Since x exp(-x^2) has an elementary antiderivative, we write
# the integrand as u dv with u=1/x and dv = x exp(-x^2) dx.
> with(student):
> J:= exp(R^2) * Int(exp(-x^2), x=R .. inf
> inity);
> J1:= intparts(J, 1/x);
# Note that for large R, the integral term here will be smaller than
# J because of the 1/x^2. This time, take u = 1/x^3.
> J1:= intparts(J1,1/x^3);
# Now a pattern is starting to emerge: each time we get another x^2
# in the denominator.
> for jj from 3 to 5 do J1:= intparts(J1,
> 1/x^(2*jj-1)) od:
> J1;
# (b) After n iterations, the integral term is c(n) exp(-x^2)/x^(2 n),
# where c(n) is a constant. Integrating that by parts:
> intparts(c(n)*Int(exp(-x^2)/x^(2*n),x=R.
> .infinity), 1/x^(2*n+1));
# So c(n+1) = -(2n+1)/2 c(n). Since c(0) = 1, we can write
> c:= n -> (-1)^n*2^(-n) * Product(2*k+1,
> k = 1 .. n-1);
# Maple actually can find a "closed form" for this product.
> value(c(n));
> c1:= simplify(");
# Alternatively, note that c(n) contains the product of the odd
# integers from 1 to 2n-1. The product of all the integers from 1 to
# 2n is (2n)!, while the product of the even ones is
# 2^n n!. So we can also write c(n) as
> c2:= (-1)^n*2^(-2*n)*(2*n)!/n!;
> simplify(c2/c1);
> c:= unapply(c2,n);
> seq(c(nn), nn = 1 .. 5);
# Whichever way we write it, the n'th term of the series is
> nthterm:=unapply( c(n-1)/2 * R^(-2*n+1),
> n);
> theseries:= Sum(nthterm(n), n=1..infinit
> y);
# Check this for the first 5 terms.
> partialsum:= n -> sum(nthterm(k), k=1..n
> );
> partialsum(5);
# To check for convergence, use the Ratio Test.
> nthterm(n+1)/nthterm(n);
> simplify(");
# For any R the absolute value of this goes to infinity as n -> infinity,
# so the series diverges.
#
# (c) The remainder after 5 terms is
> R5:= expand(J1 - partialsum(5));
# Another integration by parts yields
> intparts(R5,x^(-11));
> expand(");
# Since the integrand is positive, we conclude that
# 0 > R5 > -945/64*R^(-11). For R=2 we have
> evalf(945/64 * 2^(-11));
# We calculated nthterm(n+1)/nthterm(n) = - (2 n - 1)/(2 R^2). So
# |nthterm(n+1)| >= |nthterm(n)| when 2 n - 1 >= 2 R^2. For R = 2
# this means n >= 5.
> seq(evalf(abs(subs(R=2,nthterm(nn)))), n
> n=1..6);
# The remainder is smaller in absolute value than the next term of the
# series. This estimate for the error is smallest for the fourth partial
# sum.
# Let's check the actual errors.
> V:= value(J);
> seq(evalf(subs(R=2,V - partialsum(nn))),
> nn = 1 .. 5);
#
# 2. We'll try it with the partial sum to n=10, and up to 10 terms of
# the Euler-Maclaurin series.
> f:= n -> 1/(n^2+1):
> PartialSum:= N -> Sum(f(n), n = 1 .. N):
> readlib(eulermac): Digits:= 20:
> p:= evalf(PartialSum(10));
> ems:= [seq(eulermac(f(N),N,2*dd-1), dd=
> 1..10)]:
> ems[1];
# Note that this has a limit of Pi/2 as N -> infinity. So if the sum of
# the series is S, F(N) = PartialSum(N-1) - S + Pi/2. Our
# approximation
# for S will be PartialSum(10) - F(11) + Pi/2).
> ems:= evalf(eval(subs(N=11,O=0,ems)));
# S[3] to S[10] all agree to within 10^(-9). The best value might be
# S[9], since it and S[10] agree to the accuracy of our calculation.
# So our answer would be
> p - ems[9]+evalf(Pi/2);
# Compare to Maple's value for the sum.
> sum(f(n), n=1..infinity);
> evalf(");
# 3. (a)
> read `fib.def`:
Warning: new definition for norm
Warning: new definition for trace
mpow :=
proc(n)
options remember;
if type(n,posint) then
if type(n,even) then
evalm(mpow(1/2*n)^2)
else
evalm(mpow(1/2*n-1/2)^2 &* M
)
fi
elif type(n,negint) then
evalm(1/mpow(-n))
elif type(n,`+`) then
evalm(`&*`(op(map(mpow,n))))
elif type(n,`*`) then
evalm(mpow(n/op(1,n))^op(1,n))
else
matrix([
[F[n-1],F[n]],[F[n],F[n]+F[n-1]]
])
fi
end
zeck := proc(x)
local n;
if x = 0 then 0
else
n := invfib(x);
F[n]+zeck(x-fib(n))
fi
end
> fib(n+4)+fib(n-4);
# So the equation is F[n+4] + F[n-4] = 7 F[n].
# (b) The identity was
> fib(m+n);
# We conclude that F[m+n] is divisible by F[n] if F[m] is. Since
# F[0] = 0 is divisible by F[n], we conclude by mathematical
# induction that F[k n] is divisible by F[n] for all positive integers k.
#
# (c) Since F[4] = 3, all the Fibonacci numbers F[4 k] are divisible
# by 3. What's the first one divisible by 9?
> seq(fib(4*kk)/9, kk = 1 .. 10);
# So F[12] is divisible by 9, and then F[12 k] is divisible by 9 for all
# positive integers k. Which are divisible by 3^3 = 27?
> seq(fib(12*kk)/27, kk = 1 .. 10);
# F[36] is the first, and then of course F[36 k] for all k.
# 4, 12, 36: is there a pattern here? It looks like F[4*3^m] will be
# divisible by 3^(m+1). If F[n] is divisible by 3^m, is F[3*n]
# divisible by 3^(m+1)?
> fib(3*n);
> subs(F[n] = 3^m * r, ");
> "/3^(m+1);
> simplify(");
# Yes: if m >= 1 then this is an integer. We conclude that
# F[4*3^m*k] is divisible by 3^(m+1) for all nonnegative integers m
# and k.
#
# (d) If F[n]/F[m] = 7880997 and n and m are not too small, we
# can approximate F[n] and F[m] by fapp(n) and fapp(m) to get
> fapp(n)/fapp(m);
# So 7880997 should be approximately ((1+sqrt(5))/2)^(n-m).
# Thus n - m should be approximately
> evalf(ln(7880997)/ln((1+sqrt(5))/2));
# It looks very much like n-m = 33. If so, .
> shouldbe:= fib(m+33)/fib(m)= 7880997;
> solve(", F[m]);
# So F[m] should be 89/55*F[m-1]. In particular, F[m-1] should be
# divisible by 55. Is 55 itself a Fibonacci number?
> invfib(55);
> fib(10);
> fib(11);
# So we want m=11, making F[m] = 89 and F[m-1]=55, and
# n=11+33=44.
> fib(44)/fib(11);
#
# 4. We try a Chebyshev series.
> with(numapprox): Digits:= 15:
> f: x -> ln(1+x^3):
> ff:= x -> evalf(ln(1+x^3)):
> chebyshev(ff,x=0 .. 3, 10^(-5));
> ch:= unapply(eval(subs(T=orthopoly[T],")
> ),x);
> plot(ln(1+x^3)-ch(x), x = 0 .. 3);
# According to the graph, the maximum error is approximately
# 3.4*10^(-6) , which is acceptable.
# Note: you could expand ch(x) into a sum of coefficients times
# powers of x, but this would
# be bad for numerical accuracy when x is near 3, so it's better not
# to.
