HW 5 Solution

Be sure to evaluate the cells in order. It would be good practice to predict what each cell will evaluate to and why, before executing it. However, you may run all the cells at once by selecting all (Ctrl + A) and executing with Shift + Enter

Problem 1

Part a

Assign the number to a variable so we don' t have to type it several times

(Subscript[n, a]) =190392490709135;

Brute force search

Fibonacci[50] > (Subscript[n, a])(* too small *)
Fibonacci[100] > (Subscript[n, a])(* too large *)
Fibonacci[70] > (Subscript[n, a]) (* too small *)
Fibonacci[80] > (Subscript[n, a]) (* too large *)
Fibonacci[75]>(Subscript[n, a]) (* too large *)
Fibonacci[72] > (Subscript[n, a]) (* too large *)
Fibonacci[71] > (Subscript[n, a]) (* too large *)

(* look at the bracketing values *)
Fibonacci[70]
(Subscript[n, a])
Fibonacci[71]

Using two of the criteria provided by Wikipedia, and Mathematica's help to determine which function calls to use. It is amazing that Mathematica can reduce the first expression to 0.

Fibonacci[Floor[Log[GoldenRatio, Sqrt[5] (Subscript[n, a])]+1/2]] -(Subscript[n, a]) (* this is zero, so (Subscript[n, a]) is a Fibonacci number *)
Sqrt[{5((Subscript[n, a])^2)+4, 5 ((Subscript[n, a])^2)-4}] (* the first term is a perfect square, so (Subscript[n, a]) is a Fibonacci number *)

All these tests tell us (Subscript[n, a]) is a Fibonnaci number, and the first tells us it is (Subscript[F, 70])

Part b

Assign the number to a variable so we don' t have to type it several times

(Subscript[n, b])= 1500520536206896083287;

Brute force search

Fibonacci[100] > (Subscript[n, b])(* too small *)
Fibonacci[200] > (Subscript[n, b]) (* too large *)
Fibonacci[150] > (Subscript[n, b]) (* too large *)
Fibonacci[125] > (Subscript[n, b]) (* too large *)
Fibonacci[110] > (Subscript[n, b]) (* too large *)
Fibonacci[105] > (Subscript[n, b]) (* too large *)
Fibonacci[102] > (Subscript[n, b]) (* too small *)
Fibonacci[103] > (Subscript[n, b]) (* too small *)
Fibonacci[104] > (Subscript[n, b]) (* too large *)

(* look at bracketing values *)
{Fibonacci[103],(Subscript[n, b]), Fibonacci[104]}

Using two of the criteria provided by Wikipedia, and Mathematica's help to determine which function calls to use.

Fibonacci[Floor[Log[GoldenRatio, Sqrt[5] (Subscript[n, b])]+1/2]] -(Subscript[n, b]) (* this is nonzero, so (Subscript[n, b]) is not a Fibonacci number *)
Sqrt[{5((Subscript[n, b])^2)+4, 5 ((Subscript[n, b])^2)-4} ](* since neither of these are perfect squares, (Subscript[n, b]) is not a Fibonacci number *)

All these tests tell us (Subscript[n, b]) is not a Fibonnaci number.

Part c

(Subscript[n, c])=8222838654177922817725562880000000;

Brute force (can you find an alternative test?)

Factorial[50] > (Subscript[n, c]) (* too large *)
Factorial[30] > (Subscript[n, c]) (* too small *)
Factorial[40] > (Subscript[n, c]) (* too large *)
Factorial[35] > (Subscript[n, c]) (* too large *)
Factorial[32] > (Subscript[n, c]) (* too large *)
Factorial[31] > (Subscript[n, c]) (* too small *)

(* look at bracketing values *)
{Factorial[31], (Subscript[n, c]), Factorial[32]}

It is clear that (Subscript[n, c])=31!

Part d

(Subscript[n, d])= 334525266131638071081700620534407516651520000000000;

Brute force, using a slightly different, more compact notation

values = { 50,40,45,42,41}
Select[values, Factorial[#] > (Subscript[n, d]) &] (* we will discuss anonymous functions in Lecture 7; this expression returns the values k for which (Subscript[F, k ])> (Subscript[n, d]) *)

(* look at bracketing values *)
{Factorial[42], (Subscript[n, d]), Factorial[41]}

It is clear that (Subscript[n, d]) is not a factorial.

Part e

Be careful that in answering (a) - (d) you didn' t define the variable n, or the following expression will evaluate to a number instead of a symbolic expression.

Sum[k Cos[k], {k,1,n}]

Mathematica can simplify this expression

FullSimplify@Out[5]

You can either set the cell style to Traditional or use the TraditionalForm command to pretty print the result

Out[6]//TraditionalForm

Altogether at once (the way you might actually make this computation in practice) :

Clear[n] (* Clear any value assigned to n so the following expression is symbolic *)
TraditionalForm@FullSimplify@Sum[k Cos[k], {k,1,n}]

Problem 2

Part a

Define expressions for A and A' to reduce typing

A = ((n (π^(n/2)))/(Gamma[n/2+1]));
Ap = ((n (π^(n/2)))/((Sqrt[2π])(E^(-n/2))((((n/2)))^((n+1)/2))));

By default, the first dataset is plotted in blue, and the second in purple, but it doesn' t hurt to specify the colors using the PlotStyle directive (this is one of those commands which is best learned by looking at examples in Mathematica' s help)

ListPlot[{Table[A, {n,1,20}], Table[Ap, {n,1,20}]}, PlotStyle->{{Blue}, {Purple}}]

From the plot, it is clear that A' is larger, but it doesn' t hurt to check

Table[Ap > A, {n,1,20}]

Part b

Filling is another directive that is best learned by example. The specification of Filling is almost incomprehensible, but the provided examples show how to do the most common tasks, such as filling the area between two graphs.

Plot[{A, Ap}, {n,1,20}, Filling->{1->{2}}]

Part c

Limit[(A/Ap),n->∞]

Part d

In this case, TraditionalForm is more readable. To see that ψ is the PolyGamma function, either don't use TraditionalForm, or highlight the symbol ψ, press F1 and read the documentation for the symbol.

D[A, n] //TraditionalForm

Part e

Solve[D[A,n]==0, n] (* useless *)
Reduce[D[A,n] ==0, n] (* useless *)
NSolve[D[A,n]==0,n] (* useless *)
FindRoot[D[A,n]==0, {n,2}] (* nonsensical answer; probably a bad starting point *)
FindRoot[D[A,n]==0, {n,15}] (* nonsensical answer; note the failure notice; probably a bad starting point *)

(* Let's use the graphs we created earlier to guess a starting point *)
FindRoot[D[A,n]==0,{n,5}]

Solve and Reduce simply aren' t suited for solving highly nonlinear equations, and NSolve is only for solving systems of polynomial equations. FindRoot is pretty general purpose, but you have to provide a sensible starting point.
   
It seems that the maximum value of A occurs around x = 7.25695. Of course this isn' t guaranteed, as we didn' t justify using the derivative to determine the maximum of A, and we used a numerical method to find the root of the derivative.

Part f

(*
(* most convenient method of generating the list to be plotted, but somewhat esoteric *)
ListPlot@Table[Evaluate@D[A, n], {n,1,20}] 
*)
ListPlot@ Table[((π^(n/2))/(Gamma[1+(n/2)]))+((n (π^(n/2)) Log[π])/(2 Gamma[1+(n/2)]))-((n (π^(n/2)) PolyGamma[0,1+(n/2)])/(2 Gamma[1+(n/2)])), {n,1,20}] (* copy and paste *)

From our knowledge of the (continuous) derivative, it is clear that A reaches its maximum as a function of discrete n at n=7

Problem 3

Part a

4 Mean[Table[If[RandomReal[]^2 + RandomReal[]^2 < 1, 1, 0],{(10^5)}]] //N

Note, if you wanted multiple approximations, you could wrap this Monte Carlo calculation with a Table (see below for an example).

Part b

6/(Mean@Table[If[RandomInteger[{2,(10^5)}]~CoprimeQ~RandomInteger[{2,(10^5)}], 1,0], {(10^5)}]) //Sqrt //N (* one possible compact form *)

(* slightly less compact form *)
approxπ = Mean@Table[
If[ RandomInteger[{2,(10^5)}] ~ GCD ~ RandomInteger[{2,(10^5)}] == 1,1,0],
{(10^5)}];
(Sqrt[(6/(approxπ))])//N

Let' s see how good these two approximations are

squareπlist = Table[4 Mean[Table[If[RandomReal[]^2 + RandomReal[]^2 < 1, 1, 0],{(10^5)}]] , {20}];
integerπlist = Table[6/(Mean@Table[If[RandomInteger[{2,(10^5)}]~CoprimeQ~RandomInteger[{2,(10^5)}], 1,0], {(10^5)}]) //Sqrt, {20}];
πlist = ConstantArray[π,20];

ListPlot[{squareπlist, integerπlist, πlist}, Filling->{1->π, 2->π}, PlotStyle->{Blue,Red,Green}, Joined->{False, False, True}]

Just eyeballing it, it looks like the coprimality test gives a better approximation to π (at least it did when I ran this code)