🦠 Influenza SVEIR PDE Model

Abstract

This report addresses the numerical solution of a one-dimensional partial differential equation (PDE) model describing the spatial and temporal dynamics of an influenza epidemic with vaccination. The model follows the SVEIR compartmental structure, representing the populations of Susceptible, Vaccinated, Exposed, Infected, and Recovered individuals as functions of both space and time. Originally solved using the Method of Lines (MOL), we extend the analysis by applying a range of other methods, including several finite difference schemes, a Strang splitting method, finite volume method, spectral method and a deep learning approach. The goal is to compare the accuracy and efficiency of each method in modeling influenza spread with vaccination and diffusion.

I. Introduction

Influenza transmission is governed by multiple biological processes including infection, recovery, immunity loss, and vaccination. While classical compartmental models typically address only the temporal dynamics of disease spread, spatial heterogeneity—such as population movement and localized outbreaks—can significantly influence epidemic behavior. To capture these effects, this work models influenza using a spatiotemporal SVEIR framework, represented by a system of coupled nonlinear partial differential equations (PDEs). The model tracks the densities of Susceptible (S), Vaccinated (V), Exposed (E), Infected (I), and Recovered (R) individuals as functions of both time and space.

Originally solved using the Method of Lines (MOL), which converts the PDE system into ODEs by discretizing space, the model provides insight into how vaccination and diffusion influence epidemic spread. In this project, we extend the analysis by comparing six approaches: the Explicit Finite Difference Method, the Crank–Nicolson Method, the Strang Splitting Method, Finite Volume Method, Spectral Method and a Deep Learning model trained to approximate PDE solutions.

II. Literature Review

Recent studies have advanced epidemic modeling by incorporating spatial and temporal effects. Chen and Xu (2019) proposed a two-strain influenza model with vaccines, including spatial diffusion and time delays to reflect incubation and movement [1]. They analyzed equilibrium stability using Lyapunov methods and confirmed results numerically. Jang et al. (2020) addressed optimal vaccine allocation in a spatial SIR model, using constrained optimization to minimize infections and costs, highlighting the impact of spatial heterogeneity [2]. Building on these works, we explore numerical and data-driven methods for solving a spatial SVEIR influenza model with vaccination and diffusion.

III. Model Formulation

This study models the spread of influenza using a one-dimensional reaction–diffusion system based on the SVEIR framework. The model captures both temporal and spatial dynamics of disease transmission and vaccination. Five state variables are tracked as functions of space \( x \in [-3, 3] \) and time \( t \geq 0 \):

• S(x, t): Susceptible population
• V(x, t): Vaccinated population
• E(x, t): Exposed (latent) population
• I(x, t): Infected population
• R(x, t): Recovered population

PDE System

The system is governed by nonlinear partial differential equations representing infection, recovery, immunity loss, and diffusion:

$$ \frac{\partial S}{\partial t} = -\beta\beta_E E S - \beta\beta_I I S + \alpha I S - \phi S - rS + \delta R + \theta V + r + d_1 \frac{\partial^2 S}{\partial x^2} $$

$$ \frac{\partial V}{\partial t} = -\beta\beta_E\beta_V E V - \beta\beta_I\beta_V I V + \alpha I V - rV - \theta V + \phi S + d_2 \frac{\partial^2 V}{\partial x^2} $$

$$ \frac{\partial E}{\partial t} = \beta\beta_E E S + \beta\beta_I I S + \beta\beta_E\beta_V E V + \beta\beta_I\beta_V I V + \alpha I E - (r + \kappa + \sigma) E + d_3 \frac{\partial^2 E}{\partial x^2} $$

$$ \frac{\partial I}{\partial t} = \sigma E - (r + \alpha + \gamma) I + \alpha I^2 + d_4 \frac{\partial^2 I}{\partial x^2} $$

$$ \frac{\partial R}{\partial t} = \kappa E + \gamma I - rR - \delta R + \alpha I R + d_5 \frac{\partial^2 R}{\partial x^2} $$

Boundary and Initial Conditions

Neumann (no-flux) boundary conditions are applied at \( x = \pm 3 \) to simulate a closed domain:

$$ \left. \frac{\partial U}{\partial x} \right|_{x = -3}^{x = 3} = 0, \quad \text{for } U = S, V, E, I, R $$

Initial conditions are defined as:

$$ \begin{align*} S(x, 0) &= 0.86 \exp\left(-\left(\frac{x}{1.4}\right)^2\right), \\ V(x, 0) &= 0.10 \exp\left(-\left(\frac{x}{1.4}\right)^2\right), \\ E(x, 0) &= 0, \quad I(x, 0) = 0.04 \exp(-x^2), \quad R(x, 0) = 0 \end{align*} $$

Parameter Values

IV. Numerical Solutions

This section presents the various numerical techniques applied to solve the spatiotemporal SVEIR model for influenza transmission.

Crank–Nicolson Numerical Method

We used a semi-implicit Crank–Nicolson method to solve the SVEIR PDE model. Diffusion terms are treated implicitly, while nonlinear reaction terms are computed explicitly. This requires solving a tridiagonal system for each compartment at every time step.

Spatial and Temporal Grids

The spatial domain is defined as \( x \in [x_l, x_u] \) and divided into \( n_x \) uniform grid points:

$$ \Delta x = \frac{x_u - x_l}{n_x - 1}, \quad \Delta x^2 = (\Delta x)^2 $$

Initial Conditions

$$ \begin{aligned} S(x, 0) &= 0.86 \exp\left(-\left(\frac{x}{1.4}\right)^2\right), \\ V(x, 0) &= 0.10 \exp\left(-\left(\frac{x}{1.4}\right)^2\right), \\ E(x, 0) &= 0, \quad I(x, 0) = 0.04 \exp(-x^2), \quad R(x, 0) = 0 \end{aligned} $$

Crank–Nicolson Matrix Construction

For each compartment, define a diffusion coefficient \( d_i \) (e.g., \( d_1 = 0.05 \) for \( S \)). Compute:

$$ \lambda = \frac{d_i \, \Delta t}{2 \Delta x^2} $$

Then construct the tridiagonal matrix \( A \) as:

$$ A = \text{diag}(1 + 2\lambda) + \text{diag}(-\lambda, 1) + \text{diag}(-\lambda, -1) $$

Time-Stepping Scheme

For each time step \( k = 1 \) to \( N_t - 1 \), extract the current state vectors \( S^k, V^k, E^k, I^k, R^k \), and compute the nonlinear terms (reaction terms):

$$ \begin{aligned} F_S &= -\beta \beta_E E^k S^k - \beta \beta_I I^k S^k + \alpha I^k S^k - \phi S^k - r S^k + \delta R^k + \theta V^k \\ F_V &= -\beta \beta_E \beta_V E^k V^k - \beta \beta_I \beta_V I^k V^k + \alpha I^k V^k - r V^k - \theta V^k + \phi S^k \\ F_E &= \beta \beta_E E^k S^k + \beta \beta_I I^k S^k + \beta \beta_E \beta_V E^k V^k + \beta \beta_I \beta_V I^k V^k + \alpha E^k I^k - (r + \kappa + \sigma) E^k \\ F_I &= \sigma E^k - (r + \alpha + \gamma) I^k + \alpha (I^k)^2 \\ F_R &= \kappa E^k + \gamma I^k - r R^k - \delta R^k + \alpha I^k R^k \end{aligned} $$

Update Step

Update each compartment using the Crank–Nicolson scheme:

$$ U^{k+1} = A^{-1} \left(U^k + \Delta t \cdot F_U\right), \quad \text{for } U \in \{S, V, E, I, R\} $$

Remarks

This method is efficient due to the tridiagonal nature of matrix \( A \), allowing fast solutions. Explicit handling of the nonlinear terms avoids solving a nonlinear system, but may slightly reduce accuracy. Overall, Crank–Nicolson offers a good balance between stability and computational cost for reaction–diffusion models.

Explicit Finite Difference Method

To solve the spatiotemporal SVEIR model numerically, we first implemented the Explicit Forward Euler Method, a basic time-stepping scheme. This method discretizes both time and space and advances the solution iteratively. It is easy to implement but conditionally stable, requiring careful choice of time step size.

Space–Time Grid

The spatial domain is defined as \( x \in [x_l, x_u] \) and divided into \( n_x \) uniform grid points:

$$ x_i = x_l + (i-1)\Delta x, \quad \Delta x = \frac{x_u - x_l}{n_x - 1}, \quad i = 1, 2, \ldots, n_x $$

Time is discretized using a uniform time step \( \Delta t \), up to a final time \( t_f \):

$$ t_n = t_0 + n\Delta t, \quad n = 0, 1, 2, \ldots, N_t $$

Initial Conditions

$$ \begin{aligned} S(x, 0) &= 0.86 \exp\left(-\left(\frac{x}{1.4}\right)^2\right), \\ V(x, 0) &= 0.10 \exp\left(-\left(\frac{x}{1.4}\right)^2\right), \\ E(x, 0) &= 0, \quad I(x, 0) = 0.04 \exp(-x^2), \quad R(x, 0) = 0 \end{aligned} $$

Laplacian Approximation

The second derivative in space is approximated using central differences:

$$ \frac{\partial^2 u}{\partial x^2} \approx \frac{u_{i+1} - 2u_i + u_{i-1}}{(\Delta x)^2} $$

This approximation is applied to each compartment to capture diffusion effects.

Explicit Update Equations

At each time step \( t_n \), we compute the solution at \( t_{n+1} \) using the Forward Euler update:

$$ \begin{aligned} S_i^{n+1} &= S_i^n + \Delta t \left( -\beta\beta_E E_i^n S_i^n - \beta\beta_I I_i^n S_i^n + \alpha I_i^n S_i^n - \phi S_i^n - r S_i^n + \delta R_i^n + \theta V_i^n + r + d_1 \nabla^2 S_i^n \right) \\ V_i^{n+1} &= V_i^n + \Delta t \left( -\beta\beta_E\beta_V E_i^n V_i^n - \beta\beta_I\beta_V I_i^n V_i^n + \alpha I_i^n V_i^n - r V_i^n - \theta V_i^n + \phi S_i^n + d_2 \nabla^2 V_i^n \right) \\ E_i^{n+1} &= E_i^n + \Delta t \left( \beta\beta_E E_i^n S_i^n + \beta\beta_I I_i^n S_i^n + \beta\beta_E\beta_V E_i^n V_i^n + \beta\beta_I\beta_V I_i^n V_i^n + \alpha I_i^n E_i^n - (r + \kappa + \sigma) E_i^n + d_3 \nabla^2 E_i^n \right) \\ I_i^{n+1} &= I_i^n + \Delta t \left( \sigma E_i^n - (r + \alpha + \gamma) I_i^n + \alpha (I_i^n)^2 + d_4 \nabla^2 I_i^n \right) \\ R_i^{n+1} &= R_i^n + \Delta t \left( \kappa E_i^n + \gamma I_i^n - r R_i^n - \delta R_i^n + \alpha I_i^n R_i^n + d_5 \nabla^2 R_i^n \right) \end{aligned} $$

Here, \( \nabla^2 U_i^n \) represents the discrete Laplacian for variable \( U \) at grid point \( i \) and time \( n \).

Stability Considerations

This method is conditionally stable. The stability is especially sensitive to the diffusion terms, requiring a sufficiently small \( \Delta t \) based on the spatial resolution \( \Delta x \) and diffusivity values \( d_j \). While computationally efficient for small time intervals and moderate grid sizes, it can become unstable if the time step is too large.

Strang Splitting Method

This method decouples the reaction and diffusion components of the system. This operator splitting technique enables each physical process—reaction and diffusion—to be handled using the most suitable numerical method, enhancing both stability and efficiency.

Governing Equations

Each equation includes nonlinear reaction terms and spatial diffusion:

$$ \frac{\partial u_i}{\partial t} = \text{Reaction}_i(u_1, u_2, u_3, u_4, u_5) + d_i \frac{\partial^2 u_i}{\partial x^2}, \quad i = 1, 2, \ldots, 5 $$

where \( u_i \in \{S, V, E, I, R\} \), and \( d_i \) is the diffusion coefficient for each compartment.

Spatial and Temporal Discretization

The spatial domain is defined as \( x \in [x_l, x_u] \) and discretized into \( n_x \) grid points:

$$ x_i = x_l + (i - 1)\Delta x, \quad \Delta x = \frac{x_u - x_l}{n_x - 1}, \quad i = 1, 2, \ldots, n_x $$

Time is discretized similarly, using a uniform time step \( \Delta t \), up to a final time \( t_f \):

$$ t_n = t_0 + n \Delta t, \quad n = 0, 1, 2, \ldots, N_t $$

Initial Conditions

The variable initializations are identical to those used in the Crank–Nicolson method, following the same setup as the original model:

$$ \begin{aligned} S(x, 0) &= 0.86 \exp\left(-\left(\frac{x}{1.4}\right)^2\right), \\ V(x, 0) &= 0.10 \exp\left(-\left(\frac{x}{1.4}\right)^2\right), \\ E(x, 0) &= 0, \quad I(x, 0) = 0.04 \exp(-x^2), \quad R(x, 0) = 0 \end{aligned} $$

Splitting Strategy

Each full time step \( \Delta t \) is advanced using the following symmetric sequence:

Step 1: Half-Step Reaction

The reaction terms are integrated using a classical fourth-order Runge–Kutta method over a half-step \( \Delta t / 2 \).

Step 2: Full-Step Diffusion

The diffusion terms are updated using explicit forward Euler with substeps. The second spatial derivative is approximated using central differences:

$$ \frac{\partial^2 u}{\partial x^2} \approx \frac{u_{i+1} - 2u_i + u_{i-1}}{(\Delta x)^2} $$

The diffusion is computed with 20 substeps of size:

$$ \Delta t_{\text{sub}} = \frac{\Delta t}{20} $$

Step 3: Half-Step Reaction

A second Runge–Kutta integration is performed over the remaining half-step to complete the time advancement.

Remarks

This method allows flexibility by using different solvers suited to each component (reaction vs. diffusion). It is well-suited for nonlinear, spatially extended epidemic models, as it balances accuracy, speed, and implementation simplicity.

Finite Volume Method (FVM) with RK4

1. Overview

We solve the SVEIR reaction–diffusion PDE model for influenza spread using the Finite Volume Method (FVM) for spatial discretization. This method conserves flux across control volumes and handles diffusion and reaction terms efficiently. Time integration is performed using the classical fourth-order Runge–Kutta method (RK4).

2. SVEIR PDE Model

The equations for the five compartments are:

\[ \begin{aligned} \frac{\partial S}{\partial t} &= -\beta\beta_e E S - \beta\beta_i I S + \alpha I S - \phi S - rS + \delta R + \theta V + r + d_1 \nabla^2 S \\ \frac{\partial V}{\partial t} &= -\beta\beta_e\beta_v E V - \beta\beta_i\beta_v I V + \alpha I V - rV - \theta V + \phi S + d_2 \nabla^2 V \\ \frac{\partial E}{\partial t} &= \beta\beta_e E S + \beta\beta_i I S + \beta\beta_e\beta_v E V + \beta\beta_i\beta_v I V + \alpha I E - (r + \kappa + \sigma) E + d_3 \nabla^2 E \\ \frac{\partial I}{\partial t} &= \sigma E - (r + \alpha + \gamma) I + \alpha I^2 + d_4 \nabla^2 I \\ \frac{\partial R}{\partial t} &= \kappa E + \gamma I - rR - \delta R + \alpha I R + d_5 \nabla^2 R \end{aligned} \]

Here, \( \nabla^2 u \) represents the Laplacian term, approximated using the FVM.

3. Space–Time Grid

The spatial domain \( x \in [x_l, x_u] \) is divided into \( n_x \) uniform cells:

\[ x_i = x_l + (i-1)\Delta x, \quad \Delta x = \frac{x_u - x_l}{n_x - 1} \]

Time is discretized as:

\[ t_n = t_0 + n\Delta t, \quad n = 0,1,\dots,N_t \]

4. Initial Conditions

\[ \begin{aligned} S(x,0) &= 0.86 \exp\left(-\left(\frac{x}{1.4}\right)^2\right), \quad V(x,0) = 0.10 \exp\left(-\left(\frac{x}{1.4}\right)^2\right), \\ I(x,0) &= 0.04 \exp(-x^2), \quad E(x,0) = 0, \quad R(x,0) = 0 \end{aligned} \]

5. Laplacian Approximation (FVM)

The second spatial derivative is approximated as:

\[ \nabla^2 u_i \approx \frac{u_{i+1} - 2u_i + u_{i-1}}{(\Delta x)^2} \]

Boundary conditions (Neumann, zero-flux):

\[ \nabla^2 u_1 \approx \frac{u_2 - u_1}{(\Delta x)^2}, \quad \nabla^2 u_n \approx \frac{u_{n-1} - u_n}{(\Delta x)^2} \]

6. Time Integration (RK4)

We use RK4 to integrate the system over time:

\[ y^{n+1} = y^n + \frac{\Delta t}{6}(k_1 + 2k_2 + 2k_3 + k_4) \]

Each stage \( k \) is computed using the discretized RHS of the PDE, including diffusion (FVM) and reaction terms.

7. Update Equations

At each grid point \( i \) and time step \( n \):

\[ \begin{aligned} S_i^{n+1} &= S_i^n + \Delta t \cdot f_S(S^n, V^n, E^n, I^n, R^n) \\ V_i^{n+1} &= V_i^n + \Delta t \cdot f_V(S^n, V^n, E^n, I^n, R^n) \\ E_i^{n+1} &= E_i^n + \Delta t \cdot f_E(S^n, V^n, E^n, I^n, R^n) \\ I_i^{n+1} &= I_i^n + \Delta t \cdot f_I(S^n, V^n, E^n, I^n, R^n) \\ R_i^{n+1} &= R_i^n + \Delta t \cdot f_R(S^n, V^n, E^n, I^n, R^n) \end{aligned} \]

Where each \( f \) function includes both reaction and diffusion terms.

8. Visualization

We visualize the results:

• As functions of space at times \( t = 0, 6, 12, ..., 60 \)
• As functions of time at the spatial midpoint \( x = 0 \)

9. Conclusion

The Finite Volume Method offers a conservative and accurate framework for simulating disease spread with spatial heterogeneity. Combined with the RK4 integrator, it yields reliable solutions to the nonlinear spatiotemporal SVEIR model.

Spectral Method (DCT-Based) with lsoda

1. Overview

We present a numerical solution to the SVEIR reaction–diffusion model for influenza using a spectral method. The model consists of five coupled partial differential equations (PDEs) for the compartments: Susceptible (S), Vaccinated (V), Exposed (E), Infected (I), and Recovered (R). We use the Discrete Cosine Transform (DCT) to approximate spatial derivatives, while time integration is performed using the adaptive lsoda solver.

2. SVEIR Spectral PDE Model

The governing PDEs with spectral Laplacians are:

\[ \begin{aligned} \frac{\partial S}{\partial t} &= -\beta\beta_e E S - \beta\beta_i I S + \alpha I S - \phi S - rS + \delta R + \theta V + r + d_1 L(S) \\ \frac{\partial V}{\partial t} &= -\beta\beta_e\beta_v E V - \beta\beta_i\beta_v I V + \alpha I V - rV - \theta V + \phi S + d_2 L(V) \\ \frac{\partial E}{\partial t} &= \beta\beta_e E S + \beta\beta_i I S + \beta\beta_e\beta_v E V + \beta\beta_i\beta_v I V + \alpha I E - (r + \kappa + \sigma) E + d_3 L(E) \\ \frac{\partial I}{\partial t} &= \sigma E - (r + \alpha + \gamma) I + \alpha I^2 + d_4 L(I) \\ \frac{\partial R}{\partial t} &= \kappa E + \gamma I - rR - \delta R + \alpha I R + d_5 L(R) \end{aligned} \]

Here, \( L(u) \) denotes the Laplacian computed via the spectral approximation:

\[ L(u) \approx \text{IDCT}\left(-\left(\frac{k\pi}{2L}\right)^2 \cdot \text{DCT}(u)\right) \]

where \( k = 0, 1, \dots, N-1 \) are the spectral modes, and \( L = \frac{x_u - x_l}{2} \) is half the spatial domain width.

3. Space–Time Grid

The spatial domain \( x \in [x_l, x_u] \) is uniformly divided into \( n_x \) grid points:

\[ x_i = x_l + (i - 1)\Delta x, \quad \Delta x = \frac{x_u - x_l}{n_x - 1}, \quad i = 1, \dots, n_x \]

Time domain: \( t_n = t_0 + n\Delta t \), solved adaptively at selected output points using lsoda.

4. Initial Conditions

\[ \begin{aligned} S(x,0) &= 0.86 \exp\left(-\left(\frac{x}{1.4}\right)^2\right), \quad V(x,0) = 0.10 \exp\left(-\left(\frac{x}{1.4}\right)^2\right), \\ I(x,0) &= 0.04 \exp(-x^2), \quad E(x,0) = 0, \quad R(x,0) = 0 \end{aligned} \]

5. Spectral Approximation Steps

1. Apply the DCT: \( \hat{u}_k = \text{DCT}(u) \)
2. Scale in frequency domain: \( \hat{u}_{xx,k} = -\left(\frac{k\pi}{2L}\right)^2 \hat{u}_k \)
3. Apply the inverse DCT: \( u_{xx} = \text{IDCT}(\hat{u}_{xx,k}) \)

This yields high-accuracy approximations under Neumann (zero-flux) boundary conditions.

6. Time Integration

After transforming the PDE system into a system of ODEs via the DCT-based Laplacian, we solve it using lsoda from the deSolve R package. This solver automatically switches between stiff and non-stiff methods, offering both stability and speed.

7. Output

Solutions \( S(x,t), V(x,t), E(x,t), I(x,t), R(x,t) \) are reconstructed for:

• Fixed time snapshots \( t = 0, 6, \dots, 60 \)
• Fixed spatial point (e.g., \( x = 0 \)) across time

8. Conclusion

The spectral method delivers highly accurate spatial approximations for smooth problems. When combined with robust adaptive time-stepping (like lsoda), it efficiently simulates the dynamics of influenza spread in a spatially structured population.

V. Deep Learning Solution

A. Introduction

This part of the project uses Deep Learning (DL) techniques, specifically Physics-Informed Neural Networks (PINNs), to solve a system of partial differential equations (PDEs) modeling the spread of influenza with vaccination and spatial diffusion. Unlike traditional numerical methods (such as Finite Difference or Finite Element), DL-based methods do not require mesh generation or discretization and can generalize across space–time domains.

B. What Deep Learning Is Solving

The neural network learns to approximate the five functions:

$$ S(x, t), \quad V(x, t), \quad E(x, t), \quad I(x, t), \quad R(x, t) $$

which represent the densities of susceptible, vaccinated, exposed, infected, and recovered individuals, respectively. The network is trained to minimize the mismatch in PDE residuals and satisfy initial and boundary conditions.

C. Data and Input Representation

• Input space: a 2D grid of points \((x, t)\) generated using PyTorch functions like torch.linspace and meshgrid.
• Each network input is a 2D vector \([x, t]\).
• Only 20 labeled points from the true numerical solution are used to guide training.

D. The Neural Network (PINN)

• Inputs: \( x \) and \( t \)
• Hidden layers: Fully connected (dense) layers with activation functions like Tanh, Sigmoid, or ReLU
• Output: 5 values — \( [S, V, E, I, R] \)

The network is trained to satisfy the following conditions:

1. The PDE system (i.e., the equations governing disease spread)
2. Initial conditions (population values at \( t = 0 \))
3. Boundary conditions (Neumann/no-flux: \( \frac{\partial U}{\partial x} = 0 \))
4. Supervision on 20 data points to improve convergence

E. Loss Function

The total loss is composed of multiple mean squared error (MSE) terms:

$$ \text{Total Loss} = \text{MSE}_{\text{PDE}} + \text{MSE}_{\text{IC}} + \text{MSE}_{\text{BC}} + \text{MSE}_{\text{Data}} $$

• \( \text{MSE}_{\text{PDE}} \): Residuals of the PDEs (physics loss)
• \( \text{MSE}_{\text{IC}} \): Deviation from initial conditions
• \( \text{MSE}_{\text{BC}} \): Enforcing zero-flux boundary condition
• \( \text{MSE}_{\text{Data}} \): Error on known true points

F. Optimization Strategy

Training is performed in two stages:

1. Stage 1: Adam optimizer (adaptive, fast convergence)
2. Stage 2: LBFGS optimizer (second-order, fine-tunes weights for better precision)

G. Challenges and Notes

While PINNs offer flexibility, they also face limitations. Some compartments—particularly \( I(x, t) \) and \( E(x, t) \)—may exhibit inaccuracies or oscillations. This can be caused by:

• Imbalanced weights among loss terms
• Insufficient model complexity (too few neurons/layers)
• Inaccurate gradients or poor initialization

Adjusting weight terms and increasing the number of training data points can improve performance. Regularization and adaptive learning rates may also help stabilize learning.

VI. Results

This section presents and compares the spatial and temporal behavior of the S, V, E, I and R populations under different numerical methods. The reference solution obtained using the Method of Lines (MOL) serves as a benchmark for comparison against the Crank–Nicolson, Explicit FD, FVM, Spectral, Strang Splitting, and Deep Learning approaches.

A. Numerical Results

1. Reference Solution: Method of Lines

The Method of Lines discretizes space and solves the resulting system of ODEs using standard time integration techniques. The following plot illustrates the spatial distribution of the Susceptible and Recovered populations at regular time steps.

Each compartment displays smooth, symmetric wave-like patterns, consistent with diffusion and biological dynamics.

We also track each compartment's value over time at the spatial center (x = 0), revealing typical epidemic patterns and delays.

2. Crank–Nicolson Method

3. Explicit FD Method

4. Strang Splitting Method

5. Finite Volume Method

6. Spectral Method

B. Deep Learning Results