11 views (last 30 days)

Show older comments

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

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.

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.

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

Start Hunting!