How to code a pulse function in differential equations

10 vues (au cours des 30 derniers jours)
Rui
Rui le 10 Jan 2017
Commenté : Star Strider le 11 Jan 2017
Hello everyone,
could you help me with my problem described below?
assuming I have a simple model dX/dt = -k*X/V, where X is amount, t is time (between 0 and 24 hours), V is volume, and k is rate.
I would like to add a certain amount of mass (M) into my system every 1 hour, in other words, adding mass on 1, 2, 3...23 hour. It does not need to be exactly on every hour, but could be a small interval (tau) around 1, 2, 3, ... so equation becomes something like
dX/dt = -k*X/V + M/tau, M is zero when t is out of these small intervals.
I just cannot figure out a smart way to code this, and then solve the ODE. Could you give me some idea? Thanks in advance!
Rui

Réponses (1)

Star Strider
Star Strider le 10 Jan 2017
I would solve the ODE first, since it has a straightforward analytic solution, then ponder the coding:
syms X(t) k V M tau X0
% % dX/dt = -k*X/V + M/tau
Dx = diff(X);
Eqn = Dx == -k*X/V + M/tau;
X = dsolve(Eqn, X(0) == X0);
Xfcn = matlabFunction(X)
Xfcn = @(M,V,X0,k,t,tau) (M.*V-exp(-(k.*t)./V).*(M.*V-X0.*k.*tau))./(k.*tau);
The ‘X0’ value is the initial condition for ‘X(t)’, that is ‘X(0)’.
I would approach it as something similar to a pharmacokinetics problem, and use ‘X0’ as ‘M/tau’ for each interval:
Xfcn = @(V,X0,k,t) X0.*exp(-(k.*t)./V);
You then only need to determine when ‘M/tau’ is added.
Note I’m not certain exactly what you want to do, so much of this is a guess.
  2 commentaires
Rui
Rui le 11 Jan 2017
Hi Star, Thank you very much! I am sorry, I should make it clear.
I am actually trying to put this into a physiological pharmacokinetic model. Every a few hours, I could give my model an oral bolus dosing. The model probably has 20 to 30 ODEs. As such, may not be easy to get an analytical solution.
One thing confusing me is how to code it so that every a few hours, I could dose the system. I feel it is probably some if loop, but cannot think of an easy way to write the code.
I used to use the end of first simulation plus dosing amount as initial condition of second simulation, and then end of second simulation as initial values for the third simulation, and so on. But I am wondering if there's better way to do it.
Star Strider
Star Strider le 11 Jan 2017
If you want to do a lot of pharmacokinetic modeling, consider getting SimBiology. I don’t have it because I no longer do this sort of modeling on a large scale, so I have no experience with it.
You can probably simulate your pharmacokinetic model as a linear control system, with ‘A,B,C,D’ matrices.
This is the best I can do as a prototype of what I believe you want to do:
t = linspace(0, 1, 100)*5;
A = [0 0 1; 0 0 1; -2 -3 -5]; % Compartment Connection Matrix
L = eig(A);
B = [1 0 0];
C = eye(size(A));
D = 0;
x0 = [1; 1; 1]; % Initial Conditions
u = zeros(size(B,2),length(t));
dose = randi(90, 1, 10); % Dose Times Index Vector
for k1 = 1:length(t) % Create Input Matrix
if ismember(k1,dose)
u(:,k1:k1+9) = u(:,k1:k1+9) + [ones(1,10); zeros(2,10)];
end
end
for k1 = 1:length(t)
x = (C*expm(A*t(k1))*x0*B+D).*u(:,k1); % Simulate System
xv(:,k1) = x(:,1);
end
figure(1)
plot(t, xv)
grid
Experiment with it to see how it works, and substitute your matrix of compartment constants as the ‘A’ matrix, with appropriate changes in the sizes of the other matrices.

Connectez-vous pour commenter.

Catégories

En savoir plus sur General Applications dans Help Center et File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by