Introduction to Mathematica SS 2023.

Solve two of the following three exercises.

Exercise 1. Central limit theorem.

Exercise 1. Verify the Central Limit Theorem. Generate nx random variates distributed uniformly within - 1, 1. Call y the sum of these variables . Generate a large number ny  (eg 10000) of such y variates . Find the mean and variance for the y variates . Show graphically  (eg using Histogram) that for increasing nx  the distribution of y approximates a Gaussian with the appropriate mean and variance. Optional: repeat for different distributions (not Gaussian) of the nx variates.

y={};clt[nx_,ny_]:=Module[{vy},Do[{x=Table[RandomReal[{-1,1}],{i,nx}];AppendTo[y,Total[x]]},{j,ny}];
vy=nx*4/12;Show[Histogram[y,{.5},"PDF"],Plot[1/Sqrt[2*Pi*vy]*Exp[-1/2*x^2/vy],{x,-5*Sqrt[vy],5*Sqrt[vy]}],Frame->True]]
clt[40,10000]

y={};clt[nx_,ny_]:=Module[{vy},Do[{x=RandomVariate[UniformDistribution[{-1,1}],nx];AppendTo[y,Total[x]]},{j,ny}];
vy=nx*4/12;Show[Histogram[y,{.5},"PDF"],Plot[1/Sqrt[2*Pi*vy]*Exp[-1/2*x^2/vy],{x,-5*Sqrt[vy],5*Sqrt[vy]}],Frame->True]]
clt[10,10000]

y={};clt[nx_,ny_]:=Module[{vy,muy},Do[{x=RandomVariate[PoissonDistribution[1],nx];AppendTo[y,Total[x]]},{j,ny}];
vy=Variance[y];muy=Mean[y];Show[Histogram[y,{1},"PDF"],Plot[1/Sqrt[2*Pi*vy]*Exp[-1/2*(x-muy)^2/vy],{x,muy-3*Sqrt[vy],muy+3*Sqrt[vy]}],Frame->True]]
clt[30,10000]

Exercise 2. A Three-body problem.

Exercise 2. Three - body problem. Evolve three masses under their own gravity on a plane. Plot the resulting trajectory (Optional: you might use an Animation or a a Manipulate command to visualize the trajectories). You can use DSolve but also, perhaps more easily, perform a Do in which for every time step you calculate the force on body A from the other two bodies, then the acceleration,  then the velocity, and finally the new positions of A, with little time increments, and repeat for every body. In this way you can  build a list of positions to be plotted. Use the constants below.

(*Define the gravitational constant*)

G=6.67430*10^(-11); (*m^3 kg^(-1) s^(-2)*)

(*Define the masses of the bodies*)
m1=5*10^7; (*Mass of body 1 in kg*)
m2=3*10^7; (*Mass of body 2 in kg*)
m3=4*10^7; (*Mass of body 3 in kg*)

(*Define the initial positions and velocities of the bodies*)
pos1={1,0}; (*Initial position of body 1 in meters*)
pos2={-1,0}; (*Initial position of body 2 in meters*)
pos3={0,1}; (*Initial position of body 3 in meters*)
vel1={0,0}; (*Initial velocity of body 1 in m/s*)
vel2={0,0}; (*Initial velocity of body 2 in m/s*)
vel3={0,0}; (*Initial velocity of body 3 in m/s*)

(*Define the time step and total time*)
dt=0.01; (*Time step in seconds*)
totalTime=200; (*Total simulation time in seconds*)

(*Initialize lists to store the positions of the bodies over time*)
positions1={pos1};
positions2={pos2};
positions3={pos3};

(*Perform the simulation*)
numSteps=Floor[totalTime/dt];
For[i=1,i<=numSteps,i++,
(*Calculate the distances between the bodies*)
r1=Norm[pos2-pos1];
r2=Norm[pos3-pos1];
r3=Norm[pos3-pos2];
(*Calculate the gravitational forces*)force1=(G*m1*m2/r1^2)*Normalize[pos2-pos1]+(G*m1*m3/r2^2)*Normalize[pos3-pos1];
force2=(-G*m1*m2/r1^2)*Normalize[pos2-pos1]+(G*m2*m3/r3^2)*Normalize[pos3-pos2];
force3=(-G*m1*m3/r2^2)*Normalize[pos3-pos1]+(-G*m2*m3/r3^2)*Normalize[pos3-pos2];
(*Update the velocities of the bodies*)
vel1=vel1+force1/m1*dt;
vel2=vel2+force2/m2*dt;
vel3=vel3+force3/m3*dt;
(*Update the positions of the bodies*)
pos1=pos1+vel1*dt;
pos2=pos2+vel2*dt;
pos3=pos3+vel3*dt;
(*Store the positions in the lists*)
AppendTo[positions1,pos1];
AppendTo[positions2,pos2];
AppendTo[positions3,pos3];]

Manipulate[ListPlot[{positions1[[1;;i]],positions2[[1;;i]],positions3[[1;;i]]},PlotRange->{{-2,2},{-2,2}},Frame->True,Axes->False],{i,2,Length[positions1]}]

Exercise 3. MonteCarlo integration.

3. Define the function f (x) as below . Plot it in the interval (0, 3). Integrate f(x) inside this range. Find numerically the maximum f_max of f(x) within the range (0,3). Write a code that generate 10^N random points uniformly distributed  in the rectangular region defined by x between (0,3) and y betwen 0 and f_max and counts how many points are below f(x). Plot the points below f(x) and the function itself for a couple values of N, e.g. N=2 and N=4 or create a Manipulate command showing the plot when varying N. The ratio of points below f(x) and the total number of random points 10^N, multiplied by the total area of the rectangular region, should approximate the exact integral of f(x) between 0 and 3 (MonteCarlo integration). Plot the result of the MonteCarlo integration as a function of N and show graphically that it converges to the  true value for N increasing from 1 to 5 or more.

Clear[f];f[x_]=Sin[x]*Exp[-x^2];

Clear[x,fd];fd[x_]=D[f[x],x];root=FindRoot[fd[x]==0,{x,1}][[1,2]];f[root]

f=Compile[{{x,_Real}},Sin[x]*Exp[-x^2]];
Plot[f[x],{x,0,3},Frame->True]

fint[y_]=Integrate[f[x],{x,0,y}]

makepoints[n_]:=Module[{ic=0,out=0,x,y,points={}},Do[{y=RandomReal[{0,0.5}];x=RandomReal[{0,3}];If[y<=f[x],{++ic,AppendTo[points,{x,y}]},++out]},{i,10^n}];counts={ic,out};points]

Manipulate[ListPlot[makepoints[n],Frame->True,PlotRange->{{0,3},{0,.5}}],{n,1,5,1}]

list={};Do[{makepoints[n],AppendTo[list,{n,counts[[1]]/(counts[[1]]+counts[[2]])*(3*0.5)}]},{n,1,5.5,.25}];list
Show[ListPlot[list,Frame->True],Plot[fint[3.],{n,0,5.5}]]