如何在 Python 中解决这个与时间相关的 PDE？

Question

有没有办法在 Python 中数值求解以下 PDE？

RHS 的第二项对时间和空间都有一个导数。

我尝试在 Python 中使用 Py-PDE 包，它只解决了 dy(x,t)/dt = f(y(x,t)) 的形式，所以我尝试使用类似于 的求根算法scipy fzero 得到 dy(x,t)/dt - f(y(x,t),dy(x,t)/dt) = 0 的解（求解dy(x,t)/dt)。

class nonlinearPDE(pde.PDEBase):
    def __init__(self, bc={"derivative":0}):
        self.bc = bc #boundary conditions for operators
    
    def _make_pde_rhs_numba(self, state):
        """numba-compiled implementation of the PDE"""
        laplace = state.grid.make_operator("laplace", bc=self.bc)
        def findroot(f, df, x0, nmax):
            """Newton–Raphson method"""
            for i in range(nmax):
                x0 = x0 - f(x0)/df(x0)
            return x0
    
        @jit
        def pde_rhs(y, t):
            func = lambda dydt : dydt - a*laplace(y) - b*laplace(dydt)
            dydt = findroot(func, lambda x : 1, 0, 1)
            return dydt
    return pde_rhs

但是，当程序尝试求解 PDE 时，它会抛出错误：

  File "...\anaconda3\lib\site-packages\pde\solvers\controller.py", line 191, in run
    t = stepper(state, t, t_break)

  File "...\anaconda3\lib\site-packages\pde\solvers\scipy.py", line 82, in stepper
    sol = integrate.solve_ivp(

  File "...\anaconda3\lib\site-packages\scipy\integrate\_ivp\ivp.py", line 542, in solve_ivp
    solver = method(fun, t0, y0, tf, vectorized=vectorized, **options)

  File "...\anaconda3\lib\site-packages\scipy\integrate\_ivp\rk.py", line 94, in __init__
    self.f = self.fun(self.t, self.y)

  File "...\anaconda3\lib\site-packages\scipy\integrate\_ivp\base.py", line 138, in fun
    return self.fun_single(t, y)

  File "...\anaconda3\lib\site-packages\scipy\integrate\_ivp\base.py", line 20, in fun_wrapped
    return np.asarray(fun(t, y), dtype=dtype)

  File "...\anaconda3\lib\site-packages\pde\solvers\scipy.py", line 74, in rhs_helper
    return rhs(state_flat.reshape(shape), t).flat  # type: ignore

  File "...\anaconda3\lib\site-packages\numba\core\dispatcher.py", line 420, in _compile_for_args
    error_rewrite(e, 'typing')

  File "...\anaconda3\lib\site-packages\numba\core\dispatcher.py", line 361, in error_rewrite
    raise e.with_traceback(None)

TypingError: Cannot capture the non-constant value associated with variable 'y' in a function that will escape.

Answer 1

由于还没有人发布答案，我设法通过使用 scipy odeint 和线的方法来求解 PDE，即通过离散化拉普拉斯算子，得到了一个最小的工作示例，然后将微分方程包裹在 fsolve 中得到 dydt:

import numpy as np
from scipy.integrate import odeint
from scipy.optimize import fsolve

a=1
b=1
t = np.linspace(0,1,100)
y0 = np.asarray([1,0.5,0]) #starting values
N=len(y0)

def laplace(y):
    """ Laplace operator in cartesian coordinate system """
    d2y_0 = [2*(y[1] - y[0])] #at i=0
    d2y_N = [2*(y[-2] - y[-1])] #at i=N
    d2y_i = [y[i+2] - 2*y[i+1] + y[i] for i in range(N-2)] #elsewhere
    return np.asarray(d2y_0 + d2y_i + d2y_N)

def rhs(y, dydt, t):
    """ RHS of PDE including dydt term """
    return a*laplace(y) + b*laplace(dydt)

def pde(y, t):
    """ Function that solves for dydt using root finding """
    return fsolve(lambda dydt : dydt - rhs(y, dydt, t),np.zeros(N))

#calculate
sol = odeint(pde, y0, t)

#plot odeint with fsolve
import matplotlib.pyplot as plt
plt.figure()
plt.plot(t,sol)
plt.grid(ls='--')
plt.xlabel('$t$')
plt.ylabel('$y_i$')
plt.title('odeint with fsolve')
plt.legend(["$i=$"+str(i) for i in range(N)])

#plot theoretical
u10=y0[0]
u20=y0[1]
u30=y0[2]

E1=np.exp(-2*a*t/(1+2*b)) #1st exponential term
E2=np.exp(-4*a*t/(1+4*b)) #2nd exponential term
u1 = 1/4*((1+2*E1+E2)*u10 - 2*(E2-1)*u20 + (1-2*E1+E2)*u30)
u2 = 1/4*((1-E2)*u10      + 2*(E2+1)*u20 + (1-E2)*u30)
u3 = 1/4*((1-2*E1+E2)*u10 - 2*(E2-1)*u20 + (1+2*E1+E2)*u30)

plt.figure()
plt.plot(t,u1,t,u2,t,u3)
plt.grid(ls='--')
plt.xlabel('$t$')
plt.ylabel('$y_i$')
plt.title('Theoretical')
plt.legend(["$i=$"+str(i) for i in range(N)])

请注意，同样的拉普拉斯离散化方法允许我们求解一个 ODE 系统，以获得用于验证我们的数值计算的精确解析解（对于 N=3 大小的系统）：

>>> np.allclose(sol[:,0],u1)
True
>>> np.allclose(sol[:,1],u2)
True
>>> np.allclose(sol[:,2],u3)
True

似乎它按预期工作。