如何高效数值求解由延迟微分方程定义的Buchstab函数

明确方程定义 #

Buchstab函数由以下延迟微分方程（DDE）定义：

• 初始条件：当 (1 < x \leq 2) 时，(F(x) = 1)
• 微分关系：当 (x > 2) 时，(\frac{dF}{dx} = \frac{F(x-1)}{x})

核心计算思路：分段递推+数值积分 #

由于方程是延迟型，当前区间的函数值仅依赖前一个整数区间的结果，因此可以通过分段递推的方式高效计算：
对于任意整数 (n \geq 2)，区间 (n < x \leq n+1) 内的 (F(x)) 可通过对微分方程积分得到：
[F(x) = F(n) + \int_{n}^{x} \frac{F(t-1)}{t} dt]
其中 (F(n)) 是前一区间的终点值，(F(t-1)) 是已计算完成的 (n-1 < t-1 \leq n) 区间内的函数值。

高效实现步骤 #

• 区间划分：将计算范围按整数点划分为连续区间 ([2,3], [3,4], ..., [N, N+1])（直到目标上限 (X_{max})）
• 初始区间初始化：直接设置 (1 < x \leq 2) 时 (F(x)=1)
• 递推计算：
  1. 对当前区间的每个点，通过插值（线性/三次样条）获取前一区间的 (F(x-1)) 值
  2. 用数值积分方法（梯形法、辛普森法）计算积分项，得到当前区间的 (F(x))
• 内存优化：仅保留前一个区间的函数值，无需存储所有历史数据，大幅降低内存占用

优化技巧 #

• 自适应步长：在函数变化较明显的区间（如 (2 < x < 5)）使用小步长，平缓区间（(x > 5)）增大步长，平衡精度与计算速度
• 解析积分结合拟合：对每个区间的 (F(x)) 进行多项式拟合（如三次样条），将积分转化为解析计算，比纯数值积分更快
• 渐近值验证：当 (x \to +\infty) 时，(F(x) \to e^{-\gamma})（(\gamma \approx 0.5671)，欧拉常数），可用于验证计算结果的合理性

代码示例（Python） #

import numpy as np
import matplotlib.pyplot as plt

def compute_buchstab(x_max, step=0.01):
    # 初始化初始区间(1,2]的函数值
    x_vals = np.arange(1, 2 + step, step)
    f_vals = np.ones_like(x_vals)
    
    n = 2
    while n < x_max:
        # 当前区间的x序列
        current_x = np.arange(n, min(n+1, x_max) + step, step)
        # 对应x-1的序列，属于前一个区间
        prev_x = current_x - 1
        # 线性插值获取F(x-1)
        prev_f = np.interp(prev_x, x_vals, f_vals)
        # 计算微分方程右端项
        df_dx = prev_f / current_x
        
        # 梯形法积分计算当前区间的F(x)
        current_f = np.zeros_like(current_x)
        current_f[0] = f_vals[-1]  # 初始值为前一区间的终点值
        for i in range(1, len(current_x)):
            dx = current_x[i] - current_x[i-1]
            current_f[i] = current_f[i-1] + (df_dx[i-1] + df_dx[i]) * dx / 2
        
        # 更新存储的序列（仅保留前一区间和当前区间，节省内存）
        x_vals = np.concatenate([x_vals[-int(1/step):], current_x[1:]])
        f_vals = np.concatenate([f_vals[-int(1/step):], current_f[1:]])
        
        n += 1
    
    return x_vals, f_vals

# 计算到x=10，步长0.01
x, f = compute_buchstab(10, 0.01)

# 绘图展示
plt.plot(x, f)
plt.xlabel('x')
plt.ylabel('F(x)')
plt.title('Buchstab Function')
plt.hlines(np.exp(-np.euler_gamma), xmin=1, xmax=10, linestyles='--', colors='r', label=r'$e^{-\gamma}$')
plt.legend()
plt.show()

内容的提问来源于stack exchange，提问作者Simd