Is there a way to numerically solve the following PDE in Python?
The second term on the RHS has a derivative with respect to time as well as space.
I tried using Py-PDE package in Python, it solves only the form dy(x,t)/dt = f(y(x,t)) so I tried to use a root finding algorithm similar to scipy fzero to get the solution to dy(x,t)/dt - f(y(x,t),dy(x,t)/dt) = 0 (solving for 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
However, when the program tries to solve the PDE it throws an error:
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.
Since no one has posted an answer yet, I managed to get a minimal working example by using scipy odeint with a method of lines to solve the PDE, that is, by discretizing the Laplace operator, and then wrapping the differential equation inside fsolve to get 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)])
Note that the same Laplace discretization method allows us to solve a system of ODEs to get the exact analytical solution with which we verify our numerical calculation (for a N=3 size system):
>>> np.allclose(sol[:,0],u1)
True
>>> np.allclose(sol[:,1],u2)
True
>>> np.allclose(sol[:,2],u3)
True
Seems it works as intended.