# Readi2solve: An Online Python Toolbox for Solving Reaction-Diffusion problems

Lionel Roques

<figure><img src="image_1.png" /></figure>

## To cite this version:

Lionel Roques. Readi2solve: An Online Python Toolbox for Solving Reaction-Diffusion problems. 2024. hal-04728969

HAL Id: hal-04728969

https://hal.inrae.fr/hal-04728969v1

Preprint submitted on 9 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Readi2solve: An Online Python Toolbox for Solving Reaction-Diffusion problems

Lionel Roques

INRAE, BioSP, 84914, Avignon, France

October 9, 2024

## Abstract

"Readi2solve" is an online Python toolbox designed to solve various partial differential equations (PDEs) such as reaction-diffusion systems with nonlocal terms and principal eigenvalue problems. Hosted on Google Colab, the toolbox allows users to run code interactively without local installations. This document introduces the toolbox, presents its features, and provides examples of its use.

## Disclaimer

The purpose of this toolbox is to solve simple models and assist mathematicians in exploring theoretical problems. It helps in formulating conjectures for analytical testing and in generating figures or animations to illustrate theoretical results. It is not intended for simulations on realistic domains with complex geometry. For such tasks, please refer to the MSE project (http://mse.biosp.org).

## 1 Introduction

Partial differential equations (PDEs) are fundamental in the modeling of many physical, biological, and chemical processes. The "Readi2solve" toolbox provides a Python-based framework for solving several types of PDEs, including reaction-diffusion equations with nonlocal terms and principal eigenvalue problems. It also offers the feature of computing spreading speeds in reaction-diffusion equations, based on the Freidlin-Gärtner formula. The toolbox, available at http://readi2solve.biosp.org is designed to run in Google Colab, making it accessible to anyone with a Google account. The results can be visualized using interactive plots and animations.

This preprint describes the toolbox in detail and illustrates its usage with several examples.

## 2 Reaction-Diffusion Equations

The reaction-diffusion equations are solved using a method that combines spatial discretization with a time-integration technique. The spatial domain is discretized into a uniform grid using a finite difference scheme. The time evolution of the system is
computed using the 'odeint' solver from the SciPy library, which automatically chooses between two different methods depending on the stiffness of the system: the Adams-Bashforth-Moulton (ABM) method for non-stiff problems and the Backward Differentiation Formula (BDF) method for stiff problems.

## 2.1 Scalar equations in 1D

The toolbox solves scalar equation of the form

$$ \partial_t u = \partial_{xx}[D(t,x)u] + F \left( t, x, u, \int_{-L}^{L} u \right), \quad t > 0, x \in (-L, L), $$

with Dirichlet boundary condition $u(t, \pm L) = 0$ or Neumann condition $\partial_x(D u)(t, \pm L) = 0$ and initial condition $u_0(x)$.

### Example of use: Replicator-Mutator Equation

$$ \partial_t u = D \partial_{xx} u + u(1-u), \quad t > 0, x \in (-L, L), $$

with Dirichlet boundary conditions $u(t, \pm L) = 0$.

The following code demonstrates how to solve this problem (note that $\rho$ corresponds to the integral of $u$):

```matlab
# Example of functions
t, x, u, rho= symbols('t x u rho')
fun_u0 = exp(-(x+1)**2) * u
fun_F = (5 - x)**2 * u - u * rho
fun_D = 0.1

times = np.linspace(0, 50, 100) # times for animation

# Obtain numerical solution
numericalSolution, xx, _ = sol RD_1D(times, fun_D, fun_F, fun_u0, "Neumann")

# Plot the solution
animplot PDE1D(numericalSolution, xx, times)
```

## 2.2 Scalar equations in 2D

We consider scalar equation of the form

$$ \partial_t u = \Delta[D(x, y)u] + F \left( t, x, y, u, \int_{(-L,L)^2} u \right), \quad t > 0, (x, y) \in (-L, L)^2, $$

with Dirichlet or Neumann boundary condition and initial condition $u_0(x, y)$. Some examples of use are given below.
### Example of use: SKT Model in 2D

$$ \partial_t u = D \Delta u + r(x, y)u - u^2, \quad t > 0, \quad (x, y) \in (-L, L)^2, $$

with Neumann boundary condition.

The following code demonstrates how to solve this problem:

```matlab
# Example of functions
t, x, y, u, rho= symbols('t x y u rho')
fun_u0 = exp(-(x+1)**2 - (y+2)**2)
fun_F = 2 * cos(x)**2 * cos(y)**2 * u - u**2
fun_D = 0.1
```

```matlab
# Times to evaluate and plot
times = np.linspace(0, 50, 100) # times for animation
```

```matlab
# Obtain numerical solution
numericalSolution, xx, yy, _ = sol RD 2D(times, fun_D, fun_F, fun_u0)
```

```matlab
# Animation with the result of the simulation
animplot PDE2D(numericalSolution, xx, yy, times)
```

## 2.3 Systems of two equations in 1D

We consider systems of two equations of the form:

$$ \left\{ \begin{aligned} \partial_t u &= \partial_{xx}[D_u(x)u] + F_u \left( t, x, u, v, \int_{-L}^L u, \int_{-L}^L v \right), & t > 0, x \in (-L, L) \\ \partial_t v &= \partial_{xx}[D_v(x)v] + F_v \left( t, x, u, v, \int_{-L}^L u, \int_{-L}^L v \right), & t > 0, x \in (-L, L) \end{aligned} \right. $$

with Dirichlet or Neumann boundary conditions and initial conditions $u_0(x)$, $v_0(x)$.

### Example of use: System with mutation between two morphs

$$ \left\{ \begin{aligned} \partial_t u(t, x) &= \partial_{xx}u + r_u(x)u - u(u+v) + \mu(v-u), \\ \partial_t v(t, x) &= \partial_{xx}v + r_v(x)v - v(u+v) + \mu(u-v). \end{aligned} \right. $$

The following code demonstrates how to solve this problem:

```matlab
t, x, u, v, rho_u, rho_v = symbols('t x u v rho_u rho_v')
fun_u0 = exp(-(x+1)**2 - (y+2)**2)
fun_v0 = exp(-(x+1)**2 - (y+2)**2)
fun_F_u = (1 + cos(pi * x)) * u - u * (u + v) + 0.1 * (v - u)
fun_F_v = (1 + cos(pi + pi * x)) * v - v * (u + v) + 0.1 * (u - v)
fun_D_u = 0.1
fun_D_v = 0.1
```

```matlab
# Times to evaluate and plot
```
times = np.linspace(0, 50, 100) # times for animation

# Obtain numerical solution

U па, V па, xx, times = sol RD 1D system(times, fun_D_u, fun_D_v, fun_F_u,
fun_F_v, fun_u0, fun_v0)

# Animate the solution

animplot_PDE1D_system(U па, V па, xx, times)

2.4 Systems of two equations in 2D

We consider systems of two equations of the form:

$$
\begin{equation}
\left\{
\begin{aligned}
\partial_t u &= \Delta[D_u(x,y)u] + F_u \left( t, x, y, u, v, \int_{(-L,L)^2} u, \int_{(-L,L)^2} v \right), \\
\partial_t v &= \Delta[D_v(x,y)v] + F_v \left( t, x, y, u, v, \int_{(-L,L)^2} u, \int_{(-L,L)^2} v \right),
\end{aligned}
\right.
\quad t > 0, (x,y) \in (-L,L)^2
\end{equation}
$$

with Dirichlet or Neumann boundary conditions and initial conditions $u_0(x, y)$, $v_0(x, y)$.
An example is given below.

System with mutation between two morphs

$$
\begin{equation}
\left\{
\begin{aligned}
\partial_t u(t, x) &= D \Delta u + r_u(x, y)u - u(u+v) + \mu(v-u), \\
\partial_t v(t, x) &= D \Delta v + r_v(x, y)v - v(u+v) + \mu(u-v).
\end{aligned}
\right.
\end{equation}
$$

The following code shows how to solve this problem:

t, x, y, u, v, rho_u, rho_v = symbols('t x y u v rho_u rho_v')
fun_u0 = exp(-(x+1)**2 - (y+1)**2)
fun_v0 = exp(-(x+1)**2 - (y+1)**2)
fun_F_u = (1 + cos(pi * x)) * (1 + cos(pi * y)) * u
- u * (u + v) + 0.1 * (v - u)
fun_F_v = (1 + cos(pi + pi * x)) * (1 + cos(pi + pi * y)) * v
- v * (u + v) + 0.1 * (u - v)
fun_D_u = 0.1
fun_D_v = 0.1

# Times to evaluate and plot

times = np.linspace(0, 50, 100) # times for animation

# Obtain numerical solution

U па, V па, xx, yy, times = sol RD 2D system(times, fun_D_u, fun_D_v, fun_F_u, fun_F_v, fun_u0, fun_v0)

# Animate the solution

animplot_PDE2D_system(U па, V па, xx, yy, times)
# 3 Reaction-diffusion-advection

We consider here scalar equation of the form

$$ \partial_t u = \partial_{xx}[D(x)u] + \partial_x(v(x)u) + F \left( t, x, u, \int_{-L}^{L} u \right), \quad t > 0, x \in (-L, L) $$

with Dirichlet condition $u(t, \pm L) = 0$ and initial condition $u_0(x)$.

For instance, we solve the SKT model with advection

$$ \partial_t u = \partial_{xx}[D(x)u] + \partial_x(v(x)u) + r(x)u - \gamma(x)u^2, \quad t > 0, x \in (-L, L) $$

with Dirichlet boundary condition with the following code:

```matlab
t, x, u, rho= symbols('t x u rho')
fun_u0 = exp(-(x+0)**2)
fun_F = cos(x)**2 * u - u**2
fun_D = 0.1
fun_v = 1
```

```matlab
times = np.linspace(0, 50, 100) # times for animation
```

```matlab
numericalSolution, xx, _ = sol_RDA_1D(times, fun_D, fun_F, fun_v, fun_u0)
```

```matlab
animplot_PDE1D(numericalSolution, xx, times)
```

# 4 Principal Eigenvalues of elliptic operators

## 4.1 Elliptic operators in 1D

We consider here eigenvalue problems of the form

$$ D(x)\phi'' + q(x)\phi' + r(x)\phi = k\phi, \quad x \in (0,1), $$

with 1-periodic coefficients and periodic boundary conditions: $\phi(0) = \phi(1)$. We look for the principal eigenvalue $k$ associated with a positive eigenfunction $\phi$.

Here is an example of use of the toolbox, for the principal eigenvalue of a Schrödinger operator:

$$ \phi \mapsto \phi'' + r(x)\phi = k\phi, $$

with periodic boundary conditions.

```matlab
x = symbols('x')
fun_D = 1
fun_q = 0
fun_r = (1 + 4 * pi**2 * sp.cos(pi*x)**2 + sp.cos(pi*x)**2 - 2* pi**2)/(1+sp.cos(pi*x)**2)
# with this function r(x), we know that k=1
```
k, PHI = solve_vp1Dfun_D,fun_q,fun_r)

# Print the principal eigenvalue
print("Principal eigenvalue =", k)

## 4.2 Elliptic operators in 2D

We consider the principal eigenvalue problem

$$F(x, y)\Delta\phi + G(x, y)\phi = k\phi,$$

with 1-periodic coefficients and function $\phi$ in both directions.

Using a finite difference approach, with $A$ as the 1D diffusion matrix, we can express this problem as:

$$F \circ [AU + UA] + G \circ U = kU,$$

where $\circ$ denotes the Hadamard (element-wise) product. We then use Kronecker products to compute the eigenvalues:

$$\mathrm{vec}(F \circ [AU + UA] + G \circ U) = \mathrm{diag}(\mathrm{vec}(F))(I \otimes A + A^T \otimes I)\mathrm{vec}(U) + \mathrm{diag}(\mathrm{vec}(G))\mathrm{vec}(U),$$

and define the matrix:

$$M := \mathrm{diag}(\mathrm{vec}(F))(I \otimes A + A^T \otimes I) + \mathrm{diag}(\mathrm{vec}(G)),$$

which allows us to solve the eigenvalue problem:

$$M \mathrm{vec}(U) = k \mathrm{vec}(U).$$

Here is an example of use of the toolbox, for the principal eigenvalue of a Schrödinger operator:

$$\phi \mapsto \Delta\phi + r(x, y)\phi = k\phi,$$

with periodic boundary conditions.

x, y = symbols('x y')

fun_F = 1

fun_G = (20 * cos(x * pi)**2 * cos(y * pi)**4 * pi**2 +
3 * cos(x * pi)**2 * cos(y * pi)**4 -
28 * cos(x * pi)**2 * pi**2 * cos(y * pi)**2 -
2 * cos(y * pi)**4 * pi**2 -
6 * cos(x * pi)**2 * cos(y * pi)**2 +
8 * cos(x * pi)**2 * pi**2 +
4 * cos(y * pi)**2 * pi**2 +
3 * cos(x * pi)**2 -
2 * pi**2 + 3) / (
\cos(x * pi)**2 * \cos(y * pi)**4 -
2 * \cos(x * pi)**2 * \cos(y * pi)**2 +
\cos(x * pi)**2 + 1)

# with this function G(x,y), we know that k=3
k, PHI = solve_vp2Dfun_F,fun_G)

# Print the principal eigenvalue
print("Principal eigenvalue =", k)

## 5 Spreading speeds

We focus on the spreading speeds of the solutions of reaction-diffusion equations of the type

$$\partial_t u = \partial_x(a_L(x)\partial_x u) + b_L(x)\partial_x u + r_L(x)u - u^2, \quad t > 0, x \in \mathbb{R}$$

with L-periodic coefficients. Our numerical codes are based on the Freidlin-Gärtner formula:

$$c_L^* = \inf_{\lambda > 0} \frac{k_\lambda}{\lambda},$$

with $k_\lambda$ the principal eigenvalue of the operator:

$$a_L \phi'' + (a'_L + b_L - 2\lambda a_L)\phi' + (\lambda^2 a_L - a'_L\lambda - b_L\lambda + r_L)\phi = k_\lambda \phi,$$

with L-periodic conditions. Equivalently,

$$\frac{a}{L^2}\phi'' + \frac{1}{L}(\frac{1}{L}a' + b - 2\lambda a)\phi' + (\lambda^2 a - \frac{1}{L}a'\lambda - b\lambda + r)\phi = k_\lambda \phi,$$

with 1-periodic conditions and coefficients a, b, r.

In the example below, we compute the spreading speed of the solution of the SKT equation

$$\partial_t u = \partial_{xx} u + r(x)u - u^2.$$

x = symbols('x')
fun_a = 1
fun_b = 0
fun_r = cos(2 * pi * x)
L = 1
v = spread_speedfun_a, fun_b, fun_r, L)
# Print the principal eigenvalue
print("Speed =", v)

## 6 Usage instructions

To use the toolbox, follow these steps:

* Open the Google Colab notebook at: http://readi2solve.biosp.org

* Sign in with your Google account.

* Run the required cells to set up the environment and solve the PDEs.

* Obtaining accurate results may require adjusting certain default parameters, such as the spatial step size or the size of the domain.

The notebook allows for interactive modifications of the PDE parameters, boundary conditions, and initial conditions.
# 7 Acknowledgments

This work is supported by the ReaCh ANR Project (grant ANR-23-CE40-0023-01). The author thanks the collaborators and users for their feedback and contributions.