How should I set the value of u_tau?
6 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Hello everyone, I would like to use pdepe for solving a SIR partial differential equations in one dimension with time lag, so it looks good. But there has been a problem setting the value of u_tau, resulting in a very different result than expected. I don't know what's wrong with the code.Here is the part of the code involved:
function u_tau = interpolate_history(history, x, t, tau)
% Interpolating function to get the state at time t - tau
t_target = t - tau;
if t_target < history.t(1)
% u_tau = squeeze(history.u(1, :, :)); % Out of history range, use first value
u_tau = 0;
else
% u_tau = interp1(history.t, squeeze(history.u(:, :, :)), t_target, 'linear', 'extrap');
u_tau = 0;
end
end
function u0 = SIRInitialConditions3(x)
S0 = 1.25;
I0 = 0.85;
R0 = 0.74;
u0 = [S0; I0; R0];
end
% This is the part of the main function that deals with u_tau
function [c,f,s] = SIRPDE3(x, t, u, DuDx)
global history;
S = u(1);
I = u(2);
R = u(3);
u_tau = interpolate_history(history, x, t, tau);
S_pre = u_tau(1);
end
% This is the part of the main function that deals with u_tau
function run()
xspan = linspace(0, 10, 100);
tspan = [0, 10, 50];
% Initialising the history
history = initialize_history(xspan, tspan, @SIRInitialConditions3);
sol = pdepe(0, @SIRPDE3, @SIRInitialConditions3, @SIRBoundaryConditions3, xspan, tspan);
I originally ran the simulation with the two equations commented above neither of which resulted in a trend over time, so I set all u_tau to 0. The trend over time then turned out to be correct, but the initial values of I and R were still 0, which was not the expected result. I think there might be something wrong with the part function u_tau = interpolate_history(history, x, t, tau), how should I change it?
Thanks to the community. :)
6 comentarios
Bill Greene
el 13 de Jun. de 2024
Editada: Bill Greene
el 13 de Jun. de 2024
Apparently you are still not understanding my questions. I will try one more time.
Your initialize_history function is setting history.u to an array of size number_of_times x number_of_x_values x 3. In your commented lines in interpolate_history, u_tau will be calculated as a number_of_x_values x 3 array. NOT a scalar as you show in your un-commented line
u_tau=0;
Your SIRPDE3 function contains the line
S_pre = u_tau(1);
This is simply setting S_pre to the value of S at the left end of your domain. It is hard for me to believe that this is actually what you intended?
It seems more likely to me that either your formulation or implementation is incorrect. So, if you really want help with this issue, I STRONGLY suggest, yet again, that you either provide a link to a reference describing the formulation or provide a detailed mathematical description (shown as equations, not matlab code).
Respuestas (1)
Torsten
el 8 de Jun. de 2024
Editada: Torsten
el 8 de Jun. de 2024
The problem simply is that you cannot solve delay PDEs with "pdepe", at least not directly.
A possible remedy:
Call "pdepe" in a loop from 0 to tau, from tau to 2*tau and so on. Each time control is returned to the calling program, change the history data to the results you obtained in the preceeding step of length tau and make them accessible where needed. Use "interp1" to interpolate to the time t - tau.
10 comentarios
Torsten
el 17 de Jun. de 2024
Editada: Torsten
el 17 de Jun. de 2024
After the call to "pdepe", add the lines
TSPAN((i-1)*nt+1:i*nt) = tspan;
SOL((i-1)*nt+1:i*nt,:,:) = sol;
Then xspan, TSPAN and SOL contain all information you need for your surface plots.
Since S, I and R are constant over x for all times t, a surface plot does not add more information to the usual plot I included. You shouldn't always do what your teachers tell you to do if you have arguments that support your point-of-view.
Ver también
Categorías
Más información sobre Eigenvalue Problems en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!