Solving a PDE with 2 variables, with one variable whose derivative with respect to space is only known

Question

I am trying to solve a PDE of the form

$$\frac{\partial u}{\partial t} = D \frac{\partial^2u}{\partial x^2} + k\ \ \ (1)$$

where only $k$'s first derivative with respect to $x$ is known

$$\frac{\partial k}{\partial x}=-5u\ \ \ \ (2)$$

The boundary conditions are:

$$\frac{\partial u}{\partial x}\left(t,x=0\right)=0\ \ \ \ (3)$$

$$\frac{\partial u}{\partial x}\left(t,x=\left(J-1\right)\Delta x\right)=0\ \ \ \ (4)$$

Using the method of lines, the following discretisation for $u$ is obtained:

$$\frac{\partial u}{\partial t}\left(t,x=0\right)=2D\frac{u^1-u^0}{\Delta x^2}+k^0\ \ \ \ (5)$$

$$\frac{\partial u}{\partial t}\left(t,x=j\Delta x\right)=D\frac{u^{j+1}-2u^j+u^{j-1}}{\Delta x^2}+k^j\ \ \ \ (6)$$

$$\frac{\partial u}{\partial t}\left(t,x=\left(J-1\right)\Delta x\right)=2D\frac{u^{J-2}-u^{J-1}}{\Delta x^2}+k^{J-1}\ \ \ \ (7)$$

Using these 3 equations and scipy.integrate.solve_ivp allows me to solve the diffusion part of equation (1). However, to get k, I need to first calculate $\frac{\partial k}{\partial t}$ while I only know $\frac{\partial k}{\partial x}$. How can I do that? I've written the following Python code below to solve $\frac{\partial u}{\partial t}$ but I'm missing the code to calculate $\frac{\partial k}{\partial t}$.

import scipy.integrate as scint
import numpy as np

x = np.linspace(0, 1, 11)
dx = np.diff(x)[0]
D = 10
u0 = np.exp(-5 * x)
k0 = np.zeros(len(x))
t = np.linspace(0, 1, 101)


def get_dydt(_t, y):
    """ Return the derivative dy/dt at different points in space """

    u, k = np.split(y, 2)

    # Calculate du/dt
    dudt = np.empty_like(u)
    dudt[0] = 2 * (u[1] - u[0])
    dudt[1:-1] = np.diff(u, 2)
    dudt[-1] = 2 * (u[-2] - u[-1])
    dudt = dudt * D / dx ** 2 + k

    # Calculate dk/dt
    dkdt = None  # how can I calculate this?

    return np.concatenate([dudt, dkdt])


solver = scint.solve_ivp(get_dydt, [t[0], t[-1]], np.concatenate([u0, k0]), method='Radau')

Answer 1

The equations you have are what is called a "differential-algebraic equation" (DAE) because you have only time derivatives for one of the variables (namely, $u$) but not for the other (namely, $k$). The prototypical case of this kind of equation is the time dependent Stokes equations, which has time derivatives for the velocity but not the pressure. The way one generally interprets this is that equation (2) has to hold at all instants $t$.

There are whole classes of methods for DAEs, but in your case I would try to understand what people do for time discretizations of the Stokes equations. These techniques will be easily applicable to your case as well.

Answer 2

As Maxim Umansky mentioned in his comment if you find the $k$ from your second equation, you would end up with this differential-integral equation:

$$\frac{\partial u}{\partial t} = D \frac{\partial^{2} u}{\partial x^{2}} -5 \int_{0}^{x} u(x^{'},t) d x^{'}$$

The discretized form of this equation is:

$$\frac{u_{x}^{t+\Delta t} - u_{x}^{t}}{\Delta t} = D \frac{u_{x+\Delta x}^{t} + u_{x-\Delta x}^{t} - 2u_{x}^{t}}{\Delta x} - 5 \sum_{j=0}^{\frac{x}{\Delta x}} a_{j} u_{j \Delta x}^{t}$$

Note that $a_{j}$ are the coefficients in the integration operator that would be determined when you choose which numerical integration you want to use here such as trapezoidal rule, etc.

The discretized equations could be rewritten in vector format as:

$$| u \rangle^{t+\Delta t} = (I + \Delta t\mathcal{L}) |u \rangle^{t}$$

Where $| u \rangle^{t}$ and $| u \rangle^{t+\Delta t}$ are vectorized form of the nodal values of $u$ in your entire computational domain (e.g. a 1D line here) at timesteps $t$ and $t+\Delta t$ respectively. Also $\mathcal{L}$ is a square matrix $(N+1) \times (N+1)$ that is extracted by combining the 1D discretized Laplace operator and your integration operator.

The above discretized equation is only correct for nodes inside your computational domain ($\Omega$). For nodes on your boundary ($\partial \Omega$), you might want to add your boundary conditions. As a result, the final linear equation would be:

$$A | u \rangle^{t+\Delta t} = b$$

Where $A$ is assembled as: $$A = I + \partial \mathcal{L}$$ Where $\partial \mathcal{L}$ is a matrix that only has non-zero values on nodal points of your boundary and $b$ is assembled as: $$b = (I + \Delta t \mathcal{L}) | u \rangle^{t} + \partial b$$ where $\partial b$ is a vector that only has non-zero values on nodal points of your boundary.

As I mentioned in my comment the problem with this approach is that $A$ is not a sparse matrix because of having a integration operator in your PDE. This might not be a serious concern here if your $N$ the number of nodes would not be a large value. But, solving a dense linear equation would be really expensive computationally if $N$ becomes a large value.

Answer 3

Assuming that you mistake a $t$ for an $x$ in the right hand side of the first equation, you have the following system of equations.

\begin{align} \frac{\partial u}{\partial t} = D \frac{\partial^2u}{\partial x^2} + k\\ 0 = \frac{\partial k}{\partial x} - 5u\, . \end{align}

To get a unique solution you would need 2 boundary conditions for $u$ and 1 boundary condition for $k$. You would also need an initial condition for $u$.

Having the appropriate conditions you could then step in time with something like:

1. Compute $k^{0}_i$ for $t = 0$ from initial condition for $u$ and boundary condition.
2. Compute $u_i^t$ from $k_i^{t-1}$, $u_i^{t-1}$ and boundary conditions.
3. Compute $k_i^t$ from $u_i^t$ and boundary condition.