Handling the domain's edge with central differences

Question

Let there be the following IBVP $$ \text{Domain:}\quad (t,x)\in[0,1]\times[0,1] $$ $$ \text{IC:}\quad y(0,x)=x $$ $$ \text{BC:}\quad y(t,0)=0 $$ $$ \text{PDE:}\quad\partial_t y (t,x) = \partial_x y(t,x) $$ Suppose one chooses to address the problem with FDM with $100$ points for the x-grid the $i=1$ point being at $x=0$ and the $i=100$ at $x=1$.

In this case numerical integration requires that the $x$-derivative $\partial_x$ is approximated on every grid point except for $i=1$ where $y(t,x)$ is time-independent.

So one needs a finite different approximation for $\partial_x y$ at $i=100$. This is readily done if one chooses the backward-difference formula.

However, what if one chooses the central difference formula?

Note: It is not profound that a combination of finite different formulas would be legal. For example one may not use the forward difference formula at $i=1$ along with the backward one at $i=100$, because then numerical integration with insufficient BC would be possible.

P.S. To be more generic, my problem is that in some IBVPs, at the boundary points, one is not able to use the same finite difference formula as in the rest of the spatial grid. There are simply not enough neighboring points to do this. So one has to combine different formulas. Which is the right way to do this? Is there any bibliography?

Answer 1

I implement that problem with central difference quotients, using the backwards second order formula at $x=1$, and using a likewise second order method for the time steps, choosing 1/4th of the space step so that the time error is not dominant. The resulting python code is

def F(u):
    # at time t, u[k] approximates y((k+1)*h,t)
    du = 0*u;
    du[0] = (u[1] - 0)/(2*dx);
    du[1:-1] = (u[2:]-u[:-2])/(2*dx);
    # use second order backwards difference
    du[-1] = (3*u[-1]-4*u[-2]+u[-3])/(2*dx);
    return du

# integrate with Heun's method
N = 100;
x = np.linspace(0,1,N+1);
dx = x[1]-x[0];
u = +x[1:];
plt.plot(x[1:],u,label="t=0");
t=0; dt = dx/4;
for k in range(4):
    for j in range(N):
        du1 = dt*F(u);
        du2 = dt*F(u+du1);
        u += 0.5*(du1+du2);
        t += dt;
    plt.plot(x[1:],u,label="t=%.3f"%t);
plt.grid(); plt.legend(); plt.show()

giving the plot

As one can see, the transport resulting from the given transport equation is to the left, in principle raising the curve. However, the left boundary condition appears as the problematic one, not the right as conjectured in the question. The "arriving" values are, theoretically, steadily increasing, but you force a jump to zero from the second to the first grid point with the boundary condition.

Or in other words, the continuous problem does not have a continuous solution at $x=0$, so it would be rather unexpected if a numerical solution at that point would be regular.

---

Reversing the direction for the transport equation $∂_ty(t,x)=-∂_xy(t,x)$ gives the a regular solution, filling with the value zero from the left