How to code a pulse function in differential equations

10 次查看(过去 30 天)
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
  1 个评论
태훈
태훈 2024-4-25
If you can't get explicit solutions, I guess you could try "a series of odes" using a for loop. For example,...
tv=[0:1:100]
y0=1
for k=1:length(t)-1
tspan=[tv(k), tv(k+1)];
[t,y]=ode45(@yourfun, tspan...)
y0=y(length(y))+1 % e.g., adding 1 at every 1 hour
end

请先登录,再进行评论。

回答(1 个)

Ayush Aniket
Ayush Aniket 2024-5-9
Hi Rui,
Assuming 'tau' to be a small interval around the hour marks, you can define a custom 'odefun' to pass as argument to the 'ode45' MATLAB function. This custom function should add the mass 'M' if the evalutaion time lies in the 'tau' interval around the hour marks as shown below:
function dXdt = modelODE(t, X, k, V, M, tau, additionTimes)
% Check if t is within any of the mass addition intervals
isAddingMass = any(abs(t - additionTimes) <= tau/2);
% Compute dX/dt
if isAddingMass
dXdt = -k*X/V + M/tau;
else
dXdt = -k*X/V;
end
end
The 'additionTimes' argument specifies the hour marks around which the mass should be added to the model. For your problem it can be defined as below:
additionTimes = 1:23; % Times to add mass
% Time span
tSpan = [0 24];
% Solve the ODE
[t, X] = ode45(@(t, X) modelODE(t, X, k, V, M, tau, additionTimes), tSpan, X0);

类别

Help CenterFile Exchange 中查找有关 Ordinary Differential Equations 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by