如何在scikits.odes.dae求解器中集成中心差分法求解时空耦合方程

问题背景 #

需要求解的时空耦合DAE方程组：
$$
\begin{cases}
\frac{dy_1}{dt} = \frac{dy_1}{dz} + y_2 \
y_1 = 5 y_2
\end{cases}
$$

当前代码中直接传入常数dydz，但实际需要用中心差分格式计算空间导数：$\frac{dy_1}{dz}[i] = \frac{y_1[i+1] - y_1[i-1]}{2dz}$，需将该逻辑集成到DAE求解器中。

原始代码：

import matplotlib.pyplot as plt
import numpy as np
from scikits.odes import dae

N  = 51 #number of spacesteps
L  = 1.0 #[m] length of sorbent bed, also a guess 
dz = L/(N-1) #[m] length of space step

time = np.arange(0, 1.5, 0.1)
dydz = 1

y0 = [1, 0.2] #initial values y0[0] = y1 and y0[1] = y2
yp0 = [1, 1] #initial guess for \dot{y1} and \dot{y2}

def trial_space(t, y, ydot, result):
    result[0] = ydot[0] - 6 * dydz + y[1]
    result[1] = y[0] - 5 * y[1]

solver = dae('ida', trial_space)
solution = solver.solve(time, y0, yp0)

解决方案 #

核心思路是：将所有空间点的状态变量打包成一个一维向量，在DAE残差函数中针对每个空间点计算对应的空间导数（内部点用中心差分，边界点用一阶差分处理），再代入残差方程。

修改后的完整代码 #

import matplotlib.pyplot as plt
import numpy as np
from scikits.odes import dae

N  = 51  # 空间步数
L  = 1.0  # 床层长度
dz = L/(N-1)  # 空间步长

time = np.arange(0, 1.5, 0.1)

# 初始化状态向量：每个空间点包含y1和y2，共2*N个元素
# 假设初始时所有空间点y1=1，y2=0.2（满足y1=5y2）
y0 = np.tile([1.0, 0.2], N)
# 初始导数猜测：所有变量的导数设为1.0
yp0 = np.ones_like(y0)

def dae_residual(t, y, ydot, result):
    # 拆分状态变量：y1的所有空间点，y2的所有空间点
    y1 = y[0::2]
    y2 = y[1::2]
    # 拆分导数变量：dy1/dt的所有空间点，dy2/dt的所有空间点
    dy1_dt = ydot[0::2]
    dy2_dt = ydot[1::2]
    
    # 计算每个空间点的dy1/dz
    dy1_dz = np.zeros_like(y1)
    # 内部点：中心差分
    dy1_dz[1:-1] = (y1[2:] - y1[:-2]) / (2 * dz)
    # 左边界：向前差分（一阶）
    dy1_dz[0] = (y1[1] - y1[0]) / dz
    # 右边界：向后差分（一阶）
    dy1_dz[-1] = (y1[-1] - y1[-2]) / dz
    
    # 填充残差结果
    for i in range(N):
        # 第一个方程残差：dy1/dt - dy1/dz - y2 = 0
        result[2*i] = dy1_dt[i] - dy1_dz[i] - y2[i]
        # 第二个方程残差：y1 -5y2 =0
        result[2*i +1] = y1[i] - 5 * y2[i]

# 初始化DAE求解器
solver = dae('ida', dae_residual)
# 求解
solution = solver.solve(time, y0, yp0)

# 可选：绘制结果，比如t=1.0时y1沿空间的分布
if solution.success:
    t_idx = np.where(solution.t == 1.0)[0][0]
    y1_final = solution.y[t_idx][0::2]
    z = np.linspace(0, L, N)
    plt.plot(z, y1_final)
    plt.xlabel('z [m]')
    plt.ylabel('y1')
    plt.title('y1 distribution at t=1.0')
    plt.show()
else:
    print(f"求解失败：{solution.message}")

关键细节说明 #

1. 状态向量打包：将N个空间点的y1和y2依次存入一维数组y，结构为[y1_0, y2_0, y1_1, y2_1, ..., y1_{N-1}, y2_{N-1}]，这样DAE求解器能同时处理所有空间点的时间演化。
2. 空间导数计算：
   • 内部空间点（1 <= i <= N-2）使用中心差分，精度更高；
   • 边界点（i=0和i=N-1）无法使用中心差分（缺少邻点），因此采用一阶向前/向后差分，也可根据需求改用其他边界格式（比如二阶单边差分）。
3. 残差函数填充：对每个空间点，分别计算两个方程的残差，确保每个空间点都满足DAE约束。

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