How to use a for loop to solve ODE

37 次查看(过去 30 天)
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!

回答(1 个)

yeungor
yeungor 2018-7-23
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.
  4 个评论
Alison Michell
Alison Michell 2018-7-23
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
yeungor 2018-7-24
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.

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Loops and Conditional Statements 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by