Solving a parabolic PDE numerically with the spectral method

Question

As a homework question, I am asked to numerically solve the PDE on $x \in [-3, 3]$ and $t \in [0, 0.5]$ with the spectral method.

$$ \frac{\partial p}{\partial t} = (12x^2-4) p + \left[4x(x^2-1)+0.1\right] \frac{\partial p}{\partial x} + \frac{1}{\beta} \frac{\partial^2 p}{\partial x^2} \\ $$

with the initial condition $p(x,0)=\frac{1}{\sqrt{2\pi}} \exp(-\frac{x^2}{2})$ and zero boundary condition $p(-3,t)=p(3,t)=0$.

Attempt: By the Chebyshev expansion of the first kind $p(x, t) = \sum_{n=0}^N \tau_n(t) T_n(x)$,

\begin{equation} \begin{aligned} \frac{\partial p}{\partial t} =& (12x^2-4) p + \left[4x(x^2-1)+0.1\right] \frac{\partial p}{\partial x} + \frac{1}{\beta} \frac{\partial^2 p}{\partial x^2} \\ \sum_{n=0}^N \tau_n'(t) T_n(x) =& (12x^2-4) \sum_{n=0}^N \tau_n(t) T_n(x) + \left[4x(x^2-1)+0.1\right] \sum_{n=0}^N \tau_n(t) T_n'(x) - \frac{1}{\beta} \sum_{n=0}^N \tau_n(t) T_n''(x) \\ =& \left(2T_0(x) + 6T_2(x)\right) \sum_{n=0}^N \tau_n(t) T_n(x) + \left(0.1T_0(x) - T_1(x) + T_3(x)\right) \sum_{n=0}^N \tau_n(t) T_n'(x) - \frac{1}{\beta} \sum_{n=0}^N \tau_n(t) T_n''(x) \\ \end{aligned} \end{equation}

Ideally, this would yield one ODE for each $\tau_n(t)$ by the uniqueness of Chebyshev expansion, but the polynomials $12x^2-4$ and $4x(x^2-1)+0.1$ mess everything up.

Since I need a numerical result, I don't think I am supposed to perform the multiplication of Chebyshev polynomials by hand. What should I do next?

Answer 1

So this is the discretisation of time that I meant to suggest in the comments: \begin{equation*} \frac{p(x,t)-p(x,t-\Delta t)}{\Delta t} = \left(12 x^{2} - 4\right) p(x,t) +\left(4 x \left(x^{2}- 1\right)+ 0.1\right)\frac{\partial p(x,t)}{\partial x} +\frac{1}{\beta} \frac{\partial^{2} p(x,t)}{\partial x^{2}}. \end{equation*} That becomes \begin{equation*} Ap=b \end{equation*} with \begin{equation*} A = 1-\left(12 x^{2} - 4\right)\Delta t -\left(4 x \left(x^{2}- 1\right)+ 0.1\right)\Delta t\frac{\partial}{\partial x} -\frac{1}{\beta} \Delta t\frac{\partial^{2}}{\partial x^{2}} \end{equation*} \begin{equation*} b(x) = p(x,t-\Delta t). \end{equation*} In Python the code might look something like

import numpy as np

T = 0.5
n = 10  # The number of time steps.
dt = T / n

N = 10  # The number of collocation points in the spatial dimension.
beta = whatever

# You need to work out the vector of collocation points x and differentiation matrix D.
# Look for the routine 'lobatto' in section 4.6.4 of Numerical Recipes and the routine
# 'weights' in section 20.7.8.

A = np.diag(1 - (12 * (x ** 2) - 4) * dt) -\
    (4 * x * (x ** 2 - 1) + 0.1) * dt * D - (1 / beta) * dt * D @ D
A = A[1:N-1, 1:N-1]  # That the boundary conditions are zero makes things simple.

def step(p):
    b = p[1:N-1]
    return np.hstack([0, np.linalg.solve(A, b), 0])

p = np.exp(-(x ** 2) / 2) / np.sqrt(2 * np.pi)

for i in range(n):
    p = step(p)

Edit: it seems the Chebyshev-Gauss-Lobatto points can be worked out more easily than by transcribing that routine from Numerical Recipes. This page says :

The Chebyshev–Gauss–Lobatto points are the zeros of $\left(1-x^{2}\right)T_{n}'(x)$. Using the property $T_{n}(x)=T_{n}(\cos\theta)=\cos(n\theta)$, these can be shown to be at $x_{j}=\cos(\pi j/n)$.

Edit 2: I was asked in the comments to make a plot ... here it is

produced by

T = 0.5
n = 20  # The number of time steps.
dt = T / n

N = 50  # The number of collocation points in the spatial dimension.
beta = 1

x = 3 * np.cos(np.linspace(-np.pi, 0, N))
w = [weights(xi, x, 2) for xi in x]
D = np.array([c[:, 1] for c in w]) / 3

A = np.eye(N) - D @ (np.diag(4 * x * (x ** 2 - 1) + 0.1)
    + (1 / beta) * D) * dt
A = A[1:N-1, 1:N-1]  # That the boundary conditions are zero makes things simple.

def step(p):
    b = p[1:N-1]
    return np.hstack([0, np.linalg.solve(A, b), 0])

p = np.exp(-(x ** 2) / 2) / np.sqrt(2 * np.pi)

for i in range(n):
    p = step(p)

grid = np.linspace(-3, 3, 100)
w = np.array([weights(gi, x, 0) for gi in grid]).reshape(100, N)

fig, ax = plt.subplots()
ax.plot(grid, w @ p)

Answer 2

I oriented myself to the approach described by this post entitled "Solving 2D+1 PDE with Pseudospectral in one direction with periodic boundary condition?" and used:

With[{p = p[x, t], mid = mid[x, t]},
 eq = D[p, t] == D[mid, x];
 eqadd = mid == (4 x (x^2 - 1) + 0.1) p + D[p/b, x];]
ic = p[x, 0] == Exp[-x^2/2]/Sqrt[2 Pi];
icadd = mid[x, 0] == eqadd[[-1]] /. p -> Function[{x, t}, Evaluate@ic[[-1]]];
bc = {p[-3, t] == 0, p[3, t] == 0};

mol[n:_Integer|{_Integer..}, o_:"Pseudospectral"] := {"MethodOfLines", 
  "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n, 
    "MinPoints" -> n, "DifferenceOrder" -> o}}

tst = NDSolveValue[{eq, eqadd, ic, icadd, bc}, p, {x, -3, 3}, {t, 0, 0.5}, 
    Method -> mol[300, 2]]; // AbsoluteTiming

Plot3D[tst[i, j], {i, -3, 3}, {j, 0, 0.5}, PlotRange -> All, PlotPoints -> 200]

I owe the key hints (1, 2) to make it running to the Mathematica community. As a result we obtain the following plots when choosing $b=1$ (left plot), $b=200$ (right plot):

We obtain the following two plots using Plot[tst[i, 0.5], {i, -3, 3}, PlotRange -> All] again for $b=1$ (left plot), $b=200$ (right plot):

Update (2022-01-10):

Fix according to 2 including usage of "Pseudospectral" option.