Numerically solving a system of partial integro-differential equations in Matlab

Question

Given the following system of partial integro-differential equations -

$\frac{dS(t)}{dt}=\Lambda-\mu S(t)-\beta S(t)F(t),\\ \frac{\partial I(t,\omega)}{\partial t}+\frac{\partial I(t,\omega)}{\partial \omega}=-\delta I(t,\omega),\\ \frac{dF(t)}{dt}=-\gamma F(t)+(N-F(t))\int^{\infty}_{0}\kappa (\omega)I(t,\omega)d\omega$

with the boundary condition, $I(t,0)=\beta S(t)F(t)$ and

initial conditions $S(0)=S_{0},I(0,\omega)=0,F(0)=F_{0}$ where all parameters $\Lambda,\mu,\beta,\delta,\gamma,N,\kappa(\omega),S_{0},F_{0}$ are positive and $\kappa(\omega)$ can be a simple fourth order polynomial.

How does one numerically compute a solution for this on Matlab?

Answer 1

Well, one might try following scheme. Let $$ I^{n}_j = I(t_n, \omega_j), \quad S^n = S(t_n), \quad F^n = F(t_n), \quad t_n = n \tau, \quad \omega_j = j h $$ Note that $I(t,\omega) = 0$ when $\omega > t$. So let's limit the computational domain to $0 \leq t \leq T, 0 \leq \omega \leq L \geq T$. Then $$\begin{aligned} \frac{S^{n+1} - S^n}{\tau} &= \Lambda - \mu S^{n+1} - \beta S^{n+1} F^{n+1}\\ \frac{F^{n+1} - F^n}{\tau} &= -\gamma F^{n+1} + (N - F^{n+1}) h\sum_{j=0}^{M-1} \frac{\kappa(\omega_j) I^n_j + \kappa(\omega_{j+1}) I^n_{j+1}}{2}\\ \frac{I^{n+1}_j - I^n_j}{\tau} + \frac{I^{n}_j - I^n_{j-1}}{h} &= -\delta I^{n+1}_j, \qquad j = 1, 2, \dots, M\\ I^{n+1}_0 &= \beta S^{n+1} F^{n+1}\\ I^0_j &= 0, \quad j = 0,1, \dots, M\\ S^0 &= S_0, \quad F^0 = F_0 \end{aligned} $$ The third equation is an upwind approximation for advection equation. The scheme is stable if $\tau \leq h$. For integral term the trapeziodal rule is used. Right-hand terms of the ordinary DE are approximated implcitly, to resolve possible instabilities if system is stiff.

The algorithm:

• Assume that every value is already known at $t_n$. Initial data is given for $n = 0$.
• First, compute $$ J^n = h\sum_{j=0}^{M-1} \frac{\kappa(\omega_j)I^n_j + \kappa(\omega_{j+1})I^n_{j+1}}{2} $$
• Compute $$ F^{n+1} = \frac{F^n + \tau N J^n}{1 + \tau \gamma + \tau J^n} $$
• Next $$ S^{n+1} = \frac{S^n + \tau \Lambda}{1 + \tau \mu + \tau \beta F^{n+1}} $$
• Now $$ I^{n+1}_0 = \beta S^{n+1} F^{n+1} $$
• And finally, $$ I^{n+1}_j = \frac{(1 - \sigma)I^n_j + \sigma I^n_{j-1}}{1 + \tau \delta}, \quad j = 1, 2, \dots, M, \quad \sigma = \frac{\tau}{h}. $$ Since your advection part has constant speed, you can greatly improve its performance by taking $\sigma = 1$, i.e. $\tau = h$.

Here's Matlab code which implements algorithm above

T = 10;
L = 15;
M = 200;

h = L / M;
sigma = 1;
tau = sigma * h;

Lambda = 1;
mu = 1;
beta = 1;
gamma = 1;
delta = .1;
N = 1;
S0 = 1;
F0 = 1;
kappa = @(w) w.^3 + w;

I = zeros(M + 1, 1);
t = 0;

trapez = ones(M + 1, 1) * h;
trapez(1) = h / 2;
trapez(end) = h / 2;

omega = (0:M) * h;
omega = omega(:);

F = F0;
S = S0;
step = 0;

while t < T
    step = step + 1;
    J = dot(kappa(omega) .* I, trapez);
    F = (F + tau * N * J) / (1 + tau * gamma + tau * J);
    S = (S + tau * Lambda) / (1 + tau * mu + tau * beta * F);
    Inext = I;
    Inext(1) = beta * S * F;
    Inext(2:end) = ((1 - sigma) * I(2:end) + sigma * I(1:end-1)) ...
        / (1 + tau * delta);
    I = Inext;
    if (mod(step, 10) == 0)
        plot(omega, I, '.-');
        axis([0 L 0 1]);
        pause(.1);
    end
    t = t + tau;
end