# Numerically solving a system of partial integro-differential equations in Matlab [closed]

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?

## closed as too broad by user147263, Did, Ivo Terek, user91500, Claude LeiboviciJul 19 '15 at 4:20

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer. Avoid asking multiple distinct questions at once. See the How to Ask page for help clarifying this question. If this question can be reworded to fit the rules in the help center, please edit the question.

• Are you aware of the PDE Toolbox for Matlab? In any case this Question might be better asked at scicomp.se, with a link back to your earlier post here. – hardmath Jul 29 '15 at 3:54

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

• Thank you. Am I correct to assume that "L" is the upper limit of integration (instead of ∞)? – squeak Jul 20 '15 at 6:19
• @squeak Yes, added that to the post. We can take $L$ to anything greater than $T$. In case when $\sigma \neq 1$, the solution of the advection part would have numerical diffusion, so numerical soulution violates $I(t,\omega) = 0 \text{ when } \omega > t$ sligtly. So I took $L$ slightly greater than $T$. – uranix Jul 20 '15 at 13:10