# Boundary value method for equation $u_{tt} = u_{xxx}$

I have this funny equation

$$\frac{\partial^2 u}{\partial t^2} = \frac{\partial^3 u}{\partial x^3}, \qquad x \in [0,1], \qquad t \in (0,T]$$ with initial conditions $u(x,0) = \sin(2\pi x)$, $\frac{\partial u}{\partial t}(x,0) = 0$ and periodic boundary conditions: $u(0,t) = u(1,t)$, $\frac{\partial u}{\partial x}(0,t) = \frac{\partial u}{\partial x}(1,t)$ and $\frac{\partial^2 u}{\partial x^2}(0,t) = \frac{\partial^2 u}{\partial x^2}(1,t)$. The exact solution of this model is $$u(x,t) = \frac{1}{2}[e^{2 \pi^{\frac32}t} \sin (2\pi x - 2\pi^{\frac32} t) + e^{-2\pi^{\frac32}t} \sin (2\pi x + 2 \pi^{\frac32}t)]$$ I would like to use boundary value technique to numerically solve this equation, i.e. i would like to write the following system of $2MN$ unknowns in matrix form and solve it in Matlab.

I partition time as $0 = t_0 < t_1 < \dots < t_{N} = 0.1$. Step size is of length $h=\Delta t =\frac{0.1}{N}$, i.e. $t_n = t_0 + nh$, $n=1,2,\dots,N$.

I partition space as $0 = x_0 < x_1 < \dots < x_{M} = 1$. Step size is of length $h'=\Delta x =\frac{1}{M}$, i.e. $x_i = x_0 + ih'$, $i=1,2,\dots,M$.

I rewrite the equation as system,

\begin{align*} u_t &= v, \\ v_t &= u_{xxx} \end{align*}

and using central approximation for time derivative and second order scheme for space derivative

\begin{align*} \frac{u_i^{n+1} - u_i^{n-1}}{2\Delta t} &= v_i^n \\ \frac{v_i^{n+1} - v_i^{n-1}}{2\Delta t} &= \frac{u_{i+2}^n - 2 u_{i+1}^n +2u_{i-1}^n -u_{i-2}^n}{2\Delta x^3} \end{align*}

To add an interesting fact: this can't be solved via step by step approach. At least as far as i know. BTCS and FTCS methods do not work. Also other schemes are not stable (eigenvalues of step by step matrix form a cross in complex plane and are outside stability regions of the schemes for any $\Delta t$ and $\Delta x$). This is why i would like to use the approach of solving system of equations to get values at grid points.

Here is my try deriving system of equation when i introduce vector $v$. I define $u^{n}:=[u_1^{n},u_2^{n},\dots,u_M^{n}]^T$ and $v^{n}:=[v_1^{n},v_2^{n},\dots,v_M^{n}]^T$ and then $w^n= [(u^n)^T,(v^n)^T]^T$ and then discretizing the space yields

$$\begin{pmatrix} u^n \\ v^n \end{pmatrix}_t = \begin{pmatrix} 0 & I \\ A & 0 \end{pmatrix}w^n =: D w^n$$ where $$A=\frac{1}{2 \Delta x^3} \begin{pmatrix} 0 & -2 & 1 & 0 & \cdots & 0 & -1 & 2 \\ 2 & 0 & -2 & 1 & 0 & \cdots & 0 & -1\\ -1 & 2 & 0 & -2 & 1 & 0 & \cdots & 0 \\ 0 & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & 0 \\ 0 & \cdots & 0 & -1 & 2 & 0 & -2 & 1 \\ 1 & 0 & \cdots & 0 & -1 & 2 & 0 & -2 \\ -2 & 1 & 0 & \cdots & 0 & -1 & 2 & 0 \end{pmatrix}$$ using periodic boundary conditions. Central discretization in time yields

$$w^{n+1} - w^{n-1} - 2\Delta tD w^n = 0$$

In the last step i use backward in time derivative approximation to get

$$-w^{N-1} - (\Delta tD - I) w^{N} = 0$$

Lastly, defining $w=[(w^0)^T,(w^1)T,\dots,(w^N)^T]^T$,

\begin{align*} C &=\begin{pmatrix} I & 0 & 0 & 0 & \dots & 0 & 0 & 0 \\ -I & -2\Delta tD & I & 0 & \dots & 0 & 0 & 0 \\ 0 & -I & -2\Delta tD & I & \ddots & 0 & 0 & 0 \\ 0 & 0 & -I & -2\Delta tD & \ddots & \ddots & 0 & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots & \ddots & \ddots & \vdots \\ 0 & 0 & 0 & \ddots & \ddots & -2\Delta tD & I & 0 \\ 0 & 0 & 0 & 0 & \ddots & -I & -2\Delta tD & I \\ 0 & 0 & 0 & 0 & \dots & 0 & -I & I-\Delta tD \\ \end{pmatrix}\\ \end{align*} and $c = ((w^0)^T,0...,0)^T$ i have to solve

$$C w = c$$

Edit: It works in Matlab now and the results are nice. Just matrix inversion is not fast.

• What stops you from using a standard time-stepping scheme, with a centered second order stencil in time and the second order stencil written above in space? The approach should be identical to that for the wave equation (e.g., see here), except with a finite difference approximation of $u_{xxx}$ instead of $u_{xx}$.
– Ben
May 28 '13 at 23:42
• Thanks. But step by step approach in unstable for any grid i take.
– Uroš
May 29 '13 at 5:51

The periodic boundary condition is really not a boundary condition; it's merely a statement that the domain is infinite. What that means is that you need to discretize your PDE is space at each and every point, including $i=0$ and $i=M$ and simply use periodic indecies, i.e. $i>M: i\leftarrow i-M$ and $i<0: i\leftarrow i+M$.
As for the time discretization, how about the following? $$\frac{u_i^{n+1}-2u_i^n+u^{n-1}_i}{\Delta t^2} = \frac{u_{i+2}^{n+1}-2u_{i+1}^{n+1} + 2u_{i-1}^{n+1} - u_{i-2}^{n+1}}{2\Delta x^3}$$