I'm doing a project where I need to code a 4th-order accurate finite difference method for solving a simple wave equation: $u_{tt}-\nabla \cdot (a(x,y)\nabla u) = f(x,y,t)$ in a square domain. I know the forcing function $f$, my initial conditions $u(x,y,0)$ and $u_t (x,y,0)$, and $a(x,y)$ at each grid point in my domain. Any advice for stable 4-th order method?


  • $\begingroup$ May I ask if you want 4th order in time or in space or in both? $\endgroup$ – hypernova May 16 '18 at 0:41
  • $\begingroup$ Both. Basically, I set the spatial grid size and then calculate the required time step (CFL condition). Can I used Method of Lines to get fourth-order? So far, I've been using a 2nd-order accurate method O($h^2$), where $h$ is the grid size (same for both $x$ and $y$ dimensions. $\endgroup$ – D.B. May 16 '18 at 0:50
  • $\begingroup$ Sure. Thank you for your clarification. $\endgroup$ – hypernova May 16 '18 at 1:02
  • $\begingroup$ The main problem that I'm seeing is that for a fourth-order multi-step finite difference approximation for $u_{tt}$, I need to know the previous two values of the solution $u^{i-2}$ and $u^{i-1}$. How might I do this if I only know $u_{t}$ and $u_{tt}$ at $t=0$. $\endgroup$ – D.B. May 16 '18 at 1:58
  • $\begingroup$ Apologies for being late. I'm at a meeting. I'll get back to you in two hours. $\endgroup$ – hypernova May 16 '18 at 2:44

Spatial discretization for higher order

The explanation in this section takes the 1-D case for instance.

Suppose you have $$ u''(x)=f(x), $$ where $f(x)$ could be any given function of $x$, or any $u$-related values (e.g., $f(x)=u_t(x)$, where $t$ is regarded as an auxiliary parameter).

A usual second-order scheme reads $$ \frac{u_{j+1}-2u_j+u_{j-1}}{h^2}=f_j. $$ By contrast, a fourth-order scheme could be realized by taking $$ \frac{u_{j+1}-2u_j+u_{j-1}}{h^2}=\alpha f_{j+1}+\beta f_j+\gamma f_{j-1}, $$ where $\alpha$, $\beta$ and $\gamma$ could be determined by requiring the scheme to be a fourth-order scheme to approximate $u''=f$. You may figure out that $$ \frac{u_{j+1}-2u_j+u_{j-1}}{h^2}=\frac{f_{j+1}+10f_j+f_{j-1}}{12}, $$ which is known as the fourth-order compact scheme.

This method also applies to the $\left(au'\right)'$ case and the 2-D cases.

In addition, if the boundary condition prescribed to $u$ is the null Dirichlet or the null von Neumann (i.e., $u=0$ or $u'=0$ on the boundary), a more efficient choice would be the spectral method, which could provide you with spectral accuracy.

Time discretization for higher order

For the sake of stability, the second-order Crank-Nicholson usually suffices for most cases. Multistep methods work well for ODEs, but could perform badly for PDEs. In fact, due to the appearance of $\Delta$, a multistep scheme could be unconditionally unstable.

By contrast, the integrating factor method could be helpful. This method requires the spatial discretization at first. Say you have $u_t=\Delta u+f$. A spatial discretization gives you $$ \frac{\rm d}{{\rm d}t}\mathbf{u}=A\mathbf{u}+\mathbf{f}, $$ where $\mathbf{u}$ is the spatial discretization of $u$, $A$ is a constant matrix, while $\mathbf{f}$ is the spatial discretization of $f$ composed with boundary terms of $u$. Now this equation becomes an ODE, and is equivalent to $$ \frac{\rm d}{{\rm d}t}\left(e^{-tA}\mathbf{u}\right)=e^{-tA}\mathbf{f}. $$ As such, the numerical errors with respect to $t$ all come from the quadrature for the integration of the right-hand-side $e^{-tA}\mathbf{f}$. Specifically, if this integration could be carried out precisely, the scheme would then be exact, I repeat, be exact, with respect to time. Further, thanks to the $e^{-tA}$ term, the scheme is always unconditionally stable, as long as all eigenvalues of $A$ are negative (with is of course true for elliptic operators like $\Delta$).

Initial conditions for multistep methods

For multistep methods, more than one initial conditions are expected. Take $u_t=u_{xx}+f$ for instance, where $u|_t=0$ is given.

In this case, Taylor expansion would be helpful: $$ u^1\approx u({\rm d}t)=u(0)+u_t(0){\rm d}t+\frac{1}{2}u_{tt}(0){\rm d}t^2+\cdots=u^0+\left(u_{xx}(0)+f(0)\right){\rm d}t+\frac{1}{2}u_{tt}(0){\rm d}t^2+\cdots, $$ where $$ u_{tt}=u_{xxt}+f_t=\left(u_t\right)_{xx}+f_t=u_{xxxx}+f_{xx}+f_t, $$ for which $$ u_{tt}(0)=\left(u^0\right)_{xxxx}+f_{xx}(0)+f_t(0). $$ Using this trick, you may get $u^1$, $u^2$, ..., $u^k$, as many initial terms as you expect, with arbitrary order of accuracy you want. Plus, since these are just initial terms, there is no need to worry about the stability (stability matters only if it is applied to all time steps iteratively; if only for one or a few steps, it would ruin the result).

  • $\begingroup$ Hey, thanks a lot. From what you worte above, I think that a fourth-order compact scheme would be the most straightforward for me to implement. Do you know where I can find more information on this? $\endgroup$ – D.B. May 16 '18 at 15:06
  • $\begingroup$ Sure. You may refer to Section 3.6 in this reference. $\endgroup$ – hypernova May 16 '18 at 16:04

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.