1

I wish to solve the equation: $$\frac{\partial v}{\partial t}=-\frac{\partial v}{\partial x}$$ along with the boundary conditions: $$\frac{\partial v}{\partial x}\Bigg|_{x=0}=\frac{\partial v}{\partial x}\Bigg|_{x=1}$$

I want to use an $O(\delta x)^{2}$ stencil for the spacial derivative. So I have $$\frac{\partial v}{\partial x}\Bigg|_{t=t_{i},x=x_{j}}=-\frac{3}{2\delta x}v_{i,j}+\frac{2}{\delta x}v_{i,j-1}-\frac{1}{2\delta x}v_{i,j-2}$$, The stencil for my system is then: $$v_{i+1,j}=\left(1+\frac{3\alpha}{2}\right)v_{i,j}-2\alpha v_{i,j-1}+\frac{\alpha}{2}v_{i,j-1}$$ To deal with the boundary conditions I simply evaluate the governing equation at the boundary points to get: $$\frac{\partial v}{\partial t}\bigg|_{x=0,1}=0$$ which allows me to write: $$v_{i+1,1}=v_{i,1},\quad v_{i+1,N}=v_{i,N}$$ The issue then becomes at $j=2$, for this I evaluate the initial stencil at $j=2$ to obtain: $$v_{i+1,2}=\left(1+\frac{3\alpha}{2}\right)v_{i,2}-2\alpha v_{i,1}+\frac{\alpha}{2}v_{i,0}$$ To find the point at $j=0$ I use the boundary condition again in the form: $$\frac{\partial v}{\partial x}\Bigg|_{x=0}=\frac{v_{i,1}-v_{i,0}}{\delta x}=0$$ This allows me to write $v_{i,0}=v_{i,1}$, I then can write the stencil at $j=2$ as:

$$v_{i+1,2}=\left(1+\frac{3\alpha}{2}\right)v_{i,2}-\frac{3\alpha}{2}v_{i,1}$$

However when I run this, I get a numerical instability. Any idea what I'm doing wrong?

Matthew Hunt
  • 543
  • How did you choose your step sizes in space and time? – Sean Roberson May 15 '23 at 16:34
  • This was almost random, but thinking about it I should have used the Courant condition. I did the problem with $O(\delta x)$ and $O(\delta t)$ finite differences of the derivatives. I had an $\alpha$ of around $0.9$ and I chose 100 points in my spacial domain. – Matthew Hunt May 15 '23 at 16:50
  • You'll need to use that CFL condition. Or, if you can afford it, use an implicit solver for the space variables. – Sean Roberson May 15 '23 at 16:56
  • The wave speed, $c$ will just be $c=1$ won't it? so I will have to set $\delta x/\delta t=0.1$ or something won't I – Matthew Hunt May 15 '23 at 17:00
  • Apart from the CFL condition, was the way I handled the boundary conditions okay? I plan to use an implicit method, but I want to get this simple equation solved first. – Matthew Hunt May 15 '23 at 17:02
  • 1
    No, the boundary interpretation is wrong. What you wrote or cited first means periodic boundary conditions. – Lutz Lehmann May 15 '23 at 21:16
  • How should I have done it? Why was what I did wrong? – Matthew Hunt May 15 '23 at 21:22
  • Are you sure about the initial model? If the solution is unique, it would be just $v(t,x)=0$. – PierreCarre May 16 '23 at 10:38
  • If I have an initial disturbance, that disturbance should just propagate along a characteristic. – Matthew Hunt May 16 '23 at 10:48
  • Are the boundary conditions supposed to be periodic, i.e., $u(0) = u(1)$ instead of $u'(0) = u'(1)$? The first condition implies the second, but the second does not imply the first as values of $u$ above the characteristic $x=t$ are only specified up to an additive constant, meaning that the solution to the PDE is not unique (allowing for weak solutions). This may be causing this issue. Try using standard periodic boundary conditions ($v_{N+k} = v_k$) and see if this fixes things. – whpowell96 May 16 '23 at 16:32

1 Answers1

1

Using the NDSolve[] Module in Wolfram Mathematica, with the initial condition $u(0,x)=x(1-x)$, one gets the following numerical result:

enter image description here

Using your spatial discretization and denoting $u_i(t)= u(t,x_i)$, the semi-discretization in space leads to the system of ODEs $$ U'(t) = \frac{1}{2 \delta x}A U(t) $$

where

$$A= \begin{pmatrix} -3 &&&&\\ 4 & -3 &&&\\ -1 & 4 & -3 &&\\ & \ddots&\ddots & \ddots&\\ && -1 & 4 & -3 \end{pmatrix} $$

Using an implicit scheme for the time advance, for instance the trapezoidal rule, one gets

$$ U^{i+1} = U^i + \frac{\delta t}{4 \delta x}(A U^{i+1}+ AU^i) $$

i.e. $$ (I-\frac{\delta t}{4 \delta x} A) U^{i+1} = (I+\frac{\delta t}{4 \delta x}A) U^i $$

Regarding the boundary conditions, I used the very crude approach of substituting the first two equations in the linear system by $$ U_2^{i+1} - U_2^{i+1} = U_N^{i+1}-U_{N-1}^{i+1}, \quad U_2^{i+1}= \frac 12 (U_1^{i+1}+U_3^{i+1}) $$

and I got the following results

enter image description here

Even though there are perturbations introduced by the "stupid" discretisation of the boundary conditions, the scheme remains stable.

PierreCarre
  • 20,974
  • 1
  • 18
  • 34
  • The boundary conditions are meant to model "no flux" boundary conditions.

    What should happen is that your initial conditions should be transported for the right in time. So if your condition is $v(0,x)=f(x)$ for example, then the solution should be $$v(t,x)=f(x-t)$$

     – Matthew Hunt May 16 '23 at 19:25
  • @MatthewHunt The boundary conditions that you specified no not mean "no flux"... – PierreCarre May 17 '23 at 08:00
  • @MatthewHunt Taking $v(0,x)=x(1-x)$ you would have $v(t,x)=(x-t)(1-x+t)$ and, consequently, $$ \frac{\partial v}{\partial x} (t,0) = 1+2t, \quad \frac{\partial v}{\partial x} (t,1)= -1+2t$$, which is inconsistent with the boundary conditions you specified. – PierreCarre May 17 '23 at 10:49