2
$\begingroup$

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.

Improve this question
$\endgroup$

1 Answer 1

Reset to default
2
$\begingroup$

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.

Improve this answer
$\endgroup$
3
  • $\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$
    – ABCCHEM
    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$
    – Lutz Lehmann
    Commented Nov 27, 2019 at 10:46

Not the answer you're looking for? Browse other questions tagged or ask your own question.