I am used to solving elliptic PDEs of even order. I was wondering what would one do for odd order PDEs. Notably the discretisations of those results in unsymmetric matrices. I tried solving the following simple problem:

$$\frac{d^3u}{dx^3}(x) = 0,\, x\in \Omega \quad u(x) = f(x), \, x \in \Gamma_D, \, \partial^{l}_nu(x) = 0, x\in\Gamma_N$$

I would generally expect to get some quadratic spline as the solution to the above. But I got a very oscillatory solution which I believe is due to numerical instability. I used the discretisation: $$\frac{d^3u}{dx^3} = \begin{bmatrix} -\frac{1}{2} & 1 & 0 & -1 & \frac{1}{2} \end{bmatrix} *u+O(h^2)$$

Since the system matrix $L$ is unsymmetric, I use a CGNR solver to solve $L^TLx = L^Tb$ instead, but nevertheless things don't seem to work out. I have never dealt with odd-order PDEs numerically. Are there some references discussing this? Am I supposed to use a specific discretisation like an upwind scheme? I don't think this makes any sense as the factor in front of $\frac{d^3u}{dx^3}$ seems to be irrelevant considering that the right-hand side is zero. Would the above problem be stabilized by solving:

$$\left |\frac{d^3u}{dx^3}\right | = 0,$$

instead? Clearly I cannot write the above as $y=Lx$ anymore and would have to potentially extend it to $\partial_t u = \left|\frac{d^3u}{dx^3}\right|$, e.g. apply explicit Euler and send the time to infinity. I just want to make sure whether there are known methods to handle odd order PDEs numerically before I try some stupid ideas.

  • $\begingroup$ Since you write "solve" everywhere and not "numerically solve", I'm wondering if you realize that the problem you have written down can be solved exactly in a straightforward way? And if you are intent on solving it numerically, I think a key question that you have not answered is, how did you impose the boundary conditions? $\endgroup$ Feb 4, 2022 at 13:32
  • $\begingroup$ Oh, I see that you're trying to instead solve $u_t + u_{xxx} = 0$ by integrating in time with Euler's method. That is indeed unstable. There are many ways to integrate this stably. Since it is a BVP, there is no need to discretize in "time". $\endgroup$ Feb 4, 2022 at 13:34

2 Answers 2


The matrix for $d/dx$ is skew-symmetric for a suitable discretization scheme, and the same is true for its odd powers.

Moreover, the operator $i d/dx$ is hermitian, and again, by using a skew-symmetric discretization, all its powers are hermitian as well.

If you then have a PDE that only contains odd derivatives of $x$, you can consider the solution function $v(x)=i u(x)$, and bring the imaginary unit to the operator side and thus obtain a hermitian problem, which can be solved , e.g., by conjugate gradient again.

That's it for now, as I'm on my mobile phone and typing is exhausting. Just leave a comment and ask, I'll write more once I'm back at my laptop.

  • $\begingroup$ Thank you for the answer. It would be great if this can be solved in an uncomplicated manner. I ran several tests with explicit Euler in the meantime but they always diverged. The thing that worked was scaling back to $|\frac{du}{dx}| =\partial_t u$ and applying the Rouy-Tourin discretisation with explicit Euler. This resulted in a piecewise constant polynomial as expected. The main issue is that $u_{xxx}$ seems to have the same or worse discretisation problems, and I am not aware of a Rouy-Tourin scheme for 3rd derivatives. $u_{xxx}=0$ is preferable over $|u_{xxx}|=u_t, \, T=\infty$ too. $\endgroup$
    – lightxbulb
    Feb 4, 2022 at 1:40
  • $\begingroup$ I will try the hermitian CG approach in the following days and see what I get. I have also considered using a spectral approach, but it looks more painful handling arbitrary Dirichlet constraints with it and would be slower with the multiple FFTs required than with CG with a spatial discretisation, especially for large problem sizes. $\endgroup$
    – lightxbulb
    Feb 4, 2022 at 1:46
  • $\begingroup$ I did some experiments, most notably I formed the matrix even for the discretisation of $u_x$ with its Dirichlet conditions as a toy example, and then computed the eigenvalues of $A^TA$. Turns out that the matrix with Dirichlet conditions is singular, which explains why CGNR was not converging. I still have to figure out what this means for the problem. Whether the conditions are insufficiently many, whether they are unachievable or something else. $\endgroup$
    – lightxbulb
    Feb 4, 2022 at 17:53
  • $\begingroup$ I did some more experiments and it is non-singular for upwind schemes however, as one would expect. I need to look further into how this generalizes to $u_{xxx}$ though. $\endgroup$
    – lightxbulb
    Feb 4, 2022 at 18:06

As mentioned in my comments to davidhigh's answer, using the discretisation in my question results in a skew-symmetric matrix which when modified to add Dirichlet conditions becomes singular. To better understand the problem I looked at discretisations of:

$$\frac{du}{dx} = 0$$

and those require an upwind scheme. Applying the same idea to $\frac{d^3u}{dx^3} = 0$ problem I derived the scheme: $\begin{bmatrix} 1 & -3 & 3 & -1\end{bmatrix}$. Using said discretisation and adding Dirichlet constraints does not result in a singular system anymore. Nevertheless explicit schemes (e.g. Richardson, Jacobi, etc.) do not seem to converge for the resulting unsymmetric problem:

$$Ax = b$$

Instead I decided to symmetrize the problem by:

$$AA^Ty = b, \, x = A^Ty.$$

If CG is applied to the above form with some modification one can arrive at the CGNE algorithm. Using CGNE resulted in a convergent iteration.


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.