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 aswallisaverage := n -> 2*evalf( product( 4*k^2/(4*k^2-1),k=1..n)*(1+1/(4*n)) );

evalf( binomial(1,1/2),20 );evalf( 4/Pi, 20);

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

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)];

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

The command NB(n) finds the numberator of the 2nth Bernoulli number. Thus, the numerator of B_20 is NB[10]. 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)];

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)];

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)];

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);

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;

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);

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;

If we evaluate the first 1000 terms of this series, then we obtain a value that is betweenevalf(1/(2*1001^2), 10);andevalf(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)];

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])];

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 bRefineSearch := (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);

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);

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 bRefineHarmonicSearch := (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);