0
$\begingroup$

I want to implement the forward in time, centered in space scheme for the linear advection \begin{align} u_t+ a u_x=0 \end{align} with periodic boundary conditions and initial datum $u(x,0)$. I know that this scheme is unconditionally unstable and the theory about that, but my question is just on the implementation part

I'm doing this in the interval $x \in [0,1]$ and $t \in [0,1]$.

The scheme is

\begin{align} u_{i}^{n+1}= u_{i}^n-\frac{a dt}{2dx}(u_{i+1}^n-u_{i-1}^n) \end{align}

Say I discretize $[0,1]$ with $M$ points. So $dx=\frac{1}{M-1}$. My main problem is the implementation of the boundary conditions. I overlap the nodes $x_1=0$ and $x_M=1$.

So, at the first iteration, my scheme is (I omit the time index)

\begin{align} u_1=u_1- \frac{a dt}{2 dx} (u_2-u_0) \end{align}

Now, since $u_0$ is not known, I would use the fact that $x_1=x_M$, and hence the point at the left of $x_1$ is $x_{M-1}$ and hence $u_0=u_{M-1}$.

Now I would solve this up to node $x_{M-1}$, so I will have

\begin{align} u_{M-1}=u_{M-1} - \frac{a dt}{2 dx} (u_{1} - u_{M-2}) \end{align} where I used the fact that $u_M=u_1$.

Then I will update $u_{M}=u_{1}$.

$\endgroup$
8
  • $\begingroup$ @Harry49 If my implementation of the periodic boundary conditions is right and if there's anything wrong in what I've written, since indices are making me crazy $\endgroup$ Commented Oct 7, 2019 at 16:17
  • $\begingroup$ Many thanks! @Harry49 $\endgroup$ Commented Oct 7, 2019 at 17:52
  • $\begingroup$ The main problem is that I've seen different implementation, and I still don't understand how to treat the fact that $x_1=x_M$ (i.e. overlapping first and last node) $\endgroup$ Commented Oct 7, 2019 at 19:42
  • $\begingroup$ Why do you solve up to the point $x_{M-1}$? It should be $x_M$. And the point $i=M-1$, you should have $u_{M-1}=u_{M-1} - \frac{a dt}{2 dx} (u_{M} - u_{M-2})$. But the point that you interested is $i=M$, in this case: we have $u_{M}=u_{M} - \frac{a dt}{2 dx} (u_{M+1} - u_{M-1})$ but $u_{M+1} = u_{2}$ $\endgroup$ Commented Oct 7, 2019 at 20:02
  • $\begingroup$ @Sesame I solve up to $M-1$ and the impose $u_M=U_1$. You're right for the expression for $u_{M-1}$, it was a typo, now I fixed it $\endgroup$ Commented Oct 7, 2019 at 20:08

1 Answer 1

Reset to default
0
$\begingroup$

We give the following space discretization : $\forall i\in\{1,\dots,M\}$: $\{x_1=0,x_2,\dots,x_M=l\}$ (I started to count from 1 on purpose as you code in Matlab...). Suppose that we want to solve the above PDE by using an Euler scheme: \begin{align} u_{i}^{n+1}= u_{i}^n-\gamma u_{i+1}^n + \gamma u_{i-1}^n \end{align} where $\gamma = \frac{a dt}{2dx}$. We impose periodic boundary conditions, i.e. $u(0,t)=u(l,t)$. This means $\forall n$, \begin{align} u_0^n&=u^n_M\\ u_{M+1}^n &= u_1^n\\ \end{align} Let's look at know what happen at $i=1$ and $i=M$: \begin{align} i=1 , \quad u_{1}^{n+1} &= u_{1}^n-\gamma u_{2}^n + \gamma u_{0}^n \\ &= u_{1}^n-\gamma u_{2}^n + \gamma u_{M}^n \\ i=M, \quad u_{M}^{n+1} &= u_{M}^n-\gamma u_{M+1}^n + \gamma u_{M-1}^n \\ &= u_{M}^n-\gamma u_{1}^n + \gamma u_{M-1}^n \\ \end{align} Writting the above in matrix form, you have for $n\in \{1,\dots, N\}$ : \begin{equation} \mathbf{U}^{n+1} = \mathbf{A}\mathbf{U}^{n} \end{equation} Which can rewritten as for $M=5$ \begin{align} \begin{pmatrix} u_1^{n+1}\\ u_2^{n+1}\\ u_3^{n+1}\\ u_4^{n+1}\\ u_5^{n+1}\\ \end{pmatrix}= \begin{pmatrix} &1 &-\gamma &0 &0 &\gamma \\ &\gamma &1 &-\gamma &0 &0 \\ &0 &\gamma &1 &-\gamma &0 \\ &0 &0 &\gamma &1 &-\gamma \\ &-\gamma &0 &0 &\gamma &1 \\ \end{pmatrix}\begin{pmatrix} u_1^n\\ u_2^n\\ u_3^n\\ u_4^n\\ u_5^n\\ \end{pmatrix} \end{align}

$\endgroup$
7
  • $\begingroup$ Thanks, this makes sense. Just one more question. With this discretization, is it true that $u_1=u_M$? $\endgroup$ Commented Oct 7, 2019 at 20:54
  • $\begingroup$ Yes. I think in last row I have a coefficient of 1 for $u_1$. So that $ u_{M}^{n+1} = u_{1}^n-\gamma u_{2}^n + \gamma u_{M-1}^n $ $\endgroup$ Commented Oct 7, 2019 at 21:01
  • $\begingroup$ Sorry but I do not agree, because that expression does not imply that $u_{M}=u_1$. Moreover, in your matrix the $1$ should be in entry (m,m), as the scheme told us $\endgroup$ Commented Oct 7, 2019 at 21:14
  • $\begingroup$ @VoB: Sorry I was not thinking straight on this one. I made some modifications at my initial post. Hope it is clearer now. $\endgroup$ Commented Oct 8, 2019 at 0:03
  • $\begingroup$ so, the periodicity is imposed by using a value outside of the domain, i.e. $u_0$, not imposing that the first and last value nodes $u_1,u_M$ are equal, right? $\endgroup$ Commented Oct 8, 2019 at 7:23

You must log in to answer this question.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.