Suppose we have a fourth-order parabolic PDE

$$\partial_t u + a(t)\Delta( \Delta u) - div( b(x,t)\nabla u) =0$$

We split the equation into two second-order equations by introducing $w = \Delta u$ thus

$$\partial_t u + a(t) \Delta w - div(b(x,t)\nabla u) =0$$ $$-\Delta u + w =0$$

I discretise the first equation as $$ \dfrac{u^{n+1} - u^n}{\Delta t} + a(t^n) \Delta w^{n+1} - div(b(x,t^n) \nabla u^{n+1} =0 \ \ \ \ \ \ \text{(0 )}$$ My issue is with the second equation. Should it be \begin{equation} -\Delta u^{n+1} + w^{n+1} =0 \ \ \ \ \ \ \text{(1)} \end{equation} or this one \begin{equation} -\Delta u^{n} + w^{n+1} =0 \ \ \ \ \ \ \text{(2)} \end{equation} I was thinking of an algorithm for the numerical simulations and I think equation (2) makes more sense. Since we first solve it to obtain $w^{n+1}$ then we use $w^{n+1}$ in equation (0) to find $u^{n+1}$. If we use equation (1) I do not see how we would determine $w^{n+1}$ since $u^{n+1}$ is also not known.

Any help would be appreciated.

  • $\begingroup$ There are no initial and boundary conditions specified, so the current formulation is ill-posed. Also you can directly discretise the 4th order formulation if you wish, i.e. you may avoid the substitution if you want. $\endgroup$
    – lightxbulb
    Aug 22, 2023 at 20:24

1 Answer 1


The introduction of $w$ is just a reformulation of your initial problem.

If you use (2), this means that you make the 4th-order diffusion term explicit, which may potentially lead to stability issues. Your scheme (0) would then correspond to an implicit-explicit scheme using explicit or implicit Euler methods for the different terms.

If you use (1) instead, then your overall temporal discretisation corresponds to the implicit Euler method (aka backward Euler), which is usually a much more stable scheme.

Alternatively, you can consider having the following viewpoint on your problem. First, you discretise in space to obtain a set of coupled ODEs $u'=f(u)$ where f contains the spatial discretisation of the RHS (4-th order diffusion term and divergence). Then, you can apply any temporal discretisation scheme you want directly on $f$, e.g. implicit Euler: $$\dfrac{u^{n+1}- u_{n}}{\Delta t} = f(u^{n+1})$$

Or you can split the RHS: $u' = f_1(u) + f_2(u)$, with $f_1$ being the diffusion term, and $f_2$ the divergence. Then you can apply implicit-explicit (IMEX) schemes, or other approaches, e.g. splitting. Note that what you did is simply inserting $w(u)=\Delta u$ in $f_1$, i.e. you obtain a formulation of the form $u' = f_1^*(w(u)) + f_2(u)$.

  • $\begingroup$ Thank you! This was helpful. However, if we choose equation (1), how can we solve for $w^{n+1}$ since we have no previous step of it, while $u^{n+1}$ is also still unknown? $\endgroup$
    – Thede
    Aug 22, 2023 at 11:50
  • $\begingroup$ This depends on the solution method you use to solve (0). I guess your issue is more on the technical implementation of the scheme. Usually, the overall discretised ODEs (0) is solved via a Newton method, which iterates on $u^{n+1}$, each iteration requiring (usually) one call of $f$, until (0) is satisfied. Then, at each such call of $f$, the solver gives you its current guess of $u^{n+1}$ from which you can compute $w^{n+1}$. Here $w$ is purely an intermediate variable, e.g. some internal "cuisine" for your scheme, i.e. it should only exist within $f$. $\endgroup$
    – Laurent90
    Aug 22, 2023 at 16:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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