Numerical Study of an SIR Epidemic Model with Vaccination: RK4 versus Adams–Bashforth–Moulton in Scilab

This dissertation studies a classical SIR (Susceptible–Infected–Recovered) epidemic model extended with vaccination from two angles: the mathematics of the model and the numerics of solving it.

On the analytical side, the model's equilibria, the vaccination-adjusted basic reproduction number and the local stability of the disease-free and endemic equilibria are derived symbolically in Maxima, the computer-algebra system taught in M106P. The Jacobian is computed and its eigenvalues are examined to obtain the threshold condition R_v < 1 and the critical vaccination rate.

On the numerical side, the same system is solved in Scilab with three schemes coded by hand — explicit Euler, the classical fourth-order Runge–Kutta (RK4) method and the fourth-order Adams–Bashforth–Moulton (ABM4) predictor–corrector — and compared against a high-accuracy reference solution computed with Scilab's built-in ode solver at tight tolerances. For each method we measure global error against step size, estimate the observed order of convergence, count function evaluations as a measure of cost and look at stability for large steps. Finally, we study how the epidemic peak (its height and timing) responds to the vaccination rate.

All parameter values are hypothetical and illustrative, chosen to be in a plausible range; the model is not fitted to real data. The dissertation is typeset in LaTeX and presented with Beamer (M404P), matching Bangalore University's Semester-4 Project Work.

Abstract

We analyse an SIR model with constant recruitment, natural death and vaccination of susceptibles. Using Maxima we derive the disease-free equilibrium (DFE), the endemic equilibrium and the reproduction number R_v = βμ / ((μ + ν)(γ + μ)), and we show via the Jacobian that the DFE is locally asymptotically stable when R_v < 1. We then compare three numerical schemes — Euler, RK4 and ABM4 — implemented in Scilab, measuring global error, observed order, computational cost and step-size stability against a reference solution. A parameter study shows how the vaccination rate ν shifts the height and timing of the infection peak.

Introduction

Compartmental models reduce an epidemic to a small system of nonlinear ordinary differential equations. They rarely have closed-form solutions, so numerical integration is essential, and the choice of method affects both accuracy and cost. The SIR model with vaccination is ideal for an M.Sc project: it is simple enough to analyse fully by hand and by computer algebra, yet rich enough to show real numerical behaviour such as error accumulation and step-size limits.

Literature and gap

Textbook treatments (Brauer and Castillo-Chavez; Murray) derive thresholds analytically but usually treat the numerical solution as a black box. Numerical-analysis texts (Burden and Faires; Atkinson) analyse RK and multistep methods on toy problems. This project connects the two: the epidemiological model is analysed rigorously and then used as the test problem for a careful, reproducible comparison of one-step and multistep methods.

Objectives in brief

Symbolic analysis in Maxima; implementation of three schemes in Scilab; convergence and cost study; sensitivity of the peak to ν; report in LaTeX and slides in Beamer.

Feasibility

Scilab, Maxima and a TeX distribution are free and already used in the M.Sc practicals. The computations run in seconds on a laptop. The work fits the 16-week, 8-hours-per-week self-study allotted to Project Work, and needs no data collection or ethics approval because all parameters are hypothetical.

Vaccination changes the dynamics of an infectious disease: it removes susceptibles directly and can push the reproduction number below one. For an SIR model with vaccination, the key mathematical questions are where the equilibria lie, under what condition the infection dies out and how strongly the vaccination rate changes the size and timing of an outbreak.

These questions are answered partly analytically and partly by numerical simulation, and the numerical answer is only as trustworthy as the method used. Explicit Euler is easy but inaccurate; RK4 is accurate but needs four function evaluations per step; ABM4 needs only two per step but is less stable and needs starting values. A student who reports "the peak occurs on day 42" without an error estimate has not really answered the question.

This project therefore (i) derives the equilibria, R_v and the stability conditions symbolically, and (ii) quantifies the accuracy, observed order, cost and stability of Euler, RK4 and ABM4 on this model against a high-accuracy reference, so that the conclusions about vaccination are supported by verified numerics.

In scope

• A deterministic, homogeneous-mixing SIRV-type model in normalised form (total population constant and scaled to 1).
• Symbolic analysis: equilibria, R_v, critical vaccination rate, Jacobian and local stability.
• Three hand-coded schemes (Euler, RK4, ABM4 with RK4 start-up) plus Scilab's ode as reference.
• Error-versus-step-size study, observed order, function-evaluation count, behaviour at large step sizes.
• Parameter sensitivity of the infection peak to the vaccination rate ν, using hypothetical parameters.

Out of scope

• Fitting the model to real surveillance data or making forecasts for any real disease.
• Stochastic, age-structured or spatial models.
• Global stability proofs by Lyapunov functions (mentioned as future scope).
• Adaptive step-size control beyond what the reference solver provides.

Research design: analytical derivation followed by controlled numerical experiments (a computational study).

Model (normalised, N = 1)

dS/dt = μ − βSI − (μ + ν)S dI/dt = βSI − (γ + μ)I dR/dt = γI + νS − μR

β: transmission rate, γ: recovery rate, μ: birth = death rate, ν: vaccination rate of susceptibles. Because S + I + R = 1 is invariant, R can be recovered from S and I, and stability analysis is done on the (S, I) subsystem.

Analytical results to derive

• DFE: E0 = (μ/(μ+ν), 0, ν/(μ+ν)).
• Reproduction number: R_v = βμ / ((μ+ν)(γ+μ)); without vaccination R0 = β/(γ+μ).
• Critical vaccination rate: R_v < 1 ⇔ ν > μ(R0 − 1).
• Endemic equilibrium (exists when R_v > 1): S* = (γ+μ)/β, I* = (μ+ν)(R_v − 1)/β.

Illustrative (hypothetical) parameters

β = 0.5 per day, γ = 0.1 per day, μ = 1/(60 × 365) per day, ν ∈ {0, 0.0005, 0.001, 0.002, 0.005}, initial state S = 0.999, I = 0.001, R = 0, horizon 0–365 days. With these values R0 ≈ 5; they are chosen only to produce a visible outbreak and are not estimates for any real disease.

Numerical experiments

1. Reference solution: Scilab ode with rtol = atol = 1e−12 on a fine output grid.
2. For h ∈ {2, 1, 0.5, 0.25, 0.125, 0.0625} days, run Euler, RK4 and ABM4; compute E(h) = max over grid of ‖y_h − y_ref‖∞.
3. Observed order p ≈ log2(E(h)/E(h/2)); theory predicts p ≈ 1 (Euler) and p ≈ 4 (RK4, ABM4).
4. Cost: total right-hand-side evaluations (Euler 1, RK4 4, ABM4-PECE 2 per step) vs error, plotted on log–log axes.
5. Stability: increase h until each method's solution becomes negative or oscillates; relate to the linearised eigenvalues.
6. Peak study: for each ν, record I_max and t_peak with RK4 at a verified step size.

Timeline (16 weeks, individual)