Periodic boundary conditions and finite differences implementation

I am trying to solve the following system: $$\begin{cases} \partial_t{u} + \partial_x{v} = 0 \\ \partial_t{v} + \frac{1}{\varepsilon^2}\partial_x{u} = -\frac{1}{\varepsilon^2}(v-f(u)) \end{cases}$$ with periodic boundary conditions and $$u(x,0)=\sin{(2\pi x)}e^{-x^2}$$, $$v(x,0)=0$$.

I want to use an IMEX Runge Kutta approach as described in https://arxiv.org/pdf/1701.04370.pdf on page 6*, A unified IMEX Runge-Kutta approach, where for our particular case $$p(u) = u$$ and $$f(u) = x^2 u$$.

In order to solve it numerically and impose periodic boundary conditions, I have overlapped the first and the last real knots, i.e. those coinciding with the boundaries of the physical domain $$[0,6]$$ where I want to solve the equation.

Our enumeration of the knots goes from $$x_1=0$$ to $$x_n=6$$, and in order to implement centered finite differences, we used ghost knots $$x_0 = x_{n-1}$$ and $$x_{n+1} = x_2$$.

Here is the chunk of code related to the update of the solution:

We've used $$n=400$$ discretization space steps and as time step we have chosen $$dt = dx/n$$, where $$dx=6/(n-1)$$, and solved it until $$T=0.1$$.

clear all
close all

m = 400; %space steps
x = linspace(0,6,m)';
T = .1;
dx = x(3)-x(2);
dt = dx/m; %pay attention to CFL condition
n = floor(T/dt)+1;

%% Caso 1
epsilon=0.5;

delta = .5;
u0 = sin(2*pi*x).*exp(-x.^2/(2*delta)); %initial datum
v0 = zeros(m,1);
U = zeros(2*m,n);
U(:,1) = [u0;v0];

t =  0;
i = 1;
while t+dt<T
u = U(1:m,i);
v = U(m+1:end,i);

v(1:end-1) = epsilon^2./(epsilon^2+dt).*v(1:end-1)-dt./(epsilon.^2+dt).*(([u(2:end-1);u(1)]-[u(end-1);u(1:end-2)])/(2*dx)-x(1:end-1).^2.*u(1:end-1));
v(end) = v(1);

u(1:end-1) = u(1:end-1) - dt*(epsilon^2./(epsilon^2+dt).*([v(2:end-1);v(1)]-[v(end-1);v(1:end-2)])/(2*dx)-dt./(epsilon^2+dt).*([u(2:end-1);u(1)]-2*u(1:end-1)+[u(end-1);u(1:end-2)])/dx^2 +dt./(epsilon^2+dt).*(x(1:end-1).^2.*([u(2:end-1);u(1)]-[u(end-1);u(1:end-2)])/(2*dx)));
u(end) = u(1);

U(1:m,i+1) = u;
U(m+1:end,i+1) = v;
i = i + 1;
t = t + dt;

end

plot(x,u,'r-o',x,v,'b-o','MarkerSize',1)
xlabel('x')
ylabel('u(x,t*)')
legend('u','v')
title(sprintf('Time= %0.3f',t));

My main questions are:

• Is it ok to overlap the knots and also values of the solutions at the knots in order to impose periodic BCs?
• Doing in this way some numerical oscillation arises near the boundaries, is it a problem of the scheme we've used or something is wrong in our implementation? EDIT ---------------------

We've also tried to solve the problem with homogenous Neumann BCs to overcome the problem of varying relaxation parameter $$\varepsilon = \varepsilon(x)$$ and the code is below:

clear all
close all

m = 1000; %space steps
x = linspace(0,6,m)';
delta = .5;

u0 = cos(pi*x); %initial datum neumann
v0 = zeros(m,1);
T = .1;
dx = x(3)-x(2);

%% Case 1
epsilon =0.5*ones(m,1);

%% Case 2
% epsilon = (x<3)+0.1*(x>=3);

dt=dx/100;

n = floor(T/dt)+1;
t =  0;
i = 1;
u = u0;
v = v0;

while t+dt<T
v_old = v;
v(1:end) = epsilon(1:end).^2./(epsilon(1:end).^2+dt).*v(1:end)-dt./(epsilon(1:end).^2+dt).*(([u(2:end);u(end-1)]-[u(2);u(1:end-1)])/(2*dx)-(x(1:end).^2).*u(1:end));

u(1:end) = u(1:end) - dt*(epsilon(1:end).^2./(epsilon(1:end).^2+dt).*([v_old(2:end);v_old(end-1)]-[v_old(2);v_old(1:end-1)])/(2*dx) -...
dt./(epsilon(1:end).^2+dt).*([u(2:end);u(end-1)]-2*u(1:end)+[u(2);u(1:end-1)])/(dx^2) +...
dt./(epsilon(1:end).^2+dt).*( 2*x.*u + (x(1:end).^2).*([u(2:end);u(end-1)]-[u(2);u(1:end-1)])/(2*dx)));
i = i + 1;
t = t + dt;

end
plot(x,u,'r-',x,v,'b-','MarkerSize',1)
legend('u','v')

The problem is that it still oscillates at the boundary, even with a lot of discretization points. We have used ghost nodes and we hence substituted in the code $$u_{n+1}=u_{n-1}$$ and $$u_0=u_2$$. Another problem is that even by doing this, at the boundary the slope is not $$0$$.

EDIT 2 --------- Here there is the plot of the 2 solutions at time $$T=1$$ for a non-homogeneus $$\varepsilon(x)=\begin{cases} 1,\quad x<3\\ 0.1,\quad x\geq 3 \end{cases}$$ and with $$u_0 = \sin{4\pi(x-3)}e^{-4(x-3)^2}$$, $$v_0=0$$, periodic boundary conditions and $$f(u) = (x^2-6x)u$$.

I think this is a good result since the discontinuity of the relaxation parameter $$\varepsilon$$ emerges in the final plot but it does not explode. By the way, as highlighted by @Ruslan the solution is sensitive to space discretization step and by refining it we get the second image, while with a coarser mesh we get the first plot below.  *or on page 2090 of the published article (1)

(1) S. Boscarino, L. Pareschi, G. Russo (2017): "A unified IMEX Runge--Kutta approach for hyperbolic systems with multiscale relaxation", SIAM J. Numer. Anal. 55(4), 2085-2109. doi:10.1137/M1111449

• What are the 'knots' you are referring to? Also, the boundary issue might be due to the discretisation being too large (so just make the time and spatial steps smaller whilst still adhering to the CFL condition), or possibly because the boundary conditions are being enforced when they may not actually be physically reasonable conditions, or it might be because your problem has boundary layers in which case you might want to look at the asymptotics in different regions of the problem. See here for an introduction. – mattos Aug 1 at 18:06
• @Mattos, looking at the code, the OP solved the system from knot $x_1$ to knot $x_{N-1}$ and then he imposed periodic bd. conditions by writing u(end)=u(1), v(end)=v(1). The derivatives are approximated by centered second order finite differences using a ghost node. – VoB Aug 1 at 18:32
• @VoB Thanks, I had just never heard of the nodes being called knots. – mattos Aug 2 at 2:33

I get similar oscillations (and a couple of warnings of slow convergence) when I try passing your problem to Wolfram Mathematica's NDSolve function, with the following code:

ϵ = 0.5;
T = 0.1;
x0 = 0;
x1 = 6;
sol = NDSolve[{D[u[x, t], t] + D[v[x, t], x] == 0,
D[v[x, t], t] + 1/ϵ^2 D[u[x, t], x] == -1/ϵ^2 (v[x, t] - x^2 u[x, t]),
u[x, 0] == Sin[2 π x] Exp[-x^2],
v[x, 0] == 0,
u[x0, t] == u[x1, t],
v[x0, t] == v[x1, t]
}, {u, v}, {x, x0, x1}, {t, 0, T},
MaxStepSize -> 0.01]
Plot[{u[x, T], v[x, T]} /. sol // Evaluate, {x, x0, x1}, PlotRange -> All] This should mean that it's likely the problem itself (equations + BC) is ill-defined, not the way you implemented the solver. And indeed, as we will see below, the solution "wants" to disobey your boundary conditions, so it tries hard to be discontinuous, resulting in the oscillations you get.

Let's consider what the second equation tells us for $$t=0$$:

$$\partial_t v(x,t)=\frac1{\varepsilon^2}(x^2u(x,t)-v(x,t)-\partial_x u(x,t)).$$

With our initial conditions this becomes:

$$\partial_t v(x,t)\Bigg|_{t=0}=\frac1{\varepsilon^2}\exp(-x^2)(-2\pi\cos(2\pi x) + x (2 + x) \sin(2 \pi x)).$$

This function is strictly negative at $$x=0$$ and close to zero at $$x=6$$. Discontinuity is obvious, and this will result in poor performance of the solver, since further steps will depend on $$\partial_x v(x,t)$$, which will give you Dirac-delta singularity.

Now, what happens when you take the initial conditions for your version with Neumann BC? Similarly, you'll get such $$\partial_t v(x,t)$$ at $$x=0$$, that it'll have nonzero first spatial derivative at $$x=0$$, again going against your Neumann boundary conditions.

In general, you can't get your system to behave with the same boundary conditions at $$x=0$$ for $$u$$ and $$v$$, since if one function is e.g. even, the other wants to be odd. So you may have some luck using Neumann conditions on e.g. $$u$$ and Dirichlet ones for $$v$$, or vice-versa.

I've tried instead imposing the boundary conditions (tried both Dirichlet and periodic, with similar success) at $$x=-6$$ and $$x=6$$ instead of $$x=0$$ and $$x=6$$ — the points where the initial conditions are very close to zero. In this case there are no numerical problems. But this may be unphysical, I don't know, since I don't know what the underlying physical model is.

Your actual problem, since you're trying vastly different boundary conditions on a problem which is sensitive to them, seems to be that you don't really know what conditions should be used. This information should come from physics of the system you're trying to solve. The physical model must dictate the boundary conditions, not the try-and-guess approach.

Regarding the title of your question combined with the bounty description, I'd suggest the following. Since it seems that you're somewhat confident that you need periodic boundary conditions, an easy way to check that your solver handles them correctly is to cyclically shift the domain (à la fftshift). Namely, you can replace all the $$x$$ variables in your functions (i.e. in initial conditions, $$f$$ etc.) with something like $$\operatorname{mod}(x, 6)$$ and then change the endpoints from $$0$$ and $$6$$ to e.g. $$3$$ and $$9$$. Differentiation with periodic BCs is invariant with respect to such uniform cyclic shift, so if your solver handles the boundary conditions correctly, then the solution shouldn't change at all (up to the shift). Otherwise, you should look for mistakes.

• @VoB that totally depends on the physics underlying your equations. It should suggest how to interpret truncation of the domain at $x=0$, and thus what boundary conditions to impose. One way, if you want to preserve the behavior you had with the symmetric domain, is to set boundary conditions on $u$ and $v$ of different "parity": one Dirichlet, another Neumann (matching the initial conditions, of course). This will force the solutions to be (anti-)symmetric with respect to reflection from $x=0$ (since the weights in the equations are symmetric). – Ruslan Aug 7 at 15:00
• @VoB no, the solution should be free to change values at the boundaries (unless physics says otherwise). I was suggesting to set $\partial_x u(0,t)=0$ and $v(0,t)=0$ for this case (and whatever you want at the right border). I.e. Neumann boundary for $u$ (consistent with the initial condition) and Dirichlet for $v$. – Ruslan Aug 7 at 16:49
• @VoB discontinuity is only present with periodic boundaries, otherwise it's continuous. I suppose that the actual reason for this is that it's only for initial state that the BC I suggested is good. If you try with the symmetric boundaries, you'll see $v$ not actually cross $x=0$ as time increases, instead having zero a bit off the origin, so the problem may be the same as originally, just maybe a bit pacified by the "better" BCs. – Ruslan Aug 7 at 20:00
• @VoB why does it make you wonder? You are no longer using periodic BCs, aren't you? – Ruslan Aug 7 at 20:42
• @Ruslan By shifting and contracting the initial periodic condition defined on $[-6,6]$ for $u$ in such a way that it behaves well even in $[0,6]$ I've managed to get good results even in the asymmetric case when $\varepsilon$ is constant. Now my problem is how to manage the non-homogeneus $\varepsilon(x)$ case with periodic boundary conds on both $u$ and $v$. My attempt is to solve the problem separately in the 2 intervals where the $\varepsilon$ and imposing at the end the periodicity by saying $u(end)=u(1)$ and $v(end)=v(1)$. Could it be a good approach? – Dadeslam Aug 8 at 9:56

I have tried my own Matlab implementation, and similar results are obtained:

clear all
close all

% parameters
m = 400;
x = linspace(0,6,m);
dx = x(2) - x(1);
x = [x(1)-dx, x, x(end)+dx];
Co = 0.95;  % Courant number
Fo = 0.45;  % Fourier number
T = 0.1;
epsi = 0.5;

%initial datum
delta = .5;
u0 = sin(2*pi*x).*exp(-x.^2/(2*delta));
v0 = 0*x;

% initialization
t = 0;
u = u0;
v = v0;
dtHyp = Co*dx/max(x.^2);
dtPar = 0.5*Fo*dx^2*(1 + sqrt(1 + 4*epsi^2/(Fo*dx^2)));
dt = min([dtHyp, dtPar]);

figure;
clf;
hp = plot(x,u);
ht = title(strcat('t = ',num2str(t)));
xlim([x(1) x(end)]);
ylim([-1 1]);

%iterations
while (t+dt<T)
up = circshift(u, -1);
um = circshift(u, +1);
vp = circshift(v, -1);
vm = circshift(v, +1);
xp = circshift(x, -1);
xm = circshift(x, +1);

k = dt/(epsi^2+dt);

% scheme
v = v - k*(up-um)/(2*dx) - k*(v-x.^2.*u);
u = u - k*epsi^2*(vp-vm)/(2*dx) - k*dt*(xp.^2.*up-xm.^2.*um)/(2*dx) + k*dt*(up-2*u+um)/dx^2;

% boundary conditions
u(1) = u(end-1);
v(1) = v(end-1);
u(end) = u(2);
v(end) = v(2);

t = t + dt;

dtHyp = Co*dx/max(x.^2);
dtPar = 0.5*Fo*dx^2*(1 + sqrt(1 + 4*epsi^2/(Fo*dx^2)));
dt = min([dtHyp, dtPar]);

set(hp,'YData',u);
set(ht,'String',strcat('t = ',num2str(t)));
drawnow;
end I believe that it is caused by the fact that periodic boundary conditions are applied to a non-periodic problem. Indeed, the spatially-varying flux $$f(u) = x^2 u$$ leads to a discontinuity in speed of sound $$c = f'(u)$$ --which appears in the CFL condition-- at the domain boundaries. Thus, a (fictional) interface between neighbor domains is created.

In the case of periodic BCs, if the flux $$f(u) = u$$ is chosen instead, then everything seems to work well. Tackling such fluxes as $$f(u) = x(x-6) u$$ does not seem impossible with the prescribed method (no velocity jump at the boundaries). Nevertheless, one may consider using dedicated methods for spatially-varying fluxes. For this purpose, one may have a look at the article (1) and references therein (even if the article does not give much information on IMEX methods).

Alternatively, one could consider implementing outflow boundary conditions instead

% boundary conditions
u(1) = u(2);
v(1) = v(2);
u(end) = u(end-1);
v(end) = v(end-1);

(1) D.I. Ketcheson, M. Parsani, R.J. LeVeque (2013): "High-order wave propagation algorithms for hyperbolic systems", SIAM J. Sci. Comput. 35(1), A351-A377. doi:10.1137/110830320

• Because by setting $\varepsilon = \varepsilon(x)$ = 1 for $x\in[0,3]$ and 0 for $x\in(3,6]$ I still get oscillations. – Dadeslam Aug 2 at 12:09
• Yes, it is another problem. But do you think this variation of the problem brings oscillations because of how we discretized the domain (hence CFL) or because of the discontinuity of the $\varepsilon(x)$ function? – Dadeslam Aug 2 at 12:35
• @Harry49 I'm working with dadeslam on this problem. Could it be wrong the way in which we approximate the spatial derivatives? We approximated $u_t$ and $v_t$ forward in time, while use centred FD in the space. – VoB Aug 2 at 16:18
• @Harry49 Okay, but do you think that this could be a problem? Because in the scheme written in the paper we do not have much indication on how to discretize derivatives – VoB Aug 5 at 9:09
• @Harry49 I've just updated the question, because we've tried to implement homogeneus Neumann boundary conds to solve the problem of varying relaxation paramenter but it still oscillates even for a constant paramenter $\varepsilon$. – Dadeslam Aug 6 at 16:19