Consider the advection equation $$\frac{\partial u}{\partial t}+c(x)\frac{\partial u}{\partial x}=0.$$ With periodic boundary condtitions in $x$ with period $L$, i.e. $u(x,t)=u(x+L,t)$ and initial condition $u(x,0)=f(x)$. We can solve this numerically by discretizing in $x$ to get a set of ODEs in $t$. Let $$u_i(t)=u(x_i,t),$$ for $i=0,1,...N-1$. Hence, $$\frac{du_i}{dt}=-c(x_i)\frac{u_{i+1}-u_{i-1}}{2\Delta x}.$$ Written in matrix form this gives $$\frac{d}{dt}\begin{pmatrix} u_0 \\ u_1 \\ ... \\ u_{N-1} \end{pmatrix} = -\frac{1}{2\Delta x} \begin{pmatrix} 0 & c(x_0) & 0 & 0 & ... & 0 & -c(x_{0}) \\ -c(x_1) & 0 & c(x_1) & 0 & ... & 0 & 0 \\ ... & & & & & & \\ c(x_{N-1}) & 0 & 0 & ... & 0 & -c(x_{N-1}) & 0 \\ \end{pmatrix} \begin{pmatrix} u_0 \\ u_1 \\ ... \\ u_{N-1} \end{pmatrix}.$$ We can solve this using an ODE solver e.g. solve_ivp. Note that $$A=-\frac{1}{2\Delta x} \begin{pmatrix} 0 & c(x_0) & 0 & 0 & ... & 0 & -c(x_{0}) \\ -c(x_1) & 0 & c(x_1) & 0 & ... & 0 & 0 \\ ... & & & & & & \\ c(x_{N-1}) & 0 & 0 & ... & 0 & -c(x_{N-1}) & 0 \\ \end{pmatrix}$$ gives the Jacobian matrix of the system. It is nearly tridiagonal except for the top right and bottom left corners. It would be nice if I could give a sparse matrix for the Jacobian matrix because I assume this will save computation time. Do you know any tricks to solve this system and give the Jacobian matrix which will be as computationally efficient as possible?

  • $\begingroup$ This is not immediately related to your question, but I suggest you to use upwind discretization for the $\partial u/\partial x$ term. It will not resolve circulancy, but it is a better approximation than central difference. $\endgroup$ Sep 1, 2020 at 0:51
  • $\begingroup$ @Abdullah Ali Sivas I thought upwind is the choice for explicit time-step integration because it gives you stability; for implicit it does not matter - correct? Or the issue, aside from the stability, is having a positivity-preserving finite-difference representation? $\endgroup$ Sep 1, 2020 at 4:22
  • $\begingroup$ @MaximUmansky, stability is not much of an issue. I do not exactly remember whether central difference combined with implicit time stepping is stable, but one can always hack it and add some artificial diffusion. But definitely and more importantly, it is about preventing under/overshoots and, in general, obtaining more physically relevant solutions. For example, in the case of burgers equation the advantages of using upwind over central becomes so much more clear. $\endgroup$ Sep 1, 2020 at 4:34
  • $\begingroup$ @Abdullah Ali Sivas Do you have a reference on this? I just can't believe that using a low-order scheme is necessarily better (more accurate), if stability is not an issue. $\endgroup$ Sep 1, 2020 at 5:24
  • $\begingroup$ @MaximUmansky, I don't have a reference off the top of my head - especially on accuracy. I bet when central difference works, it works better than upwind, but I remember that for purely hyperbolic flows, central difference may introduce oscillations due to under/overshoots in approximation. I will look for a reference later today, or I will prepare a demo. $\endgroup$ Sep 1, 2020 at 17:45

1 Answer 1


If your advection problem had Dirichlet of Neumann boundary conditions, the linear system would be tridiagonal and you could apply the Thomas algorithm. With periodic boundary conditions, however, we lose this. If c(x) is a constant independent of x, the matrix would be circulant and linear systems could be solved efficiently using FFTs. An even better solution exists for this periodic 1D problem, though. The key is to realize your matrix is a rank one update from a tridiagonal matrix:

$$ A=-\frac{1}{2\Delta x} \begin{pmatrix} c(x_0) & c(x_0) & 0 & 0 & ... & 0 & 0 \\ -c(x_1) & 0 & c(x_1) & 0 & ... & 0 & 0 \\ ... & & & & & & \\ 0 & 0 & 0 & ... & 0 & -c(x_{N-1}) & -c(x_{N-1}) \\ \end{pmatrix} -\frac{1}{2\Delta x} \begin{pmatrix} -c(x_0) \\ 0 \\ \vdots \\ 0 \\c(x_{N-1}) \end{pmatrix} \begin{pmatrix} 1 & 0 & \cdots & 0 & 1 \end{pmatrix} $$

Applying the Sherman–Morrison formula allows you to use the Thomas algorithm with a small amount of extra work. Asymptotically, it still costs $\mathcal{O}(N)$. The Thomas algorithm Wikipedia page has a section that describes this and also provides the algorithm for periodic BCs.

Note you can use sparse storage formats for $A$ even with the periodic BCs (e.g. spdiags in Matlab). If you had to do things manually, you would just need an $N \times 3$ matrix for the 3 diagonals where you include the corner elements of $A$ in the first and third columns.

  • 1
    $\begingroup$ I was coming to answer "FFT then you can solve this easy", but I saw your Woodbury matrix formula and I knew I had to upvote. That is pretty smart, my kudos. $\endgroup$ Sep 1, 2020 at 0:48
  • 1
    $\begingroup$ @Abdullah Ali Sivas The FFT approach would be limited to linear problems, but if $c=c(u)$ this solution shown here still works. $\endgroup$ Sep 1, 2020 at 4:24
  • $\begingroup$ Thank you. Do you know how to calculate the solution using FFTs if $c=c(x)$? They show how to use FFTs for Laplace's equation in the link below. But that is for constant coefficients. Do FFTs work if $c=c(x)$? en.wikipedia.org/wiki/…. $\endgroup$
    – Peanutlex
    Sep 1, 2020 at 10:20
  • 1
    $\begingroup$ For linear problems with variable coefficients you'd have to expand the coefficients in FFT as well. A single harmonic is a solution only for constant coefficients, for variable coefficients different harmonics are coupled. In my view, the spectral approach is usually not the best way to go, but if $c(x)$ happens to be a simple harmonic function, e.g., $c(x) = a + b \sin(x)$ that could be a good choice. You can use spectral method for nonlinear problems too; but that's not an easy way to solve it. You would need a reason why use a spectral method there, e.g., the high rate of convergence. $\endgroup$ Sep 1, 2020 at 14:38

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.