appendixa.mw

Appendix A

Maple code for exercises in section A.1

1.

 > wallis := n -> evalf( product( 4*k^2/(4*k^2-1),k=1..n));

 > Show that the average value of the upper and lower bounds on Pi can be computer as

 > wallisaverage := n -> 2*evalf( product( 4*k^2/(4*k^2-1),k=1..n)*(1+1/(4*n)) );

7.

 > evalf( binomial(1,1/2),20 );

 > evalf( 4/Pi, 20);

13.

 > evalf( binomial(1,1/3),20);

 > evalf(1/int((1-x^(3/2))^(1/3), x=0..1),20);

Maple code for exercises in section A.2

1.

 > The nth Bernoulli number is stored in Maple as bernoulli(n). The following command generates a list of the first forty Bernoulli numbers

 > [seq(bernoulli(n),n=1..40)];

 > We can now use the recursive formula for Bernoulli polynomials:

 > B := (n,x) -> if n=1 then x-1/2 else n*int(B(n-1,t),t=0..x)+bernoulli(n) end if;

 > [seq(B(n,x),n=1..10)];

2.

 > for n from 1 to 8 do plot(B(n,x),x=-1..2) end do;

8 & 9

 > The command NB(n) finds the numberator of the 2nth Bernoulli number. Thus, the numerator of B_20 is NB. The command DB(n) finds the denominator.

 > NB := n -> numer(bernoulli(2*n)); DB := n -> denom(bernoulli(2*n));

 > The following commands list the factorizations of the numerators and of the denoninators of the Bernoulli numbers. The first column gives the value of n (it is the Bernoulli number for 2n), the second column list the prime divisors, and the third column gives the powers of each of those prime divisors.

 > [seq([n,ifactor(NB(n))],n=1..20)];

 > [seq([n,ifactor(DB(n))],n=1..20)];

Maple code for exercises in section A.3

1.

 > The following command calculates the difference between the sum of the reciprocals of the first n perfect squares and Pi^2/6.

 > df := n -> evalf( sum(1/k^2,k=1..n) - Pi^2/6,10);

 > [seq( [100*n, df(100*n)],n=1..10)];

 > The next command adds 1/n to the sum of the reciprocals of the first n squares.

 > dfplus := n -> evalf( 1/n + sum(1/k^2,k=1..n) - Pi^2/6,10);

 > [seq( [100*n, dfplus(100*n)],n=1..10)];

2.

 > The following command calculates the difference between the sum of the reciprocals of the first n fourth and sixth powers and the exact values of the infinite series.

 > df4 := n -> evalf( sum(1/k^4,k=1..n) - Pi^4/90,10);

 > [seq( [100*n, df4(100*n)],n=1..10)];

 > df6 := n -> evalf( sum(1/k^6,k=1..n) - Pi^6/945,10);

 > [seq( [100*n, df6(100*n)],n=1..10)];

8.

 > P := (n,x) -> 1 + sum( bernoulli(k)*x^k/k!,k=1..n);

 > The following command superimposes the the polynomials that approximate x/(exp(x)-1).

 > plot([P(4,x),P(6,x),P(8,x),P(10,x),P(12,x),x/(exp(x)-1)],x=-8..8);

11.

 > The command given below will find these derivatives, but notice that the denominators are 0. Show that in the limit as x approaches 0, these all approach the correct values.

 > for n from 1 to 3 do simplify(diff(x/(exp(x)-1),x\$n)) end do;

14.

 > Q := (n,x) -> 1 + sum((-1)^k * bernoulli(2*k)*(2*x)^(2*k)/(2*k)!, k=1..n);

 > plot([Q(2,x),Q(4,x),Q(6,x),Q(8,x),x+cot(x)],x=-4..4);

17.

 > R := (n,z) -> 1/z + 2*sum( z/(z^2-k^2*Pi^2), k=1..n);

 > plot([R(3,x),R(6,x),R(9,x),R(12,x),cot(x)],x=-4..4,y=-1..1);

 > The following command plots the differences between the cotangent and each of thes seies. Use these graphs to give a reasonable approximation of the error function cot(x)-R(n,x)

 > for n from 1 to 4 do plot(cot(x) - R(3*n,x),x=-4..4) end do;

20.

 > If we evaluate the first 1000 terms of this series, then we obtain a value that is between

 > evalf(1/(2*1001^2), 10);

 > and

 > evalf(1/(2*1000^2), 10);

 > below the true value of zeta(3), and therfore zeta(3) lies in the interval

 > [evalf(sum(1/k^3,k=1..1000)+1/(2*1001^2),12),evalf(sum(1/k^3,k=1..1000)+1/(2*1000^2),12)];

Maple code for exercises in section A.4

2.

 > The command StirlingTriple produces three numbers: n, n! (to 10-digit accuracy), and Stirling's formula for n with the first two terms of the asymptotic series (to 10-digit accuracy).

 > StirlingTriple := x -> [x,evalf(x!,10),evalf( (x/exp(1))^x*sqrt(2*Pi*x)*exp(1/(12*x) - 1/(30*x^3)), 10)];

 > [seq(StirlingTriple(x),x in [5, 10, 20, 50, 100])];

3.

 > The following command finds the aboluate value of the mth summand in the asymptotic series, evaluated at n.

 > Summand := (m,n) -> abs(bernoulli(2*m)/((2*m)*(2*m-1)*n^(2*m-1)));

 > The next command lists the absolute values of the summands, taking 20 summands between 1 and 10n.

 > ListSummands := n -> [seq([floor(m*n/2), evalf(Summand(floor(m*n/2),n),10)],m=1..20)];

 > ListSummands(10);

 > To refine your search, the following commmands lists all the summands in the range a to b

 > RefineSearch := (n,a,b) -> [seq([m,evalf(Summand(m,n),10)],m=a..b)];

 > RefineSearch(10,25,40);

 > This command allows to compare the value of the approximation to n! with the best asymptotic estimate (to m terms):

 > Cf := (n,m) -> [n!,evalf( (n/exp(1))^n *sqrt(2*Pi*n)*exp(sum(bernoulli(2*k)/(2*k*(2*k - 1)*n^(2*k - 1)), k=1..m))),ceil(log10(n!))];

 > Cf(10,3);

5.

 > The first command calculates 1/2 + the sum of the first n terms. The second produces a table of these values as n goes from 1 to 10. The third command gives the actual value of Euler's gamma to 10-digit accuracy for purposes of comparison.

 > GammaApprox := n -> evalf( 1/2 + sum(bernoulli(2*k)/(2*k), k=1..n), 10);

 > [seq([n,GammaApprox(n)],n=1..10)];

 > evalf(gamma,10);

6.

 > The following command finds the aboluate value of the mth summand in the asymptotic series, evaluated at n.

 > HarmonicSummand := (m,n) -> abs( bernoulli(2*m)/(2*m*n^(2*m)) ) ;

 > The next commmand lists the absolute values of the summands, taking 20 summands between 1 and 10n.

 > ListHarmonicSummands := n -> [seq( [floor(m*n/2), evalf( Summand(floor(m*n/2),n), 10)],m=1..20)];

 > ListHarmonicSummands(10);

 > To refine your search, the following commmands lists all the summands in the range a to b

 > RefineHarmonicSearch := (n,a,b) ->  [seq([m,evalf( Summand(m,n),10)],m=a..b)];

 > RefineHarmonicSearch(10,25,40);

 > This command allows you to compare the value of the approximation to the harmonic series with the best asymptotic estimate (to m terms):

 > HarmonicCf := (n,m) -> [ evalf(sum((1/k),k=1..n),n), evalf( ln(n) + gamma + 1/(2*n) - sum(bernoulli(2*k)/(2*k*n^(2*k)), k=1..m) , n) ];

 > HarmonicCf(10,3);