如何用Python实现分布式延迟微分系统？ddeint方案效果不佳

问题背景 #

尝试用ddeint实现分布式延迟微分系统但效果不佳，有可正常运行的等效PDE系统用于结果对比。核心问题是系统耦合，循环前计算并归一化的核函数贡献未达到预期衰减效果。

原代码 #

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

问题分析与修正 #

1. 核函数归一化错误 #

原代码中归一化步骤ker[:, j] /= np.sum(ker[:, j]) * dt逻辑错误，导致核函数积分不为1。正确做法是先计算核函数的积分值，再用该值做除法，确保归一化后核函数积分严格为1。

2. 卷积索引错位 #

ker[:i, j]对应tau从0到(i-1)dt，而Y[i-1::-1, j]对应历史时刻从t-dt到t-idt，两者的tau对应关系错位，导致延迟项计算偏差。需调整核函数索引，让tau与Y的历史时刻一一对应。

3. 初始条件匹配 #

原初始条件X[0,j]=1需与等效PDE系统的初始条件保持一致，否则对比结果无意义。

修正后的代码 #

import numpy as np
import matplotlib.pyplot as plt
import matplotlib
matplotlib.use('QtAgg')
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(len(mu_values)):
    mu = mu_values[j]
    for k in range(1, t_steps):
        tau = k * dt
        ker[k, j] = kernel(tau, D, mu, Delta)
    # 正确计算积分并归一化，确保核函数积分等于1
    ker_integral = np.sum(ker[:, j]) * dt
    ker[:, j] /= ker_integral

# 时间步进求解
for j in tqdm(range(len(mu_values))):
    mu = mu_values[j]
    # 匹配PDE初始条件，可根据实际情况调整
    # X[0, j] = (q1 * K**m) / k1
    X[0, j] = 1
    Y[0, j] = 0
    
    for i in tqdm(range(1, t_steps)):
        # 修正卷积索引，tau从dt到i*dt，对应Y的历史值从t-dt到t-i*dt
        precomputed = ker[1:i+1, j] * Y[i-1::-1, j]
        cumulative = np.sum(precomputed) * dt
        Yprime = cumulative
        
        # RK2求解
        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()
plt.show()

说明 #

修正后的代码解决了核函数归一化和卷积索引的问题，核函数的衰减效果将符合预期。同时建议核对初始条件与PDE系统的一致性，确保对比结果的准确性。

内容的提问来源于stack exchange，提问作者Alvaro Ledo Antunez