Bateman equation using dsolve

I am attempting to use the bateman equation to follow a decay series. When I have a non-zero initial concentration in isotope 2 and/or 3 the program fails to function. It also appears to "instantly" decay isotope 1 when 2 and 3 are zero. Any ideas what I did wrong with applying dsolve and diff(eqn) in q1:q3?

Code Below::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

clear all;
clc;

%start the symbolic solver variables
syms N1(t) N2(t) N3(t) l1 l2 l3 l4 l5 N0 N20 N30

%Constants & initial values
tmax =40;
N0=1e6;
N20=0;              
N30=0;              
% *I should be able to change these values to the IBV problem and still get a solution*

tmax = 40; 
%set up the differential equations
*I should be able to solve each individually.  See Wikipedia <https://en.wikipedia.org/wiki/Bateman_Equation bateman equation> *

q1=dsolve(diff(N1)==-l1*N1, N1(0) == N0);   
q2=dsolve(diff(N2)==-l2*N2+l1*q1, N2(0) == N20);
q3=dsolve(diff(N3)==l2*q2, N3(0)==N30);         
%the time interval is so short this does not start to decay- simplification to calculation

 *%below this line everything works as intended* 

t1 = (1.2e-4);   
t2 = (37);       
t3 = (12*24*60*60);  

% change to lamba values
k1=log(2)/t1;
k2=log(2)/t2;
k3=0;
N0val = 1;  % value for N1(0) 

% Substitutions
eq1=subs(q1, [l1,l2,l3, N0], [k1,k2,k3, N0val]);
eq2=subs(q2, [l1,l2,l3, N0], [k1,k2,k3, N0val]);
eq3=subs(q3, [l1,l2,l3, N0], [k1,k2,k3, N0val]);

% Set graph limits
tmin=0;
ymin=0;
ymax=1;

% Plot easily via ezplot
figure
hold on;
h1=ezplot(char(eq1), [tmin,tmax,ymin,ymax]);
h2=ezplot(char(eq2), [tmin,tmax,ymin,ymax]);
h3=ezplot(char(eq3), [tmin,tmax,ymin,ymax]);

title('Isotope fraction')
legend('N1', 'N2','N3')
xlabel('t in seconds')
set(gca, 'XScale', 'log');

Answer 1

[Star Strider]

Formatting your code would have helped for a start. See: TUTORIAL: how to format your question with markup.

I’m not familiar with the Bateman equations, so doing a bit of research (see Introduction to Nuclear Science - GSI, slides 23-4), it seems to me that these equations:

h1=ezplot(char(eq1), [tmin,tmax,ymin,ymax]);
h2=ezplot(char(eq2), [tmin,tmax,ymin,ymax]);
h3=ezplot(char(eq3), [tmin,tmax,ymin,ymax]);

need to be solved and then summed rather than plotted individually. (It doesn’t appear that they be solved as a system.) See if that works for you.

I’ve not done anything with radioactive decay since undergraduate physical chemistry, so a relevant free reference on the Bateman equations would help.

[Derek Haws]

Thank you for taking the time to reply and helping me make the question better. h1:h3 are the final solutions to the ode when the variables have been applied in eq1:eq3 and do exactly what they are suppose to. The problem is with q1:q3 when the IBV problem has a non-zero condition for the second and third isotope.

Everything beyond line 15 works as intended.

[Star Strider]

I need to know more about the Bateman equation, how you’re using it, and what you want to get from it.

[Derek Haws]

The Bateman equation is just an exponential decay ODE where N`(t) = -N*Lambda whose solution is just N = N(0) exp (-lambda*t).* Lambda is the decay constant for that particular isotope. N(0) is the initial number of the first isotope/particle. When the isotope decays into a second it obeys the rate law adding the decays from the first to the concentration of the second while allowing the second to decay. The rate law is N`(t)=N1*Lambda1-N2*lambd2.

The code as written will solve for the initial isotope decaying into the second and the second decaying into the third. But only if the initial concentration of 2 and 3 is zero. I need to be able to solve the ode with the initial values of 2 and/or 3 non-zero. Currently the code returns zeros for all three isotopes as soon as I change the initial values for 2 and/or 3 to a non-zero entry. I left the equations as an ode and not just a numerical solution hoping that it would provide greater flexibility and since the solution is well known and explicit I wasn't concerned with numerical integration.

[Star Strider]

I’m familiar with first-degree ODEs describing radioactive decay, so that’s not the problem. In my recent research, there were a number of different implementations of the Bateman equation with respect to daughter particles and such. I don’t know what your code represented.

How do you want to define ‘N20’ and ‘N30’? Right now, you’re defining them as zero when you solve the original ODEs, and that seems to be the problem. If you want to define them as different isotope amounts depending on the previous decay equation ‘q1’ at a specific time, you have to specify that.

I didn’t get the impression that the Bateman equations were solved as a system of simultaneous differential equations (I don’t have that much experience with them), but that could be one approach. Define ‘N10’ and let everything else solve itself.

[Derek Haws]

I have tried to define N20 and N30 as non-zero entries. When they are zero, as the code currently shows, it works. When they are non-zero it fails to solve. It will not even solve the first ODE which has no dependence on the second isotope concentration. The code should be able to solve the decay for the first and then add that solution to the solution for the second, regardless of the starting amount for number 2. The reason I posted was to find help to this problem.

Without branching, each parent daughter decay is a separate problem which is why I don't understand the problem with the code as written. It handles the zero entries for the second isotopes but nothing else.

[Derek Haws]

Changing it to something like this doesn't work at all:

• N0=1e6;
• N20=4e5;
• N30=2e3;
• q1=dsolve(diff(N1)==-l1*N1, N1(0) == N0);
• q2=dsolve(diff(N2)==-l2*N2+l1*N1, N2(0) == N20, N1(0) ==N0);
• q3=dsolve(diff(N3)==l2*N2-l3*N3, N2(0)==N20, N3(0)==N30);

[Star Strider]

I’m thinking that this is similar to a compartmental analysis problem and a system of differential equations describing it.

What process are you simulating? What are the three isotopes, and how are they related to the decay process? (What becomes what else?)

[Derek Haws]

number 1 decays into two ... decays into three ... that's all ... actual isotope is not a concern. It should be able to do any set of initial concentrations with a given half-life or decay constant.

[Star Strider]

It can. You have a system of first-order ODEs, so you have to solve them as such. That was the principal problem.

I had problems getting ezplot to display correctly, so I used the matlabFunction function to convert the solved equations to anonymous functions and plotted them conventionally using subplot in figure(2).

I made some changes (I kept your original code but commented-out the parts I changed), and I got it to work. From my experience with such systems, it gives appropriate results:

syms N1(t) N2(t) N3(t) l1 l2 l3 l4 l5 N0 N20 N30

%Constants & initial values
tmax =40;
N0=1e6;
N20=4e5;
N30=2e3;             
% *I should be able to change these values to the IBV problem and still get a solution*

tmax = 40; 
%set up the differential equations
% *I should be able to solve each individually.  See Wikipedia <https://en.wikipedia.org/wiki/Bateman_Equation bateman equation> *

% q1=dsolve(diff(N1)==-l1*N1, N1(0) == N0);   
% q2=dsolve(diff(N2)==-l2*N2+l1*q1, N2(0) == N20);
% q3=dsolve(diff(N3)==l2*q2, N3(0)==N30); 

q = dsolve(diff(N1)==-l1*N1, diff(N2)==-l2*N2+l1*N1, diff(N3)==l2*N2, N1(0) == N0, N2(0) == N20,  N3(0)==N30);

t1 = (1.2e-4);   
t2 = (37);       
t3 = (12*24*60*60);  

% change to lamba values
k1=log(2)/t1;
k2=log(2)/t2;
k3=0;
N0val = 1;  % value for N1(0) 

% Substitutions
% eq1=subs(q.N1, [l1,l2,l3, N0], [k1,k2,k3, N0val]);
% eq2=subs(q.N2, [l1,l2,l3, N0], [k1,k2,k3, N0val]);
% eq3=subs(q.N3, [l1,l2,l3, N0], [k1,k2,k3, N0val]);
eq1=subs(q.N1, [l1,l2,l3], [k1,k2,k3]);
eq2=subs(q.N2, [l1,l2,l3], [k1,k2,k3]);
eq3=subs(q.N3, [l1,l2,l3], [k1,k2,k3]);

% Set graph limits
tmin=0;
tmax = 0.5;
ymin=0;
ymax=1;

% Plot easily via ezplot
figure(1)
hold on;
% h1=ezplot((eq1), [tmin,tmax,ymin,ymax]);
% h2=ezplot((eq2), [tmin,tmax,ymin,ymax]);
% h3=ezplot((eq3), [tmin,tmax,ymin,ymax]);
h1=ezplot((eq1), [tmin,tmax]);
h2=ezplot((eq2), [tmin,tmax]);
h3=ezplot((eq3), [tmin,tmax]);
hold off
title('Isotope fraction')
legend('N1', 'N2','N3')
xlabel('t in seconds')
set(gca, 'XScale', 'log', 'YScale','log');

eq1fcn = matlabFunction(eq1);
eq2fcn = matlabFunction(eq2);
eq3fcn = matlabFunction(eq3);

tv = linspace(1E-4, 1E-1, 1E+5);

figure(2)
subplot(3,1,1)
plot(tv, eq1fcn(tv));
title('Isotope fraction N_1')
set(gca, 'XScale', 'log', 'YScale','log');
grid

subplot(3,1,2)
plot(tv, eq2fcn(tv));
title('Isotope fraction N_2')
set(gca, 'XScale', 'log', 'YScale','log');
grid

subplot(3,1,3)
plot(tv, eq3fcn(tv));
title('Isotope fraction N_3')
xlabel('t in seconds')
set(gca, 'XScale', 'log', 'YScale','log');
grid

[Derek Haws]

Thank you, I was working with my advisor on this today. Trouble shooting the code showed that the code was solving correctly and that the problem was in the ezplot. We also used the matlabFunciton and plot to see that it was doing the calculations correctly. Sorry I was not able to pose the question correctly. I appreciate your help. Your code will also solve the problem and produces the same results.

[Star Strider]

My pleasure.

Once I understood your system you’re studying, I knew the problem was that you needed to integrate it as a system, and not as individual equations.

The easiest way would be to create one anonymous function from the equations in your system, and then use ode15s to solve it. (It’s a ‘stiff’ system because the magnitudes of the parameters vary across a wide range.) The anonymous function will pick up the values of the parameters from your workspace, so you could probably integrate and plot your equations in at most a couple dozen lines of code.