How to use Eulers Algorithm to solve half life decay on matlab

clc; clear all;
n = 100; % Number of iterations / decays
N = 1:n+1; % Array containing 100 elements
N0=1000; % Initial number of atoms
N(1)= N0; % First element in array N or initial number of atoms
time(1)=0;
tau = 1.0;
Delta = 0.02;
for i=1:n+1 % Compute the remaining atoms after each time
 N(i+1)=N(i)-Delta*N(i)/tau % Euler algorithm for decay process.
 time(i+1)=i*Delta
end
ms = 1.5;
figure(1);
plot(time,N,'-ob','Markersize',ms);xlabel('Time(s)'); ylabel('Number of atoms');

Im trying to use the algorithm above to write a Matlab routine to determine the half-life of the radioactive element and compare it against the expected/theoretical value

[DGM]

What exactly is the goal? You do seem to be generating an exponential decay curve. What would be the point of calculating lambda when you already knew tau to start with? Are you supposed to be pretending that the generated N represents measured data with an unknown time constant?

For finding the time constant, are you intended to use standard Matlab tools, or are you supposed to somehow rewrite the above code to do the task?

[Rookiee]

Thanks @DGM Im required to write a code to derive the time constant and the half life. Thanks again

Answer 1

[DGM]

I posted a comment already, but I'll make some assumptions and throw down an attempt at an answer.

Let's assume that you're only using this code to generate vectors that are to be treated as a stand-in for real-world measured data. I'm going to assume that this code is done and we don't need to do anything with it or invent any convoluted means of using it to back-calculate tau or lambda. At this point, we just have two vectors (time, N) in the workspace and we need to find tau or lambda.

If we weren't living in homework-land where everything tends to be done the hard way, we might make use of the tools provided by Matlab.

thisfit=fit(time',N','exp1'); % do a single-term exponential curve fit
N0calc=thisfit.a % this is the initial value
taucalc=-1/thisfit.b; % this is the back-calculated tau
N2=N0calc*exp(-1/taucalc*time); % you can plot this against your original curve

If you wanted a more simplistic approach, remember that we already know the form of the equation:

We know that at t0, N=N0. At a given time t1:

We can just solve for tau:

... and of course, lambda is 1/tau.

t1=time(end); % pick some time
tauest=-t1/log(N(end)/N0) % calculate tau
lambdaest=1/tauest % convert to lambda if desired (it's the same for 1, obviously)

which would give us an estimate of tauest=0.98997

EDIT: I forgot you actually wanted the half-life. That's easy enough.

halflife=tauest*log(2)

[Rookiee]

@DGM thanks for your assistance, please how do I calculate the time constant?

[DGM]

tau is the time constant

lambda is the decay rate

[Rookiee]

@DGM I get this error when I try to run the code

Check for missing argument or incorrect argument data type in call to function 'fit'.

[DGM]

You probably don't have curve fitting toolbox. There are multiple versions of the function fit(), and it's probably using the one from the stats toolbox, which doesn't do the same thing.

which fit -all

If you still want to use some curve fitting tools, try this:

% since we don't have fit() for doing exponential fits
% linearize by taking the natural log of the data
Nlog = log(N);
% now we can do a first-order polynomial fit
c = polyfit(time,Nlog,1);
N0calc = exp(c(2)); taucalc = -1/c(1);
N2 = N0calc*exp(-1/taucalc*time);
plot(time,N,'o',time,N2)