Method of Lines: How to simplify Jacobian with periodic BCs?

Question

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}}{2Delta x}.$$
Written in matrix form this gives
$$frac{d}{dt}begin{pmatrix}
u_0 
u_1 
... 
u_{N-1}
end{pmatrix}
=
-frac{1}{2Delta 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}{2Delta 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?

Steven Roberts · Accepted Answer

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}{2Delta 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}{2Delta 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.

Method of Lines: How to simplify Jacobian with periodic BCs?

One Answer

Add your own answers!

Ask a Question