2
$\begingroup$

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?

$\endgroup$
  • $\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$ – Abdullah Ali Sivas Sep 1 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$ – Maxim Umansky Sep 1 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$ – Abdullah Ali Sivas Sep 1 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$ – Maxim Umansky Sep 1 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$ – Abdullah Ali Sivas Sep 1 at 17:45
3
$\begingroup$

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.

| cite | improve this answer | |
$\endgroup$
  • 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$ – Abdullah Ali Sivas Sep 1 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$ – Maxim Umansky Sep 1 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 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$ – Maxim Umansky Sep 1 at 14:38

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.