A numerical boundary conditions paradox

Question

For $(t,z)\in[0,1]\times[-1,0]$

zmin = -1; tmax = 1;

and some fields $w(t,z)$ and $y(t,z)$

n = 100; h = -zmin/(n-1);
W[t_] = Table[w[i][t], {i, n}]; 
Y[t_] = Table[y[i][t], {i, n}];

let there be the following PDE's system $$\partial_tw=\partial_zy+w\partial_zw$$ $$\partial_ty=\partial_zw+w\partial_zy$$

For the implementation of the Method of Lines derivatives $\partial_zw$ and $\partial_zy$ are numerically approximated as

Wz[t_] = Join[{(w[2][t] - w[1][t])/h}, 
Table[(w[i + 1][t] - w[i - 1][t])/(2h), {i, 2, n - 1}], {(
w[n][t] - w[n - 1][t])/h}];

and

Yz[t_] = Join[{(y[2][t] - y[1][t])/h}, 
Table[(y[i + 1][t] - y[i - 1][t])/(2h), {i, 2, n - 1}], {(
y[n][t] - y[n - 1][t])/h}];

Notice that the above derivative formulas change for $i=1$ (i.e. $z=-1$) and $i=100$ (i.e. $z=0$). This is a way to handle the fact that numerical integration for $z$ is confined in $[-1,0]$ and does not imply any boundary condition.

Then the above PDE's system can be written as

wall[t_] = Yz[t] + W[t]*Wz[t];
eqw = Thread[
D[W[t], t] == wall[t] - PadLeft[{ wall[t][[n]]}, n]];
eqy = Thread[D[Y[t], t] == Wz[t] + W[t]*Yz[t]]; 

The only boundary condition that is implied by the above equations is that $$w(t,0)=0$$ and this is the reason for the cumbersome statement of the dynamical equation for $w$ ( further explanation: press ctrl+F and type "here is the answer for point 3").

The boundary condition is accompanied by the following initial conditions

w0[z_] = -0.01*Sin[z*Pi]^2;
y0[z_] = 1;
initw = Thread[W[0] == Table[w0[zmin + (i-1)*h], {i, n}]];
inity = Thread[Y[0] == Table[y0[zmin + (i-1)*h], {i, n}]];

and then NDSolve is called to implement the method of lines

lines = NDSolveValue[{eqw, eqy, initw, inity}, {W[t], Y[t]}, {t, 0, 
tmax}];

So there arise the following questions:

1. Except $w(t,0)=0$ is there any other boundary condition implicit in the finite difference equations? If it does then which? If it doesn't then why does the code run? The problem seems underdetermined.

2. Can one call the Method of Lines as internal routine so as to increase the accuracy of the above code?

I am working on these but would appreciate any help.

Answer 1

The solution you've observed is the artifact of 1st order one-sided difference formula

$$f' (x_n)\simeq \frac{- f (x_{n}-h)+ f (x_n)}{ h}$$

for approximating the PDE at the boundary. This can be confirmed by replacing it with 2nd order one-sided formula

$$f' (x_n)\simeq \frac{f (x_{n}-2h)-4 f (x_{n}-h)+3 f (x_n)}{2 h}$$

If you're not familiar with one-sided formula, start from page 6 of this book.

zmin = -1; tmax = 1;

n = 100; h = -zmin/(n - 1);
W[t_] = Table[w[i][t], {i, n}]; 
Y[t_] = Table[y[i][t], {i, n}];


help[var_] := With[{w = var}, Join[{-{1, -4, 3}.{w[3][t], w[2][t], w[1][t]}/(2 h)}, 
   Table[(w[i + 1][t] - w[i - 1][t])/(2 h), {i, 2, 
     n - 1}], {{1, -4, 3}.{w[n - 2][t], w[n - 1][t], w[n][t]}/(2 h)}]]

Wz[t_] = help@w ;
Yz[t_] = help@y;


wall[t_] = Yz[t] + W[t]*Wz[t];
eqw = Thread[
   D[W[t], t] == wall[t] - PadLeft[{ wall[t][[n]]}, n]];
eqy = Thread[D[Y[t], t] == Wz[t] + W[t]*Yz[t]]; 

w0[z_] = -0.01 Sin[z π]^2;
y0[z_] = 1;
initw = Thread[W[0] == Table[w0[zmin + (i - 1)*h], {i, n}]];
inity = Thread[Y[0] == Table[y0[zmin + (i - 1)*h], {i, n}]];

lines = NDSolveValue[{eqw, eqy, initw, inity}, {W[t], Y[t]}, {t, 0, 
    tmax}];

{testw, testy} = 
 ListInterpolation[
    Developer`ToPackedArray@#[[0]]["ValuesOnGrid"] & /@ # // 
     Transpose, {#[[1, 0]]["Coordinates"][[1]], Array[# &, n, {zmin, 0}]}] & /@ lines

Plot3D[testw[t, z], {t, 0, tmax}, {z, zmin, 0}, AxesLabel -> {t, z, f}]

But if we use 1st order one-sided formula instead:

help[var_] := With[{w = var}, Join[{(w[2][t] - w[1][t])/h}, 
   Table[(w[i + 1][t] - w[i - 1][t])/(2 h), {i, 2, n - 1}], {(
      w[n][t] - w[n - 1][t])/h}]]

The solution will be

The difference is obvious i.e. the solution depends on how we approximate the differential term at the boundary!

Further check by varying n shows both solutions are stable. This behavior never shows up when b.c. is enough AFAIK. For example, when dealing with the initial-boundary value problem

tend = 1/10; xl = 0; xr = 1;
With[{u = u[t, x]}, eq = D[u, t] == D[u, x, x];
 ic = u == Exp[-100 (x - (xl + xr)/2)^2] /. t -> 0;
 bc = {u == 0 /. x -> xl, D[u, x] == 0 /. x -> xr};]

sol = NDSolveValue[{eq, ic, bc}, u, {t, 0, tend}, {x, xl, xr}]

Both 1st and 2nd order approximation for the b.c. lead to the same solution, when the grid is dense enough:

Clear@dx
formula = eq /. {D[u[t, x], t] -> u[x]'[t], 
    D[u[t, x], x, x] -> (u[x - dx][t] - 2 u[x][t] + u[x + dx][t])/dx^2};

points = 25;
dx = (xr - xl)/(points - 1);
ode = Table[formula, {x, xl + dx, xr - dx, dx}];
odeic = Table[ic /. u[t_, x_] :> u[x][t] // Evaluate, {x, xl, xr, dx}];

bcnew1 = bc[[1]] /. u[t_, x_] :> u[x][t];

bcnew2 = bc[[2]] /. 
   D[u[t, x_], x_] :> (u[x - 2 dx][t] - 4 u[x - dx][t] + 3 u[x][t])/(2 dx);
bcnew3 = bc[[2]] /. D[u[t, x_], x_] :> (- u[x - dx][t] + u[x][t])/(dx);
mid[bc_] := (sollst = 
   NDSolveValue[{ode, odeic, bcnew1, bc}, 
    u /@ Array[# &, points, {xl, xr}], {t, 0, tend}];
  ListInterpolation[
     Developer`ToPackedArray@#["ValuesOnGrid"] & /@ # // 
      Transpose, {#[[1]]["Coordinates"][[1]], Array[# &, points, {xl, xr}]}] &@sollst)
soltest1 = mid[bcnew2];
soltest2 = mid[bcnew3];
Manipulate[Plot[{soltest1[t, x], soltest2[t, x]}, {x, xl, xr}, 
  PlotStyle -> {Automatic, {Thick, Dashed}}, PlotRange -> {0, 2}], {t, 0, tend}]

OK, then how to explain this behavior? Is the difference formula actually equivalent to a hidden b.c.? This is exactly what I've asked in this and this post, but sadly nobody has found a satisfactory answer so far.

Answer 2

@bbgodfrey commented above that boundary conditions (b.c.) are needed for both $y$ and $w$. These b.c. can be time-independent, i.e. of the form $y_b=const.$ and $w_b=const.$ But it is reasonable that b.c. could be time-independent as well, i.e. of the form $y_b=y_b(t)$ and $w_b=w_b(t)$. Fuctions $y_b(t)$ and $w_b(t)$ can have an analytic expression, but it is conceivable that this won't be always the case. They could be solutions of some differential equation that does not have any analytic solution. And in particular, I think here comes the answer, functions $y_b(t)$ and $w_b(t)$ could be solutions of the given PDE system.

In short I think that the code I posted gives a solution that by construction is characterised by b.c. $y(t,-1), y(t,0), w(t,-1)$ that are solutions of the PDE's system in question. These are complemented by b.c. $w(t,0)=0$ and together with initial conditions (i.c.) lead to a unique solution-the one to which 1st and 2nd order approximation mentioned by @xzczd converge as grid density increases.

I am not a methematician so I cannot be entirely sure that my answer is correct,though apparently reasonable. An expert's confirmation would be important.

Also if the answer is correct here comes the question if one can call method of lines as an internal routine and therefore avoid explicit discretization of the problem.