I’d like here to share the reoprt I prepared for my final project for the course of Numerical Analysis of PDE here at KAUST. All the code to reproduce the image here presented can be found on this dedicated repo. We are interested in solving numerically the PDE system proposed by Alan Turing in 1952 to model the pattern formation on the tissue of some animals such as zebras, ladybugs etc. In particular we will focus our attention on the following diffusion-reaction coupled PDE system:

(1)

and we can separate the diffusion terms from the reactions term, obtaining,

(2)

where the functions and are clearly defined as:

(3)

To begin we will focus our attention on solving the above coupled system of autonomous ordinary differential equations representing the reaction part of the Turing model. In particular we will solve these implementing the 4th order Runge-Kutta method. First we investigate numerically the stability of the method, to do this we kept fixed the initial conditions , and vary the number of points in the temporal discretization for the values of the parameters and in the first line of Parameters Table. In particular from Figure 1 one observe that the method appear to be stable if we chose more then points in the temporal discretization. Now we investigate the long term behaviour of the quantitates and for different values in Parameters Table. As we can observe from Figure 2, for different values of initial data and of the parameters all the ODEs in the long term engage in an alternating oscillatory behaviour, similar to the one observed in simpler non linear ODE model such as Lotka-Volterra equations. Studying the linearization of the systems in the equilibrium points, one will found that the Jacobian has only negative eigenvalues. In particular numerically we notice that all the equilibrium points of the system are unstable.

0.00225 | 0.0045 | 0.02 | 0.2 | 0.899 | -0.91 | -0.899 |

0.001 | 0.0045 | 0.02 | 0.2 | 0.899 | -0.91 | -0.899 |

0.00225 | 0.0045 | 0.02 | 0.2 | 1.9 | -0.91 | -1.9 |

0.00225 | 0.0045 | 2.02 | 0.0 | 2.0 | -2.0 | -2.0 |

0.00225 | 0.0021 | 3.5 | 0.0 | 0.899 | -0.91 | -0.899 |

0.00225 | 0.0045 | 0.02 | 0.2 | 1.9 | -0.85 | -1.9 |

0.00225 | 0.0045 | 2.02 | 0.0 | 2.0 | -0.91 | -2.0 |

**Parameters Table**The table contains the parameters for which we simulate the Turing patterns model.

**Figure 2**The following figures show the behaviour of the reaction ODE \eqref{eq:ReactionODE}, for different initial values and different choice of the parameters and . The number of points in the temporal discretization was fixed equal to . In particular from right to left we have initial values equal to , and . Each row correspond to the equally number row in Parameters Table.

The next step will be to solve the diffusion equation, in particular we implemented a second order five points finite difference scheme to approximate the second derivative. In particular the matrix:

(4)

approximate the one dimensional second derivative with periodic boundary condition, then the approximation of the Laplacian is obtain by the Kronecker product , since the same mesh is taken in the horizontal and vertical direction. Now we would like to compute theoretically how small to chose the time marching step in order to obtain a stable integration in time. In order to answer this question we notice that is a circulant matrix and therefore its eigenvalues can be easily computed,

(5)

Furthermore since circulant matrix commute one has that the eigenvalues of the Kronecker product are given by:

(6)

and this tell us that the eigenvalue of lie on the real line and attains as maximum value and as minimum value . Now for the numerical method to be stable we wont , the marching step in time, to be such that lies with in the stability region for the 4-th order Runge-Kutta method. Now from Figure 3 we notice that this imply that the following inequality has to be verified,

(7)

where is the where the negative real line meets the border of the stability region, marked in red in Figure 3. If we take we get that for the first value in the Table ?? that we should take a time marching step . The simulation of the diffusive phenomena obtain with this numerical scheme can be observed in Figure 4.

**Figure 3**Stability region for the 4-th order Runge-Kutta method in red and the eigenvalue for the matrix marked in blue, when we use points in the spatial mesh and .

**Figure 4**The fist row in the figures show the initial data that were taken for all the simulation. The second row show the result of the diffusion only model at . The last row show the result of the diffusion only model at . Respectively from right to left we have the values for and .

Now is time to simulate both the reaction and diffusion term together, the result of this simulation can be seen in Figure 5. Using the 4th order Runge-Kutta method for the running time is approximatively 3 minutes.

**Figure 5**The first row show the result of the reaction diffusion model at . The last row show the result of the reaction diffusion model at . Respectively from right to left we have the values for and . As parameters for the Turing model we chose thre one in the first line of Parameters Table.

To over come the CFL condition and have a faster method we resort to a splitting technique, i.e. we solve first the diffusion part of the PDE using an implicit Euler method. The implicit Euler scheme is absolutely stable therefore we can use much larger time steps. Then we use the solution obtain as initial data to the reaction part of the system and this is solved using the 4-th order Runge-Kutta method previously implemented. This is advantageous because as we have seen in the beginning the time step required to solve in a stable manner \eqref{eq:ReactionODE} is much larger then the one needed to solve in a stable manner the diffusive part of the model. The result of this simulation at , obtained using a time marching step are shown in Figure 6. From Figure 6 we observe that even starting from initial data without any sort of pattern such as the one shown in Figure 4 one obtains solutions with patterns such as the one presented in different animals tissue.

**Figure 6**Each couple shows the numerical integration obtained with 100 points in the spatial mesh grid and a time step , for the Turing model with parameter corresponding to the same row in Parameters Table. Respectively from right to left we have the values for and .

Here will follow the code to obtain all figures presented in this report.

```
import numpy as np
from math import sqrt
import scipy.integrate
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from IPython.display import HTML
from scipy import sparse
from scipy.sparse.linalg import eigs
from scikits.umfpack import spsolve #UMFPACK If you work in serial.
from tqdm import trange, tqdm
from nodepy.runge_kutta_method import *
class ODESol:
def __init__(self,timesteps,timestep,U):
self.t = timesteps;
self.h = timestep;
self.y = U;
def RK4(F,T,U0,arg,N):
tt = np.linspace(T[0],T[1],N);
h = tt[1]-tt[0];
U = np.zeros([len(U0),N]);
U[:,0] = U0;
for i in trange(0,N-1):
Y1 = U[:,i];
Y2 = U[:,i] + 0.5*h*F(tt[i],Y1,arg);
Y3 = U[:,i] + 0.5*h*F(tt[i]+0.5*h,Y2,arg);
Y4 = U[:,i] + h*F(tt[i]+0.5*h,Y3,arg);
U[:,i+1] = U[:,i]+(h/6)*(F(tt[i],Y1,arg)+2*F(tt[i]+ 0.5*h,Y2,arg)+2*F(tt[i]+ 0.5*h,Y3,arg)+F(tt[i]+h,Y4,arg))
sol = ODESol(tt,h,U);
return sol;
Table = []
Table = Table + [{"delta1":0.00225,"delta2":0.0045,"tau1":0.02,"tau2":0.2,"alpha":0.899, "beta":-0.91,"gamma":-0.899}]
Table = Table + [{"delta1":0.001,"delta2":0.0045,"tau1":0.02,"tau2":0.2,"alpha":0.899, "beta":-0.91,"gamma":-0.899}]
Table = Table + [{"delta1":0.00225,"delta2":0.0045,"tau1":0.02,"tau2":0.2,"alpha":1.9, "beta":-0.91,"gamma":-1.9}]
Table = Table + [{"delta1":0.00225,"delta2":0.0045,"tau1":2.02,"tau2":0.0,"alpha":2.0, "beta":-0.91,"gamma":-2}]
Table = Table + [{"delta1":0.00105,"delta2":0.0021,"tau1":3.5,"tau2":0.0,"alpha":0.899, "beta":-0.91,"gamma":-0.899}]
Table = Table + [{"delta1":0.00225,"delta2":0.0045,"tau1":0.02,"tau2":0.2,"alpha":1.9, "beta":-0.85,"gamma":-1.9}]
Table = Table + [{"delta1":0.00225,"delta2":0.0005,"tau1":2.02,"tau2":0.0,"alpha":2.0, "beta":-0.91,"gamma":-2}]
def Reaction(t,x,parameters):
#The ODE is autonomus so we don't really need
#the depende on time.
u = x[0]; v = x[1]; #We grab the useful quantity.
d1 = parameters["delta1"]; d2 = parameters["delta2"];
t1 = parameters["tau1"]; t2 = parameters["tau2"];
a = parameters["alpha"]; b = parameters["beta"]; g = parameters["gamma"];
du = a*u*(1-t1*(v**2))+v*(1-t2*u);
dv = b*v*(1+(a*t1/b)*(u*v))+u*(g+t2*v);
return np.array([du,dv])
t0 = 0. # Initial time
u0 = np.array([0.5,0.5])# Initial values
tfinal = 100. # Final time
dt_output=0.1# Interval between output for plotting
N=int(tfinal/dt_output) # Number of output times
print(N)
tt=np.linspace(t0,tfinal,N) # Output times
ODE = RK4(Reaction,[t0,tfinal],u0,Table[6],N);
uu=ODE.y
plt.plot(tt,uu[0,:],tt,uu[1,:])
plt.legend(["u","v"])
plt.show();
```

```
import numpy as np
from math import sqrt
import scipy.integrate
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from IPython.display import HTML
from scipy import sparse
from scipy.sparse.linalg import eigs
from scikits.umfpack import spsolve #UMFPACK If you work in serial.
from tqdm import trange, tqdm
from nodepy.runge_kutta_method import *
class ODESol:
def __init__(self,timesteps,timestep,U):
self.t = timesteps;
self.h = timestep;
self.y = U;
def RK4(F,T,U0,arg,N):
tt = np.linspace(T[0],T[1],N);
h = tt[1]-tt[0];
U = np.zeros([len(U0),N]);
U[:,0] = U0;
for i in trange(0,N-1):
Y1 = U[:,i];
Y2 = U[:,i] + 0.5*h*F(tt[i],Y1,arg);
Y3 = U[:,i] + 0.5*h*F(tt[i]+0.5*h,Y2,arg);
Y4 = U[:,i] + h*F(tt[i]+0.5*h,Y3,arg);
U[:,i+1] = U[:,i]+(h/6)*(F(tt[i],Y1,arg)+2*F(tt[i]+ 0.5*h,Y2,arg)+2*F(tt[i]+ 0.5*h,Y3,arg)+F(tt[i]+h,Y4,arg))
sol = ODESol(tt,h,U);
return sol;
Table = []
Table = Table + [{"delta1":0.00225,"delta2":0.0045,"tau1":0.02, "tau2":0.2,"alpha":0.899,"beta":-0.91,"gamma":-0.899}]
Table = Table + [{"delta1":0.001,"delta2":0.0045,"tau1":0.02, "tau2":0.2,"alpha":0.899,"beta":-0.91,"gamma":-0.899}]
Table = Table + [{"delta1":0.00225,"delta2":0.0045,"tau1":0.02, "tau2":0.2,"alpha":1.9,"beta":-0.91,"gamma":-1.9}]
Table = Table + [{"delta1":0.00225,"delta2":0.0045,"tau1":2.02, "tau2":0.0,"alpha":2.0,"beta":-0.91,"gamma":-2}]
Table = Table + [{"delta1":0.00105,"delta2":0.0021,"tau1":3.5, "tau2":0.0,"alpha":0.899,"beta":-0.91,"gamma":-0.899}]
Table = Table + [{"delta1":0.00225,"delta2":0.0045,"tau1":0.02, "tau2":0.2,"alpha":1.9,"beta":-0.85,"gamma":-1.9}]
Table = Table + [{"delta1":0.00225,"delta2":0.0005,"tau1":2.02, "tau2":0.0,"alpha":2.0,"beta":-0.91,"gamma":-2}]
def laplacian_1D(m):
em = np.ones(m)
e1=np.ones(m-1)
A = (sparse.diags(-2*em,0)+sparse.diags(e1,-1)+sparse.diags(e1,1))/((2/(m+1))**2);
A[0,-1]=1/((2/(m+1))**2);
A[-1,0]=1/((2/(m+1))**2);
return A;
def laplacian_2D(m):
I = np.eye(m)
A = laplacian_1D(m)
return sparse.kron(A,I) + sparse.kron(I,A)
def Diffusion(t,x,parameters):
B = sparse.bmat([[parameters["delta1"]*A,None],[None,parameters["delta2"]*A]]);
return B*x;
m=100
x=np.linspace(-1,1,m+2); x=x[1:-1]
y=np.linspace(-1,1,m+2); y=y[1:-1]
X,Y=np.meshgrid(x,y)
A=laplacian_2D(m)
#Generating the initial data
mu, sigma = 0, 0.5 # mean and standard deviation
u0 = np.random.normal(mu, sigma,m**2);
v0 = np.random.normal(mu, sigma,m**2)
#Plotting the initial data
plt.figure()
U0=u0.reshape([m,m])
plt.pcolor(X,Y,U0)
plt.colorbar();
plt.figure()
V0=v0.reshape([m,m])
plt.pcolor(X,Y,V0)
plt.colorbar();
TIndex = 0;
t0 = 0.0 # Initial time
tfinal = 100 # Final time
dt_output=0.018 # Interval between output for plotting
N=int(tfinal/dt_output) # Number of output times
print("CFL suggested time step",(2.5*(0.0198)**2)/(4*(Table[TIndex]["delta1"]+Table[TIndex]["delta2"])));
#ODE = scipy.integrate.solve_ivp(Diffusion,[t0,tfinal],np.append(u0,v0,axis=0),args=[Table[TIndex]],method='RK45',t_eval=tt,atol=1.e-10,rtol=1.e-10);
ODE = RK4(Diffusion,[t0,tfinal],np.append(u0,v0,axis=0),Table[TIndex],N);
uu=ODE.y
print("Step size of the Method",ODE.h);
print(ODE.h*eigs(sparse.bmat([[Table[TIndex]["delta1"]*A,None],[None,Table[TIndex]["delta2"]*A]]),k=10)[0])
RK44 = loadRKM('RK44')
RK44.plot_stability_region(bounds=[-5,1,-5,5])
plt.plot(ODE.h*eigs(sparse.bmat([[Table[TIndex]["delta1"]*A,None],[None,Table[TIndex]["delta2"]*A]]),k=10)[0].real,ODE.h*eigs(sparse.bmat([[Table[TIndex]["delta1"]*A,None],[None,Table[TIndex]["delta2"]*A]]),k=10)[0].imag,"b*")
ut = uu[0:m**2,-1];
vt = uu[m**2:,-1];
Ut=ut.reshape([m,m])
plt.figure()
plt.pcolor(X,Y,Ut)
plt.colorbar();
plt.figure()
Vt=vt.reshape([m,m])
plt.pcolor(X,Y,Vt)
plt.colorbar();
plt.show();
```

```
import numpy as np
from math import sqrt
import scipy.integrate
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from IPython.display import HTML
from scipy import sparse
from scipy.sparse.linalg import eigs
from scikits.umfpack import spsolve #UMFPACK If you work in serial.
from tqdm import trange, tqdm
from nodepy.runge_kutta_method import *
class ODESol:
def __init__(self,timesteps,timestep,U):
self.t = timesteps;
self.h = timestep;
self.y = U;
def laplacian_1D(m):
em = np.ones(m)
e1=np.ones(m-1)
A = (sparse.diags(-2*em,0)+sparse.diags(e1,-1)+sparse.diags(e1,1))/((2/(m+1))**2);
A[0,-1]=1/((2/(m+1))**2);
A[-1,0]=1/((2/(m+1))**2);
return A;
def laplacian_2D(m):
I = np.eye(m)
A = laplacian_1D(m)
return sparse.kron(A,I) + sparse.kron(I,A)
def Reaction(t,x,parameters):
u = x[0:m**2]; v = x[m**2:]; #We grab the useful quantity.
d1 = parameters["delta1"]; d2 = parameters["delta2"];
t1 = parameters["tau1"]; t2 = parameters["tau2"];
a = parameters["alpha"]; b = parameters["beta"]; g = parameters["gamma"];
#Reaction
du = a*u*(1-t1*(v**2))+v*(1-t2*u);
dv = b*v*(1+(a*t1/b)*(u*v))+u*(g+t2*v);
b = np.append(du,dv,axis=0);
return b;
def SplitSolver(F,T,U0,arg,N):
tt = np.linspace(T[0],T[1],N);
h = tt[1]-tt[0];
U = np.zeros([len(U0),N]);
U[:,0] = U0;
for i in trange(0,N-1):
B = sparse.bmat([[arg["delta1"]*A,None],[None,arg["delta2"]*A]]);
UStar = spsolve((sparse.identity(B.shape[0])-h*B),U[:,i])
Y1 = UStar;
Y2 = UStar + 0.5*h*F(tt[i],Y1,arg);
Y3 = UStar + 0.5*h*F(tt[i]+0.5*h,Y2,arg);
Y4 = UStar + h*F(tt[i]+0.5*h,Y3,arg);
U[:,i+1] = UStar+(h/6)*(F(tt[i],Y1,arg)+2*F(tt[i]+ 0.5*h,Y2,arg)+2*F(tt[i]+ 0.5*h,Y3,arg)+F(tt[i]+h,Y4,arg))
sol = ODESol(tt,h,U);
return sol;
Table = []
Table = Table + [{"delta1":0.00225,"delta2":0.0045,"tau1":0.02, "tau2":0.2,"alpha":0.899,"beta":-0.91,"gamma":-0.899}]
Table = Table + [{"delta1":0.001,"delta2":0.0045,"tau1":0.02, "tau2":0.2,"alpha":0.899,"beta":-0.91,"gamma":-0.899}]
Table = Table + [{"delta1":0.00225,"delta2":0.0045,"tau1":0.02, "tau2":0.2,"alpha":1.9,"beta":-0.91,"gamma":-1.9}]
Table = Table + [{"delta1":0.00225,"delta2":0.0045,"tau1":2.02, "tau2":0.0,"alpha":2.0,"beta":-0.91,"gamma":-2}]
Table = Table + [{"delta1":0.00105,"delta2":0.0021,"tau1":3.5, "tau2":0.0,"alpha":0.899,"beta":-0.91,"gamma":-0.899}]
Table = Table + [{"delta1":0.00225,"delta2":0.0045,"tau1":0.02, "tau2":0.2,"alpha":1.9,"beta":-0.85,"gamma":-1.9}]
Table = Table + [{"delta1":0.00225,"delta2":0.0005,"tau1":2.02, "tau2":0.0,"alpha":2.0,"beta":-0.91,"gamma":-2}]
TIndex = int(input("What pattern you want to simulate ? "))
m=100
x=np.linspace(0,1,m+2); x=x[1:-1]
y=np.linspace(0,1,m+2); y=y[1:-1]
X,Y=np.meshgrid(x,y)
A=laplacian_2D(m)
#Generating the initial data
mu, sigma = 0, 0.5 # mean and standard deviation
u0 = np.random.normal(mu, sigma,m**2);
v0 = np.random.normal(mu, sigma,m**2)
#Plotting the initial data
plt.figure()
U0=u0.reshape([m,m])
plt.pcolor(X,Y,U0)
plt.colorbar();
plt.title("Initial Data u");
plt.figure()
V0=v0.reshape([m,m])
plt.pcolor(X,Y,V0)
plt.colorbar();
plt.title("Initial Data v");
t0 = 0.0 # Initial time
tfinal = 150 # Final time
dt_output=0.3 # Interval between output for plotting
N=int(tfinal/dt_output) # Number of output times
ODE = SplitSolver(Reaction,[t0,tfinal],np.append(u0,v0,axis=0),Table[TIndex],N);
uu=ODE.y
ut = uu[0:m**2,-1];
vt = uu[m**2:,-1];
Ut=ut.reshape([m,m])
plt.figure()
plt.pcolor(X,Y,Ut)
plt.colorbar();
plt.figure()
Vt=vt.reshape([m,m])
plt.pcolor(X,Y,Vt)
plt.colorbar();
plt.show();
```