# A numerical boundary conditions paradox

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];
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.

• 1. h = -zmin/n should be h = -zmin/(n-1). 2. /h, {i, 2, n - 1}] should be /(2h), {i, 2, n - 1}]. 3. What's the meaning of PadLeft[{ wall[t][[n]]}, n]? If you think this will impose $w(t,0)=0$, you're wrong. – xzczd Dec 23 '18 at 12:37
• I believe that boundary conditions are needed for both w and y. I am unaware of any way to call Method of Lines as internal routine, but you could use NDSolve directly to solve the coupled PDEs. – bbgodfrey Dec 23 '18 at 13:36
• @xzczd here is the answer for point 3. PadLeft[{ wall[t][[n]]}, n] creates a list with n elements all of which is zero except the nth one that equals wall[t][[n]] . As a result the subtraction results in a list of n elements the first n-1 equal to Yz[t] + W[t]*Wz[t] and the last one equal to zero. This means that $\partial_tw(t,0)=0$ and results in $w(t,0)=w(0,0)=\sin(0)^2=0$. – dkstack Dec 23 '18 at 13:36
• @xzczd I edited my question according to point 2. – dkstack Dec 23 '18 at 13:47
• @bbgodfrey boundary condition about y concerns its value or could be stated in terms of its z-derivative? Also can you see why the code runs I mean if there is an implicit boundary condition in my code? – dkstack Dec 23 '18 at 13:50

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];
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[
DeveloperToPackedArray@#[[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[
DeveloperToPackedArray@#["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.

• xzczd can you help me with this question? – dkstack Jan 9 at 8:42
• xzczd, what do you think of this post? It may lead to a general method for first order IBVP problems that mix spatial and temporal derivatives. – dkstack Jan 17 at 16:31

@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.

• "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." Notice that given the highest differential order of $w$ and $y$ respect to $z$ is $1$, usually you only need $2$ boundary conditions. (Related: math.stackexchange.com/q/450367/58219) There exist more involved cases e.g. Inviscid Burgers equation though. – xzczd Dec 26 '18 at 14:29
• Ok then forget about $y(t,-1)$ and $w(t,-1)$. My point is that the solution that occurs is not an artifact but a fare and square solution with b.c. that are time dependent, without any analytic expression, and which by construction are numerical solutions of the above PDE system. – dkstack Dec 26 '18 at 14:34
• Notice I've never denied that lines etc. are solutions for the PDE. I use the word "artifact" because those solutions aren't determined before a certain one-sided formula is chosen. – xzczd Dec 26 '18 at 15:10
• Yes but you have proven that any one-sided formula in the limit of very dense grid points to the same solution. – dkstack Dec 26 '18 at 15:53
• So we can reverse our logic and say: yes b.c. can be time dependent. Yes time dependent b.c. may be solutions of equations that do not have analytical solutions. And yes time dependent b.c. can be solutions of the specific PDE system we try to solve. – dkstack Dec 26 '18 at 15:56