6
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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}$. $\endgroup$ – Ben May 28 '13 at 23:42
  • $\begingroup$ Thanks. But step by step approach in unstable for any grid i take. $\endgroup$ – Uroš May 29 '13 at 5:51
2
$\begingroup$

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} $$

Be warned though, this matrix has constant null space that you need to remove before solving.

What sort of problems are you running into? Maybe if you post them here, we could be of more help.

$\endgroup$
  • $\begingroup$ we tried doing this but could not produce results in Matlab. $\endgroup$ – Uroš May 29 '13 at 5:52

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.