[Python] How to python code a distributed delayed differential system in python?

Already tried with ddeint, does not yield good results. I have an equivalent PDE system that works just fine and I can compare results. The key of the problem is that the system is coupled, so I compute the kernel and normalize it before the for loop, however it seems that the contributions of the kernel are such that the expected decay is very much not enough. The code I am using is

import numpy as np
import matplotlib.pyplot as plt
import matplotlib
matplotlib.use('QtAgg') # or 'QtAgg' if available
from tqdm import tqdm



def X_equation(x, y_delta, m, k1, K, q1):
return k1 / (K**m + (y_delta)**m) - q1*x

def Y_equation(x, y, q2, k2):
return k2*x - q2*y

def kernel(tau, D, mu, Delta):
return (1/(np.pi*D*tau)**0.5)*np.exp(-Delta**2 /(4*D*tau) - mu*tau)

k1 = 0.1
K = 1.0
m = 6.0
q1 = 0.03
q2 = 0.03
k2 = 0.1
D = 0.06
mu = 0.0
Delta = 7.5

mu_values = [0.0]

t_span = (0, 1000)
dt = 0.01
t_steps = int(t_span[-1]/dt)


X = np.zeros([t_steps, len(mu_values)])
Y = np.zeros([t_steps, len(mu_values)])
Yprima = np.zeros_like(Y)

time = np.linspace(0,t_span[-1],t_steps)

ker = np.zeros([t_steps, len(mu_values)])

for j in range(0, len(mu_values)):
mu = mu_values[j]
for k in range(1, t_steps): # Start from k=1 to avoid division by zero
tau = k*dt
ker[k, j] = kernel(tau, D, mu, Delta)
ker[:, j] /= np.sum(ker[:, j]) * dt

for j in tqdm(range(0, len(mu_values))):
mu = mu_values[j]
# X[0, j] = (q1 * K**m) / k1
X[0, j] = 1
Y[0, j] = 0
for i in tqdm(range(1,t_steps)):
##Here I compute the convolution integral at each step
precomputed = ker[:i, j] * Y[i-1::-1, j]
cumulative = np.sum(precomputed) * dt
Yprime = cumulative
##This is a RK scheme to compute the rest of the variables once I got the ##convolution integral for this time step
k1_x = X_equation(X[i-1, j], Yprime, m, k1, K, q1)
k1_y = Y_equation(X[i-1, j], Y[i-1, j], k2, q2)

aux_x = X[i-1, j]+dt*k1_x
aux_y = Y[i-1, j]+dt*k1_y

k2_x = X_equation(aux_x, Yprime, m, k1, K, q1)
k2_y = Y_equation(aux_x, aux_y, k2, q2)

X[i, j] = X[i-1, j]+0.5*dt*(k1_x+k2_x)
Y[i, j] = Y[i-1, j]+0.5*dt*(k1_y+k2_y)
Yprima[i, j] = Yprime

plt.plot(time, X[:, 0], label=f"X(t), μ = 0", color='b')
plt.plot(time, Yprima[:, 0], label=f"Ydelta(t), μ = 0", color='g')
plt.plot(time, Y[:,0], label=f"Y0(t), μ = 0", color='r')

plt.xlabel("t")
plt.ylabel("X, Y")
plt.title("Time series")
plt.legend()
plt.grid()

Continue reading...