MATLAB Answers

How to code a pulse function in differential equations

11 views (last 30 days)
Rui
Rui on 10 Jan 2017
Commented: Star Strider on 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

Answers (1)

Star Strider
Star Strider on 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 Comments
Star Strider
Star Strider on 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.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!