# Solvers for odd order PDE finite difference discretisation

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.

• 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? Feb 4, 2022 at 13:32
• 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". Feb 4, 2022 at 13:34

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.

• 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. Feb 4, 2022 at 1:40
• 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. Feb 4, 2022 at 1:46
• 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. Feb 4, 2022 at 17:53
• 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. 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.