How to use a for loop to solve ODE

Hi I am trying to write a code to simulate the distribution of a drug in someone's body. I am having some trouble creating an iterative process so that the ODE's are solved per a time interval and at each time interval my initial conditions are in a way re-applied. In context I basically want the loop to simulate someone taking another dose of the medication a certain amount of times over the course of time the code is run.

Here is my code

function DrugDelivery
prompt = 'Time between dosages?'
tint = input(prompt);
prompt = 'Total time?'
tfin = input(prompt)
tspan = [0:tint:tfin]
for k = 0:tint
    y0 = [500; 0; 0]
          [t, y] = ode45(@ODEfun, k, y0)
    end         
Co = y(:,1); 
Cc = y(:,2);
Cp = y(:,end);
for i = 1: size(y,3)
    disp(['     Solution for Dependant Variable y'   int2str(i)]);
    disp(['           t            y'   int2str(i)]);
    disp([t  y(:,i)]);
    plot(t,Co);

      title(' Plot of Dependant Variables (Co,Cc,Cp) vs Time ');
      xlabel('  Independant Variable t  ')
      ylabel(  '  Dependant Variable Concentration  ');

      hold on
      plot(t,Cc);
      hold on
      plot(t,Cp);

      legend('Co', 'Cc', 'Cp')

      hold off
  end 

function dydt = ODEfun(t, yfunvec)
Co = yfunvec(1);
Cc = yfunvec(2);
Cp = yfunvec(3);
ka = 0.07;
kr = 0.03;
ke = 0.02;
kc = 0.004;
kp = 0.006;
dCodt = -ka.*Co;
dCcdt = -(ke+kr+kc).*Cc + kp.*Cp + ka.*Co;
dCpdt = -kp.*Cp + kc.*Cc;
dydt = [dCodt; dCcdt; dCpdt];

I know my loop is currently very wrong but if you could please help me understand how to achieve the kind of loop that I described above I would greatly appreciate it!

Answer 1

[yeungor]

There are two options: use a prebuilt ode integrator like ode45, or you can explicitly compute at each step. The latter sounds more like what you intend to do.

Let's say y is a vector (e.g. y = [Co;Cc;Cp]) and Y(:,i) is y at time t=i. Usually we can write something like dy/dt = f(t,y), the ODE, so that we can write something like

y(i+1) = y(i) + dy(i); % inside your for loop will look something like this

There are a few methods for computing dy(i) and in general it depends on y(i) and t(i), a common class of methods are called runge-kutta methods, and ODE45 implements one that is called RK45, a 4th order runge-kutta with 5th order error estimation. But RK4 is very easy to implement. An example is below for the following problem

tyy' + 4t^2+y^2 = 0, t>0 from here the solution is y(t)=\pm\sqrt( (a-2x^4)/(x^2) )

y' = (-4t^2-y^2)/(ty) = f(t,y)

Y = [y0 zeros(length(y0), length(tspan))];
for i=1:length(tspan)
  ti = tspan(i); yi = Y(:,i);
  k1 = f(ti, yi);
  k2 = f(ti+dt/2, yi+dt*k1/2);
  k3 = f(ti+dt/2, yi+dt*k2/2);
  k4 = f(ti+dt  , yi+dt*k3);
  dy = 1/6*(k1+2*k2+2*k3+k4);
  Y(:,i+1) = yi +dy;
end

plot(tspan, Y(1,:), tspan, Y(2,:), tspan, Y(3,:))

The other option is using ode45, it will figure out how big dt on the fly, but it will vary the size of it, so t would be unevenly spaced. Use it this way.

[t,Y] = ode45(@ODEFUN, [0 tfin], y0);
yy(:,1) = interp1(t, Y(:,1), tspan);

The interp1 line will poll the data at all of the times in tspan and it will interpolate between times close to each value in tspan.

Instead of doing the second part, instead of giving a 2 element vector for the second argument you can give a vector of all the times you want to evaluate the function. Note that this means ode45 will use the supplied times and not use adaptive step methods.

[t,Y] = ode45(@ODEFUN, tspan, y0);

Last thing, you have an input query that specifies tint as the time between dosages. I'm not sure about your system, but unless dosages are very frequent, then you may want a smaller number.

[Alison Michell]

 Okay, This is the code using interpolation, but each iteration my tspan restarts from zero, how would I fix that? 

 for k = 0:((tfin./tint)-1);
      x = 0;
      tspan = [x:tint]
      y0 = y0;
      [t, y] = ode45(@ODEfun, tspan, y0);
      for h = 1;
          x = tint
          tint = tint + 120;
          t1 = x
          yt1 = interp1(t,y,t1)
          y0 = yt1 + [500, 0, 0]
      end  
  end

[yeungor]

The purpose of the for loop I mentioned is to compute y(i+1), and ode45 implements this for loop, so you only need the for loop for each time you reset the initial condition (another dose). Below is what I would do. One call of ode45 for each time you dose, just giving start and end times. When it's all finished, interpolate to evenly spaced times and plot.

Y = []; % store the solution for all times over all doses
T = []; % store all times over all doses
tstart = (0:numDoses-1)*timeBetweenDoses % times you reset IC (new dose)
for i=1:numDoses
  tspan = tstart(i) + [0, timeBetweenDoses]; 
  [t, y] = ode45(@ODEfun, tspan, y0); 
  Y = [Y; y];
  T = [T; t];
end

% make evenly spaced time to interpolate to
tt = linspace(0, (numDoses-1)*timeBetweenDoses);
YY = zeros(length(tt), 3); % preallocate
for i=1:size(Y,2) % interpolate each of your concentration variables
  YY(:,i) = interp1(T, Y, tt);
end
plot(tt,YY);

Since this way is resizing the array on the fly (Y=[Y;y];), it's not the fastest, but there probably aren't that many dosages either, so it's not an issue. If there are a lot of dosages, an alternative would be using a cell array to stack each of these vectors on top of each other and cell2mat afterwards will shove it all together.

Y = cell(numDoses,1); % store the solution for all times over all doses
T = cell(numDoses,1); % store all times over all doses
tstart = (0:numDoses-1)*timeBetweenDoses % times you reset IC (new dose)

for i = 1:numDoses
  tspan = tstart(i) + [0, timeBetweenDoses]; 
  [t, y] = ode45(@ODEfun, tspan, y0); 
  Y{i} = y; 
  T{i} = t;
  y0  = y0 + [500,0,0];
end

T = cell2mat(T);
Y = cell2mat(Y);
% then interpolate and plot

One thing to notice in your code is that you set

tspan = [0:tint]; % because x is always 0

And that's why it always started at 0. Note this isn't necessarily a problem, but it is kinda weird to work with, especially since you'll want to interpolate the data.

[Alison Michell]

When I try to store the values I get this error

Error using griddedInterpolant
The grid vectors must contain unique points.
Error in interp1 (line 148)
        F = griddedInterpolant(Xext,V,method);
Error in DrugDelivery (line 33)
      YY(:,i) = interp1(T, Y, tt)

[yeungor]

I do have a bug and I found it. But you should be able to debug on your own. What the error means is that the first argument in interp1 had repeat values, which means it's not a proper function, so it can't interpolate.

Why are there repeat values in T? Even if you don't know why, where would they be and how would you get rid of them? Note that T= [t1;t2;t3...] where ti is time from the ith dose.