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')
  • 2
    $\begingroup$ What is the question here? Are you asking for help to debug your code or check for correctness? $\endgroup$ Apr 28 at 15:36
  • 3
    $\begingroup$ I suspect you meant second $x$ derivatives in the right hand side of the first equation? $\endgroup$ Apr 28 at 15:48
  • $\begingroup$ @WolfgangBangerth Thanks, I've corrected it. $\endgroup$ Apr 28 at 17:32
  • $\begingroup$ @AbdullahAliSivas I've rephrased my question to make it a bit more clear. Basically, i do not know how to calculate dk/dt such that I can numerically solve the Equation (1). $\endgroup$ Apr 28 at 17:37

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.


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.

  • 1
    $\begingroup$ That looks fine if the integral sum in the discretized form of the equation runs not to the end of the grid but to a specific j corresponding to the given point x. Also, for explicit time step it should not be a problem if the matrix is not sparse, you don't need to invert it but only multiply it by a state vector. $\endgroup$ Apr 29 at 23:36
  • $\begingroup$ @MaximUmansky Well I think $A$ would be a pseudo lower triangular matrix and at least it’s sparsity is low. $\endgroup$ Apr 30 at 3:26

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.
  • $\begingroup$ The RHS, combining $𝑒_{π‘₯π‘₯}$ and βˆ«π‘’π‘‘π‘₯, is just a linear operator acting on u, L(u), so it can be represented by a matrix multiplied by 𝑒⃗ , where 𝑒⃗ is a set of values of u on the grid points. So performing a time step is just solving a linear problem, $𝑒_{𝑛𝑒𝑀} = 𝑒_{π‘œπ‘™π‘‘}+𝑑𝑑𝐿(𝑒_{π‘œπ‘™π‘‘})$, or implicitly, $𝑒_{𝑛𝑒𝑀}=𝑒_{π‘œπ‘™π‘‘} + 𝑑𝑑𝐿(𝑒_{𝑛𝑒𝑀})$. $\endgroup$ Apr 29 at 6:19
  • $\begingroup$ @MaximUmansky, yes. That can be done that way as well. $\endgroup$
    – nicoguaro
    Apr 29 at 13:42
  • $\begingroup$ @MaximUmansky Well, I think a possible problem with your approach is that the integral operator acting on $u$ needs to have access to all nodal values of $u$ across the computational domain. As a result, the $L$ operator would not be a sparse matrix and would be a dense matrix. For 1D problems such as this, this might not be a problem at all but I don't think it would be a good idea to use it for 2D or 3D problems. $\endgroup$ Apr 29 at 18:02

Your Answer

By clicking β€œPost Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.