如何结合Danckwerts边界条件，用Scipy odeint实现线法(MOL)求解

问题描述 #

我正在求解管式反应器的对流-扩散-反应问题，采用线法（Method of Lines, MOL）——通过空间离散将问题转化为时间域的常微分方程组来求解。目前最大的难点是处理Danckwerts边界条件，我已经通过二阶向前近似推导得到C₁的表达式，通过二阶向后近似推导得到C_N+1的表达式。

但我写的代码存在问题：根据Scipy官方文档，odeint函数中不能直接修改输入的y数组，仅能修改其导数。现有代码直接修改了I和mu0的边界点值，不符合要求，求正确的实现方法。

现有代码 #

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator
import time

# 参数定义
L = 1.0           # 反应器长度（米）
v_x = 3.3e-4      # 流体平均流速（m/s）
D = 0.0000033     # 扩散系数（m^2/s）
Pe = L*v_x/D      # 佩克莱数

kd = 0.000055         # 解离常数（1/s）
factor = 0.6          # 效率因子
Da = L*kd*factor/v_x  # 达姆科勒数

kt = 6e7
Da_mu0 = L*kt/v_x

# 网格与时间步设置
m = 101
n = 100001
x = np.linspace(0, 1, m) # 空间网格
t = np.linspace(0, 1, n) # 时间网格
h = x[1] - x[0] # 空间步长
k = t[1] - t[0] # 时间步长

from scipy.integrate import odeint

# 初始条件
init_I = np.ones([m])
init_mu0 = np.zeros([m])

def odefunc(y,t):

    I, mu0 = np.split(y,2)

    dIdt = np.zeros(x.shape)
    dmu0dt = np.zeros(x.shape)

    # 直接修改边界点值，此处不符合odeint要求
    I[0] = (2*Pe/m + 4*I[1] - I[2])/(2*Pe/m + 3)
    I[-1] = (4*I[-2] - I[-3])/3

    mu0[0] = (4*mu0[1] - mu0[2])/(3 + 2*Pe/m)
    mu0[-1] = (4*mu0[-2] - mu0[-3])/3

    # 内部点的导数计算
    for i in range(1, m-1):
        dIdt[i] = -(I[i+1] - I[i-1])/(2/m) + (1/Pe)*(I[i+1] - 2*I[i] + I[i-1])/(1/m**2) - Da*(I[i+1] + I[i-1])/2
        dmu0dt[i] = -(mu0[i+1] - mu0[i-1])/(2/m) + (1/Pe)*(mu0[i+1] - 2*mu0[i] + mu0[i-1])/(1/m**2) + Da*(I[i+1] + I[i-1]) - (Da_mu0/4)*(mu0[i+1] +mu0[i-1])**2

    return np.ravel([dIdt, dmu0dt])

# 初始条件展开
init = np.ravel([init_I, init_mu0])
# 求解ODE
sol = odeint(odefunc, init, t, atol=1e-10, rtol=1e-10)
sol_I = sol[:,:m]
sol_mu0 = sol[:,m:]

# 绘制I的3D结果图
fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
M, N = np.meshgrid(x, t)

surf = ax.plot_surface(M, N, sol_I, cmap=cm.viridis, linewidth=0, antialiased=True)

ax.set_zlim(0, 1.0)
ax.zaxis.set_major_locator(LinearLocator(5))
ax.set_xlabel('x\'')
ax.set_ylabel('t\'')
ax.set_zlabel('c\'')

ax.zaxis.set_major_formatter('{x:.02f}')
ax.set_title(f'{Pe =:.02f} and {Da =:.02f}')

fig.colorbar(surf, shrink=0.5, aspect=5)

plt.show()

# 绘制mu0的3D结果图
fig, ax = plt.subplots(subplot_kw={"projection": "3d"})

surf = ax.plot_surface(M, N, sol_mu0, cmap=cm.viridis, linewidth=0, antialiased=True)

ax.set_zlim(0, 1.0)
ax.zaxis.set_major_locator(LinearLocator(5))
ax.set_xlabel('x\'')
ax.set_ylabel('t\'')
ax.set_zlabel('c\'')

ax.zaxis.set_major_formatter('{x:.02f}')
ax.set_title(f'{Pe =:.02f} and {Da =:.02f}')

fig.colorbar(surf, shrink=0.5, aspect=5)

plt.show()

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