# Upwind finite difference: Matrix Implementation

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

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.

• 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. – ABCCHEM Nov 24 '19 at 16:54
• 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. – cfdlab Nov 25 '19 at 2:32
• 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. – Lutz Lehmann Nov 27 '19 at 10:46