How to code a pulse function in differential equations
10 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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
0 commentaires
Réponses (1)
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
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.
Voir également
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!