How should I set the value of u_tau?

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. :)

[Bill Greene]

Why don't you define the PDE system you are trying to solve?

[鑫彤 纪]

Sorry, I forgot to write up the PDE equation part because I only attached part of the code.

Thanks for the heads up!

    dSdt = Lambda - gamma * S_pre * I - mu * S;
    dIdt = theta * gamma * S_pre * I - mu * I - beta*I/(1+alpha*I);
    dRdt = (1 - theta) * gamma * S_pre * I - mu * R + beta*I/(1+alpha*I);

    c = [1; 1; 1];
    f = [d1; d2; d3] .* DuDx;
    s = [dSdt; dIdt; dRdt];

[Bill Greene]

I was actually expecting a mathematical description of the system you are solving. Do you have a reference for this PDE-form of the SIR equations you are trying to solve?

In particular, your description of "history" is incomplete. Exactly how do you initialize it and then update it as the solution progresses (this is not straightforward to do in pdepe)? It looks to me like history.u, might contain the complete solution at various previous time points. If you interpolate this at some particular time, you get a vector of length=3*xspan. But you are not referencing u_tau that way later.

[鑫彤 纪]

Thank you for responding！

Here is the code for history：

function history = initialize_history(xspan, tau, initial_conditions)
    % Initialising the history
    history.x = xspan;
    history.t = linspace(-tau, 0, 5)'; % Initialise from -tau to 0 and make sure it's a column vector.
    history.u = zeros(length(history.t), length(xspan), 3);
    for i = 1:length(xspan)
        initial_u = initial_conditions(xspan(i));
        for j = 1:length(history.t)
            history.u(j, i, :) = initial_u;
        end
    end
end

But I don't think this place should be a problem.

[Bill Greene]

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).

[鑫彤 纪]

Thank you for your reply, I am very sorry for my lack of presentation and for my failure to understand what you meant!

Here is the link to the literature I referenced, which contains the partial differential equation I am trying to implement：

Rumor propagation dynamics in social networks under the influence of time delay and diffusion

This reproduction has been bugging me for a while, here is the code after I changed it:

function run3()
    global history;

    xspan = linspace(0, 10, 100);
    tau = 2.6;
    total_time = 1000;
    num_steps = ceil(total_time / tau);

    % Initialising the history
    history = initialize_history(xspan, tau, @SIRInitialConditions3);

    for step = 1:num_steps
        % Define the refinement time step within each time period
        tspan = linspace((step-1)*tau, step*tau, 5);
        sol = pdepe(0, @SIRPDE3, @SIRInitialConditions3, @SIRBoundaryConditions3, xspan, tspan);
        update_history(sol, tspan);
    end

    t = linspace(0, total_time, size(history.u, 1));
    plot_results(t, xspan, history);
end

function [pl, ql, pr, qr] = SIRBoundaryConditions3(xl, ul, xr, ur, t)
    pl = [0; 0; 0];
    ql = [1; 1; 1];
    pr = [0; 0; 0];
    qr = [1; 1; 1];
end

function u0 = SIRInitialConditions3(x)
    S0 = 1.25;
    I0 = 0.85;
    R0 = 0.74;
    u0 = [S0; I0; R0];
end

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
    else
        u_tau = interp1(history.t, squeeze(history.u(:, :, :)), t_target, 'linear', 'extrap');
    end
end

function history = initialize_history(xspan, tau, initial_conditions)
    % Initialising the history
    history.x = xspan;
    history.t = linspace(-tau, 0, 5)'; % Initialise from -tau to 0 and make sure it's a column vector.
    history.u = zeros(length(history.t), length(xspan), 3);
    for i = 1:length(xspan)
        initial_u = initial_conditions(xspan(i));
        for j = 1:length(history.t)
            history.u(j, i, :) = initial_u;
        end
    end
end

function update_history(sol, tspan)
    global history;
    % Get the solution for the last time step
    new_u = sol(end, :, :);
    % Add new solution to history
    history.t = [history.t; tspan(2:end)'];
    history.u = cat(1, history.u, sol(2:end, :, :));
    % Ensure that the length of the history does not exceed tau
    max_len = length(history.x) * 5;
    if size(history.u, 1) > max_len
        history.t = history.t(end-max_len+1:end);
        history.u = history.u(end-max_len+1:end, :, :);
    end
end

function plot_results(t, xspan, history)
    S = squeeze(history.u(:, :, 1));
    I = squeeze(history.u(:, :, 2));
    R = squeeze(history.u(:, :, 3));

    figure;
    subplot(1, 3, 1);
    surf(xspan, t, S);
    title('Susceptible (S)');
    xlabel('Distance (x)');
    ylabel('Time (t)');   
    shading interp;

    subplot(1, 3, 2);
    surf(xspan, t, I);
    title('Infectious (I)');
    xlabel('Distance (x)');
    ylabel('Time (t)'); 
    shading interp;

    subplot(1, 3, 3);
    surf(xspan, t, R);
    title('Recovered (R)');
    xlabel('Distance (x)');
    ylabel('Time (t)'); 
    shading interp;
end

function [c,f,s] = SIRPDE3(x, t, u, DuDx)
    global history;
    
    Lambda = 0.8; 
    gamma = 0.6; 
    mu = 0.3;
    alpha = 0.5;
    beta = 0.2;
    theta = 0.7;
    tau = 2.65;
    d1 = 0.1; d2 = 0.1; d3 = 0.1;

    S = u(1);
    I = u(2);
    R = u(3);
    u_tau = interpolate_history(history, x, t, tau);
    S_pre = u_tau(1);

    dSdt = Lambda - gamma * S_pre * I - mu * S;
    dIdt = theta * gamma * S_pre * I - mu * I - beta*I/(1+alpha*I);
    dRdt = (1 - theta) * gamma * S_pre * I - mu * R + beta*I/(1+alpha*I);

    c = [1; 1; 1];
    f = [d1; d2; d3] .* DuDx;
    s = [dSdt; dIdt; dRdt];
end

Can you help me see what the problem is? How should I modify it?

Answer 1

[Torsten]

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.

[鑫彤 纪]

Thank you very much for your reply and help! I will try the suggestions you provided right away, and thank you again for your enthusiastic answers!

[鑫彤 纪]

Thanks again for your reply, I added a loop to it here, but for some reason it didn't come out quite the same as expected.

Here is the expected result:

This is the actual result:

Here is part of the code after the change:

function history = initialize_history(xspan, tau, initial_conditions)
    % Initialising the history
    history.x = xspan;
    history.t = linspace(-tau, 0, 5)'; % Initialise from -tau to 0 and make sure it's a column vector.
    history.u = zeros(length(history.t), length(xspan), 3);
    for i = 1:length(xspan)
        initial_u = initial_conditions(xspan(i));
        for j = 1:length(history.t)
            history.u(j, i, :) = initial_u;
        end
    end
end

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
    else
        u_tau = interp1(history.t, squeeze(history.u(:, :, :)), t_target, 'linear', 'extrap');
    end
end

function run3()
    global history;
    xspan = linspace(0, 10, 100);
    tau = 2.6;
    total_time = 1000;
    num_steps = ceil(total_time / tau);

    % Initialising the history
    history = initialize_history(xspan, tau, @SIRInitialConditions3);

    for step = 1:num_steps
        % Define the refinement time step within each time period
        tspan = linspace((step-1)*tau, step*tau, 5);
        sol = pdepe(0, @SIRPDE3, @SIRInitialConditions3, @SIRBoundaryConditions3, xspan, tspan);
        update_history(sol, tspan);
    end

    % Plotting results
    t = linspace(0, total_time, size(history.u, 1));
    plot_results(t, xspan, history);
end

I'm not sure what's wrong, can you fix it for me?

[Torsten]

I don't know what you do in the function "update_history".

It would be easiest if you included executable code to reproduce the results together with a mathematical description of your problem. Further I don't think that 5 interpolation points will generally be enough.

[鑫彤 纪]

Thank you for your reply!

I thought I couldn't paste the full code before so I only posted part of it up.

Here is the full code after I added the loop:

function run3()
    global history;
    xspan = linspace(0, 10, 100);
    tau = 2.6;
    total_time = 1000;
    num_steps = ceil(total_time / tau);
    % Initialising the history
    history = initialize_history(xspan, tau, @SIRInitialConditions3);
    for step = 1:num_steps
        % Define the refinement time step within each time period
        tspan = linspace((step-1)*tau, step*tau, 5);
        sol = pdepe(0, @SIRPDE3, @SIRInitialConditions3, @SIRBoundaryConditions3, xspan, tspan);
        update_history(sol, tspan);
    end
    t = linspace(0, total_time, size(history.u, 1));
    plot_results(t, xspan, history);
end
function [pl, ql, pr, qr] = SIRBoundaryConditions3(xl, ul, xr, ur, t)
    pl = [0; 0; 0];
    ql = [1; 1; 1];
    pr = [0; 0; 0];
    qr = [1; 1; 1];
end
function u0 = SIRInitialConditions3(x)
    S0 = 1.25;
    I0 = 0.85;
    R0 = 0.74;
    u0 = [S0; I0; R0];
end
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
    else
        u_tau = interp1(history.t, squeeze(history.u(:, :, :)), t_target, 'linear', 'extrap');
    end
end
function history = initialize_history(xspan, tau, initial_conditions)
    % Initialising the history
    history.x = xspan;
    history.t = linspace(-tau, 0, 5)'; % Initialise from -tau to 0 and make sure it's a column vector.
    history.u = zeros(length(history.t), length(xspan), 3);
    for i = 1:length(xspan)
        initial_u = initial_conditions(xspan(i));
        for j = 1:length(history.t)
            history.u(j, i, :) = initial_u;
        end
    end
end
function update_history(sol, tspan)
    global history;
    % Get the solution for the last time step
    new_u = sol(end, :, :);
    % Add new solution to history
    history.t = [history.t; tspan(2:end)'];
    history.u = cat(1, history.u, sol(2:end, :, :));
    % Ensure that the length of the history does not exceed tau
    max_len = length(history.x) * 5;
    if size(history.u, 1) > max_len
        history.t = history.t(end-max_len+1:end);
        history.u = history.u(end-max_len+1:end, :, :);
    end
end
function plot_results(t, xspan, history)
    S = squeeze(history.u(:, :, 1));
    I = squeeze(history.u(:, :, 2));
    R = squeeze(history.u(:, :, 3));
    figure;
    subplot(1, 3, 1);
    surf(xspan, t, S);
    title('Susceptible (S)');
    xlabel('Distance (x)');
    ylabel('Time (t)');   
    shading interp;
    subplot(1, 3, 2);
    surf(xspan, t, I);
    title('Infectious (I)');
    xlabel('Distance (x)');
    ylabel('Time (t)'); 
    shading interp;
    subplot(1, 3, 3);
    surf(xspan, t, R);
    title('Recovered (R)');
    xlabel('Distance (x)');
    ylabel('Time (t)'); 
    shading interp;
end
function [c,f,s] = SIRPDE3(x, t, u, DuDx)
    global history;
    
    Lambda = 0.8; 
    gamma = 0.6; 
    mu = 0.3;
    alpha = 0.5;
    beta = 0.2;
    theta = 0.7;
    tau = 2.65;
    d1 = 0.1; d2 = 0.1; d3 = 0.1;
    S = u(1);
    I = u(2);
    R = u(3);
    u_tau = interpolate_history(history, x, t, tau);
    S_pre = u_tau(1);
    dSdt = Lambda - gamma * S_pre * I - mu * S;
    dIdt = theta * gamma * S_pre * I - mu * I - beta*I/(1+alpha*I);
    dRdt = (1 - theta) * gamma * S_pre * I - mu * R + beta*I/(1+alpha*I);
    c = [1; 1; 1];
    f = [d1; d2; d3] .* DuDx;
    s = [dSdt; dIdt; dRdt];
end

Perhaps you could give me some help and suggestions for changes? I would be grateful!

[Torsten]

The initial conditions in the calls to "pdepe" have to be set to the values of u at time t = (step-1)*tau. You always restart the solver at time t = 0.

And you don't need to concatenate the history data - you can simply replace them.

[鑫彤 纪]

Thank you very much for your reply!

After learning what the problem was, I modified it a few more times and found that in the end the modifications were not correct, and I couldn't even get it to run.

Can you help me to see how to modify this code? I would be grateful if you could help me solve this problem!

[Torsten]

Since you start with flat profiles for S, I and R and set the boundary conditions as dS/dx = 0, dI/dx = 0 and dR/dx = 0, the profiles will remain flat in the x-coordinate for all time t. Thus you don't need to solve a partial differential equation using "pdepe", but you only need MATLAB's delay solver "dde23"

( f = [d1; d2; d3] .* DuDx; can be neglected).

S0 = 1.25;
I0 = 0.85;
R0 = 0.74;
tau = 2.6;

xstart = 0;
xend = 10;
nx = 100;

tstart = 0;
tend = tau;
nt = 40;
ncycles = 265;

xspan = linspace(xstart,xend,nx);
tspan0 = linspace(tstart,tend,nt);
tspan = tspan0;

u_history = S0*ones(nt,nx);
SIRInitialConditions3 = @(x)[S0;I0;R0];

hold on
for i = 1:ncycles-1
  sol = pdepe(0, @(x, t, u, DuDx)SIRPDE3(x, t, u, DuDx, xspan, tau, tspan0, u_history, tspan(end)), SIRInitialConditions3, @SIRBoundaryConditions3, xspan, tspan);
  plot(tspan,sol(:,50,1),'r')
  plot(tspan,sol(:,50,2),'b')
  plot(tspan,sol(:,50,3),'m')
  tspan = tspan + tau;
  u_history = sol(:,:,1);
  SIRInitialConditions3 = @(x)[interp1(xspan,squeeze(sol(end,:,1)),x);interp1(xspan,squeeze(sol(end,:,2)),x);interp1(xspan,squeeze(sol(end,:,3)),x)];
end
hold off
grid on

function [c,f,s] = SIRPDE3(x, t, u, DuDx, xspan, tau, tspan0, u_history, tend)    
    Lambda = 0.8; 
    gamma = 0.6; 
    mu = 0.3;
    alpha = 0.5;
    beta = 0.2;
    theta = 0.7;
    d1 = 0.1; d2 = 0.1; d3 = 0.1;
    S = u(1);
    I = u(2);
    R = u(3);
    tref = min(max(0,t-tend+tau),tau);
    %[~,idx] = min((x-xspan).^2);
    %S_pre = interp1(tspan0,u_history(:,idx),tref);
    S_pre = interp2(xspan,tspan0,u_history,x,tref);
    dSdt = Lambda - gamma * S_pre * I - mu * S;
    dIdt = theta * gamma * S_pre * I - mu * I - beta*I/(1+alpha*I);
    dRdt = (1 - theta) * gamma * S_pre * I - mu * R + beta*I/(1+alpha*I);
    c = [1; 1; 1];
    f = [d1; d2; d3] .* DuDx;
    s = [dSdt; dIdt; dRdt];
end
function [pl, ql, pr, qr] = SIRBoundaryConditions3(xl, ul, xr, ur, t)
    pl = [0; 0; 0];
    ql = [1; 1; 1];
    pr = [0; 0; 0];
    qr = [1; 1; 1];
end

Edit:

I had to make a change to the original code because SIRPDE3 is not called at the spatial coordinates defined in "xspan", but at the corresponding midpoints.

[鑫彤 纪]

Thank you very much for your answer!

I understand what you are saying. But how should I write the code if I insist on wanting to see the result of a 3D image. My teacher wants to see a 3D result, but I've been running out of ideas and am very anxious.

I would appreciate if you could teach me how to change it!

[Torsten]

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.

[鑫彤 纪]

Ok, thank you very much! Your reply really helped me a lot!