How to python code a distributed delayed differential system in python?

Question

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()

Break it down into parts. Is your RK scheme working? Can it solve a trivial integration problem correctly? There is a lot of dense impenetrable code above which without knowing the actual differential equation(s) you are trying to solve makes it very difficult to make any meaningful comments. Adding the differential equation would help a lot. — Martin Brown
– Martin Brown, Commented May 27 at 9:25
Have you considered using docs.scipy.org/doc/scipy/tutorial/integrate.html? — Reinderien
– Reinderien, Commented May 27 at 13:01
The equations are ``` \begin{align} \dv{X(t)}{t} &= \frac{k_1}{K^m + \qty(\int_0^t \dd \tau G(\tau) Y(t-\tau))^m} - q_1X(t); \\ \dv{Y(t)}{t} &= k_2 X(t) - q_2 Y(t), \end{align} ``` The RK scheme has been tested for multiple other equations and it yields good results. I am guessing maybe the convolution integral is the troubling stuff, since the scheme is working, and probably has something to do with the distributed delay part of the equations. Most of the code is just for context — Alvaro Ledo Antunez
– Alvaro Ledo Antunez, Commented May 27 at 14:24
Please add it to the question and in a form that is easily legible for a human. Latex markup text escapes is not all that easy to read. A rendering of the equations as an image would be better. BTW Where does the convolution kernel fit into the scheme? — Martin Brown
– Martin Brown, Commented May 27 at 17:19

fiktor · Accepted Answer · 2025-05-28 00:35:09Z

Overall, your code appears to be correct. Since you don't describe the problem fully, it is not clear to see if your code matches the equations you are trying to solve. Below are a observations about the code (some are minor).

`ker` can be computed using vectorized operations

Generally, if you are using numpy, vectorized operations are preferred to for loops. Your for nested for loop

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

Can be written in a vectorized form as follows:

ker[1:, :] = kernel(
    np.arange(1, t_steps)[:, None] * dt,
    D,
    np.array(mu_values)[None, :],
    Delta)
ker /= np.sum(ker, axis=0, keepaxes=True) * dt

Normalization makes kernel dependent on `t_span`

You normalize the kernel s.t.

$int {0}^{T} G(tau) d tau = 1$

(up to errors introduced by numerical integration)

Is this intentional? This clearly depends on T (i.e. t_span[-1] in your code). E.g. if you change T from 1000 (as in your current code) to, e.g., 1500, then the part of your solution for t < 1000 will change because normalization of the kernel will change. Another reason this does not look intentional is that you have (1/(np.pi*D)**0.5) coefficient in your code which is later removed by the normalization.

Mix of second order Runge-Kutta (RK2) and first order integration

You are attempting to implement the second order Runge-Kutta method which may allow higher precision, but the errors you are getting are still first order in dt for the following reasons:

Your integration method is first order.
You do not recompute Yprime when computing k2_x. k2_x is supposed to be the expected derivative of X at time t_i (contrary to k1_x, which is the derivative at time t_i - dt), which means that the integrals in k1_x and k2_x are supposed to be different.

`time` is slightly inconsistent with `dt`

You write

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

You probably meant

time = np.arange(t_steps) * dt

If not, then the time axis is slightly incorrect (a relative error of the order of 1 / t_steps). This is probably not noticeable given t_steps is around 10^5 in your code.

Collectives™ on Stack Overflow

How to python code a distributed delayed differential system in python?

1 Answer 1

`ker` can be computed using vectorized operations

Normalization makes kernel dependent on `t_span`

Mix of second order Runge-Kutta (RK2) and first order integration

`time` is slightly inconsistent with `dt`

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

ker can be computed using vectorized operations

Normalization makes kernel dependent on t_span

Mix of second order Runge-Kutta (RK2) and first order integration

time is slightly inconsistent with dt

Comments

Your Answer

Sign up or log in

Post as a guest

Related

`ker` can be computed using vectorized operations

Normalization makes kernel dependent on `t_span`

`time` is slightly inconsistent with `dt`