Chapter 1 
A First Numerical Problem

Radioactive Decay

Sample code 1.1.1. Use the command Plot to plot the curve Nu(t)=N0 e^(-t/tau)

N0:=10000;
(* Initial number of nuclei present *)
tau=10^9;(*half life of U-235*)
Nu[t_]:=N0 Exp[-t/tau];
Plot[Nu[t],{t,0,10tau}]

Sample code 1.1.2. Use the Plot command to plot the function f(x) = e^(-x);

Plot[ Exp[-x],{x,0,10},PlotRange->All]

Sample code 1.1.3. Use the command Solve to solve the equation e^(-t/tau)=1/2;

Solve[Exp[-t/tau]==1/2,t]

Sample code 1.1.4. Use the command Solve to solve the equation e^(-t/tau)=frac, where frac = 0.35;

Solve[Exp[-t/tau]==0.35,t]

Sample code 1.1.5. Extract the numerical value of the solution and print it out. The code in 1.1.3 is improved by defining the solution tprime = t /. Solve[ . .. ] [[1]].

frac=0.35;
Print[Solve[Exp[-t/tau]==frac,t]];
tprime=t/.Solve[Exp[-t/tau]==frac,t][[1]];
Print[tprime]

Sample code 1.1.6. New command DSolve. Use it to ANALYTICALLY solve the first order differential equation dN/dt = -N/tau.

DSolve[Nu'[t]==-Nu[t]/τ,Nu[t],t]

Sample code 1.1.7. New command NDSolve. It is used to NUMERICALLY solve the FIRST order differential equation dN/dt = -N/tau. The numerical values of range of t to plot has to be specified, as in {t, 0, 10tau}. Note also ONE initial condition must be supplied, as in Nu[0]=N0. The command Evaluate[Nu[t] /. Nsol] is used to plot the solution of NDSolve.

tau=10^9;N0=10000;
Nsol=NDSolve[{Nu'[t]==-Nu[t]/tau,Nu[0]==N0},Nu,{t,0,10tau}];
Print[Nsol];
Plot[Evaluate[Nu[t]/.Nsol],{t,0,10tau},PlotRange->All]

Sample code 1.1.8. A For-loop to evulate the equation dN/dt = -N/tau numerically using Euler’s method.

(* Begin of initialisation*)
tau=10^9;Deltat=tau/100;tfinal=10tau;Nu0=10000;
(*End of Initialisation*)

(*Begin calculation*)
NStep=IntegerPart[tfinal/Deltat];Nu[0]=Nu0;
For[i=1,i≤NStep,i++,
Nu[i]=N[Nu[i-1]-(Nu[i-1]/tau)*Deltat];
];(*end for i*)
(*End of calculation*)

(*Output*)
toutputi=Table[Nu[i],{i,0,NStep}];(*end for i*)
ListPlot[toutputt,PlotRange->All]
Print[toutputi];

toutputt=Table[{i*Deltat,Nu[i]},{i,0,NStep}];(*end for i*)
ListPlot[toutputt,PlotRange->All]
Print[toutputt];

(*End of output*)