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$
    – VoB
    Oct 7 '19 at 16:17
  • $\begingroup$ Many thanks! @Harry49 $\endgroup$
    – VoB
    Oct 7 '19 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$
    – VoB
    Oct 7 '19 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$
    – Sesame
    Oct 7 '19 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$
    – VoB
    Oct 7 '19 at 20:08
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$
    – VoB
    Oct 7 '19 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$
    – Sesame
    Oct 7 '19 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$
    – VoB
    Oct 7 '19 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$
    – Sesame
    Oct 8 '19 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$
    – VoB
    Oct 8 '19 at 7:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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