I want to implement the upwind finite difference scheme for the 1D linear advection equation using a finite difference matrix in python:

$$ A =\begin{pmatrix} 1-a\cfrac{\Delta t}{\Delta x} & 0 & 0 & \cdots & a\cfrac{\Delta t}{\Delta x}\\ a\cfrac{\Delta t}{\Delta x} & 1-a\cfrac{\Delta t}{\Delta x} & 0 & \cdots & 0 \\ \vdots & \vdots& \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & 1-a\cfrac{\Delta t}{\Delta x} \end{pmatrix}$$

with $1 -a\cfrac{\Delta t}{\Delta x}$ on the diagonal and $a\cfrac{\Delta t}{\Delta x}$ on the lower subdiagonal and periodic boundary conditions as indicated by the top right entry of the matrix ($a$ being constant). To generate the velocities at the next time step $n$ I used the following equation:

$$\ U_{n+1} = AU_n $$

My idea is that given inital conditions $ U_1 $, I could apply $ A^{n+1} $ to $U_1$ to get the $U_{n+1}$:

$$\ U_{n+1} = A^{n+1} U_1 $$

I am not entirely sure if this is a valid approach as I am complete novice to the subject. Any help and hints would be greatly appreciated.


1 Answer 1


Your idea of taking matrix powers is not a good one.

You should code it for step by step updating $$ u_0 = \textrm{initial condition} $$ $$ u_0 \rightarrow u_1 \rightarrow u_2 \rightarrow $$ so at any time you at most store two levels of solution $u_n, u_{n+1}$.

Even when updating $u_n \rightarrow u_{n+1}$ I would avoid using a matrix since you are multiplying with a lot of zeroes, unless you use a sparse matrix, which is possible in Python.

  • $\begingroup$ Thank you for your reply. Just as a follow up: what should my analytical solution look like? I imagine that if my inital condition is say a Gaussian pulse, then at every new time step I should still have my initial Gauss wave, but translated along the x-direction. $\endgroup$
    Commented Nov 24, 2019 at 16:54
  • $\begingroup$ The initial wave will move but will change shape due to dissipation/dispersion errors. Also, the last row in your "A" matrix has some error I think. $\endgroup$
    – cfdlab
    Commented Nov 25, 2019 at 2:32
  • $\begingroup$ As $A=(1-p)I+pR$ where $p=\frac{aΔt}{Δx}$, you can apply the binomial theorem to get $$A^n=(1-p)^nI+n(1-p)^{n-1}pR+\tbinom{n}2(1-p)^{n-2}p^2R^2+...+n(1-p)p^{n-1}R^{n-1}+p^nR^n.$$ With $R^n=I$ this then folds back. With $p$ sufficiently small you could also approximate $$A^m=\exp(-mp(I-R))=e^{-mp}\exp(mpR),$$ which shows that you will always get dissipation of the original sharp impulse. $\endgroup$ Commented Nov 27, 2019 at 10:46

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.