Why do I still obtain a unique solution with one-sided formula when b.c. isn't enough?

Question

Let me illustrate the issue with an simplified example. Suppose we want to solve the following problem with finite difference method (FDM):

$$\frac{\partial u(t,x)}{\partial t}=\frac{\partial^2 u(t,x)}{\partial x^2}$$ $$u(0,x)=x(1-x),\ u(t,0)=0$$ $$t\in[0,1],\ x\in[0,1] $$

One boundary condition (b.c.) is missing, but let's ignore it for a while and discretize the equation, initial condition (i.c.) and b.c. with difference formula. Here we use 2nd order difference formula. Clearly if we only use

$$f' (x_i)\simeq \frac{f (x_{i}+h)-f (x_{i}-h)}{2 h}$$ $$f'' (x_i)\simeq\frac{f (x_{i}-h)-2 f (x_i)+f (x_{i}+h)}{h^2}$$

we cannot generate equation for $t=1,\ x=1$, so we need one-sided formula (If you're not familiar with one-sided formula, start from page 6 of this book ):

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

OK, now here comes the problem: we surprisingly find that, even if b.c. is not enough, we still obtain a closed algebraic equation system!

One may suspect the system will be kind of ill-posed e.g. "This system won't be solvable!" or "This system will have infinite many solutions!" etc., but once again, surprisingly, further check shows the system determines a unique solution and seems to be stable i.e. with smaller grid size or higher order difference formula, the solution doesn't have significant change.

How to explain this solution? If a hidden b.c. has been imposed, what's it?

At the very least, if the solution is really nothing more than numeric error, I'd like to know the exact cause of the error.

---

Here's the source code for solving the problem above in Mathematica. I've used pdetoae for the generation of difference formula:

eq = D[u[t, x], t] == D[u[t, x], x, x];
ic = u[0, x] == x (1 - x);
bc = u[t, 0] == 0;
points = 25;
grid = tgrid = Array[# &, points, {0, 1}];
ptoafunc = pdetoae[u[t, x], {tgrid, grid}, 2];
ae = Rest /@ (ptoafunc[eq] // Rest);
aeic = ptoafunc[ic] // Rest;
aebc = ptoafunc[bc];
var = Outer[u, tgrid, grid] // Flatten;
{b, m} = CoefficientArrays[{ae, aeic, aebc} // Flatten, var];
solarray = Partition[LinearSolve[N@m, -b], points];
solfunc = ListInterpolation[solarray, {tgrid, grid}];

style = Style[#, 16] &;
Plot3D[solfunc[t, x], {t, 0, 1}, {x, 0, 1}, AxesLabel -> style /@ {"t", "x"}]

Plot[solfunc[t, #] & /@ grid // Evaluate, {t, 0, 1}]

The following is the norm of solution $\sqrt{\int_0^1 u(t,x)^2 \, dx}$:

norm[t_?NumericQ, order_: 2, func_: solfunc] := 
 Power[NIntegrate[func[t, x]^order, {x, 0, 1}, 
   Method -> {Automatic, "SymbolicProcessing" -> 0}], 1/order]

Plot[norm[t], {t, 0, 1}, AxesLabel -> {t, norm}, PlotRange -> All]

Remark

1. Notice that the value of i.c. is used at the corner $t=0,\ x=1$.

2. LinearSolve is not a iterative solver but a symbolic solver for linear algebraic equation system.

Also, it's actually enough to reproduce the issue by discretizing in $x$ direction only, and the resulting ordinary differential equation (ODE) system can be solved symbolically. The following is the corresponding Mathematica code. I've used pdetoode for the generation of ODE system:

eq = D[u[t, x], t] == D[u[t, x], x, x];
ic = u[0, x] == x (1 - x);
bc = u[t, 0] == 0;
points = 4;
grid = Array[# &, points, {0, 1}];
ptoofunc = pdetoode[u[t, x], t, grid, 2];
ode = ptoofunc[eq] // Rest;
odeic = ptoofunc[ic] // Rest;
odebc = ptoofunc[bc];
var = u /@ grid;

asol = DSolve[{ode, odeic, odebc}, var[t] // Through, t] // Simplify // First

$$ \begin{array}{l} u(0)(t)=0 \\ u\left(\frac{1}{3}\right)(t)=\frac{1}{18} \left(-18 t+e^{-18 t}+3\right) \\ u\left(\frac{2}{3}\right)(t)=\frac{2}{9}-2 t \\ u(1)(t)=-3 t-\frac{e^{-18 t}}{18}+\frac{1}{18} \\ \end{array} $$

style = Style[#, 16] &;
ParametricPlot3D[
 Flatten@{t, #} & /@ ({grid, asol[[All, -1]]}\[Transpose]) // Evaluate, {t, 0, 1}, 
 BoxRatios -> {1, 1, 0.4}, AxesLabel -> style /@ {"t", "x", "u"}]

I only used 4 grid points, or else the symbolic calculation will be too slow, but as you can see the result agrees well with the FDM solution.

---

For completeness, here are implementations in some other programming languages. Notice I've only discretized $x$ in the following.

Maple

points := 24: y[0](t):=0: h := 1/points:  
eq := seq((D(y[i]))(t) = (y[i-1](t)-2*y[i](t)+y[i+1](t))/h^2, i = 1 .. points-1), 
          (D(y[points]))(t) = 
            (-y[points-3](t)+4*y[points-2](t)-5*y[points-1](t)+2*y[points](t))/h^2:  
ic := seq(y[points*x](0) = -x^2+x, x = h .. 1, h):  
var := [seq([t, y[i](t)], i = 1 .. points)]:  
p := dsolve([eq, ic], numeric, range = 0 .. 1):  
with(plots):  
odeplot(p, var, size = [default, .618]);

Notice Maple can also calculate the symbolic solution of the sytem with the following line (remember to choose a small points first):

dsolve([eq, ic]);

The solution is consistent with that of Mathematica.

Octave

The code isn't tested in, but should be compatible with MATLAB.

Notice $u(t,0)=0$ isn't plotted in the resulting graph.

points=24;
dx=1/points;
one=ones(points,1);
mat1=spdiags([one -2*one one],[-1 0 1],points-1,points);
mat2=zeros(1,points);
mat2(points-3:points)=[-1 4 -5 2];
mat=[mat1;mat2]/dx^2;
span=dx:dx:1;
[T,Y]=ode45(@(t,y) mat*y,[0,1],span.*(1-span));
plot(T,Y)

Python

from numpy import zeros,linspace
from scipy.integrate import odeint

points=24;h=1/points;

def func(y,t0):
    dydt=zeros(points+1)
    dydt[0]=0
    for i in range(1,points):
        dydt[i]=(y[i-1]-2*y[i]+y[i+1])/h**2
    dydt[points]=(-y[points-3]+4*y[points-2]-5*y[points-1]+2*y[points])/h**2
    return dydt

y0=[x*h-x*x*h*h for x in range(0,points+1)]
t=linspace(0,1,200)
sol=odeint(func,y0,t)

from matplotlib.pyplot import plot,show

plot(t,sol)
show()

As you can see, all these implementations lead to the same solution.

---

Background Information

OK, let me explain a bit about why I insist on finding the meaning of this strange difference scheme. Actually what I've reproduced in this question is the behavior of Mathematica function NDSolve when insufficient b.c. is added to it for solving time dependent PDE. For more information, check this post. As I wrote there:

I know the best countermeasure is to add the missing boundary condition, I'm just curious. Anyway, if the output with insufficient boundary condition is completely meaningless, why doesn't NDSolve simply stop calculating and return the input?

So, perhaps we'll have to admit the solution is really meaningless in the end, but I think this should be the last thing to do.

Answer 1

So the main issue with this is that you are trying to take uniqueness in your approximation, which is based on solving a difference equation, to infer something about your ill-posed PDE problem, like perhaps some unknown physics (i.e. unknown boundary condition). Unfortunately, the approximation cannot reveal such information because it does not imply anything is hidden in the first place.

Given you discretize your PDE into the form $u_{k+1} = A u_k$ where $u_k$ is the solution at time $t_k$, you can start to treat it as a difference equation and see what you can understand.

We know the solution to this difference equation is $u_k = A^k u_0$ and it is trivial to use this result to show uniqueness in the solution for this difference equation. This uniqueness holds basically no matter what form $A$ takes, so you could have some matrix $A$ that represents an ill-posed physics problem and it would still produce a unique solution. The fact it produces a solution does not mean the problem is not still ill-posed or that there's hidden information you can extract, it just means you lost some of that information as you tried to approximate the problem and that loss in information led the approximation to be unique.

Answer 2

This is not the ideal answer; better rigorous analysis can be done here. However, I hope that this has some value.

We all agree that partial differential equation applies to all points in the domain without boundary. Here we apply PDE as well on some points on boundary, and it has consequences.

Indeed, using one side formula on the grid 4x3 (space x time) gives close analytical solution, starting numeration of nodes from zero, numerating by rows:

\begin{equation} u_0:=0\\ u_1:=h(h-1)\\ u_2:=(2h)(2h-1)\\ u_3:=(3h)(3h-1)\\ u_4:=0\\ u_8:=0\\ u_5={\frac {{h}^{2} \left( {h}^{2}+4\,h+6 \right) }{{h}^{2}+3\,h+4}}\\ u_6=4\,{h}^{2}\\ u_7={\frac {{h}^{2} \left( 9\,{h}^{2}+26\,h+34 \right) }{{h}^{2}+3\,h+4}}\\ u_9={\frac { \left( {h}^{3}+6\,{h}^{2}+9\,h+4 \right) h}{{h}^{2}+3\,h+4}}\\ u_{10}=2\, \left( 2\,h+1 \right) h\\ u_{11}={\frac {h \left( 9\,{h}^{3}+28\,{h}^{2}+43\,h+12 \right) }{{h}^{2}+3\, h+4}} \end{equation}

Since we apply one side formula to enforce PDE on right side, for example at node 7, we can equivalently write, \begin{equation} (u_6-2u_7+u^*_6)/h^2 = (u_3-u_{11})/(2h), \end{equation} and solve for $u^*_6$ to satisfy PDE equation using this central finite diffrence scheme, \begin{equation} u^*_6=4\,{\frac {{h}^{2} \left( 3\,{h}^{2}+8\,h+10 \right) }{{h}^{2}+3\,h+4} } \end{equation} with that at hand slope can be evaluated, \begin{equation} u'_7 \approx -2\,{\frac {h \left( 2\,{h}^{2}+5\,h+6 \right) }{{h}^{2}+3\,h+4}} \end{equation} And that what you implicitly apply at that node for this particular case.

In other words, applying on side rule on top-right boundary you assume that PDE governs values at that nodes, so this PDE has to be satisfied as well when central finite difference scheme is applied. That enables you to calculate derivatives implicitly at that nodes. Values of those derivatives are hidden boundary conditions enforced at that nodes. Here PDE itself become boundary condition.

An interesting question to ask, is a wave is reflected from such boundary?

Answer 3

Following the same procedure that for the first problem below, one would have the operator you proposed: $$ \tilde{D}=\frac{1}{\Delta x^2}\begin{bmatrix}-2& 1 & 0 & 0 & \dots & 0 & 0 & 0\\ 1 & -2 & 1 & 0 & \ddots& 0 & 0 & 0\\ \vdots & \ddots & \ddots & \ddots & \ddots & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 1 & -2 & 1 & 0\\ 0 & 0 & 0 & 0 & 0 & 1 & -2 & 1\\ 0 & 0 & 0 & 0 & -1 & 4 & -5 & 2\end{bmatrix} \tag{2}$$

Notice that your matrix $\tilde{D}$ is now applied to the vector $\vec{u} = (u_1,...,u_N)^T$. Note that $\tilde{D}$ has the BC imposed on the first but on the last row there is not any Bc as it should be. Note also that the last row: $r_N$ can be obtained with the following operation: $$r_N = 2r_{N-1} -r_{N-2}$$

Therefore $\tilde{D}$ is non invertible and your problem does not have any definite solution.

This was for the first problem

Let us begin defining the following problem (without any BC as you do): $$ D\,u(x) = f(x) \qquad x\in[0,L] \tag{P}$$ where $D$ is the operator derivative: $$D = \frac{\partial }{\partial x} \tag{1}$$

that acts on some continuous function $u(x)$ in order to provide the derivative: $u_x = D\,u(x)$.

The discrete version of $D$, namely $\tilde{D}$ is the matrix: $$ \tilde{D}=\frac{1}{\Delta x}\begin{bmatrix}-1& 1 & 0 & 0 & \dots\\ 0 & -1 & 1 & 0 & \ddots\\ \vdots & \ddots & \ddots & \ddots & \ddots \end{bmatrix} \tag{2}$$ that acts on the discrete version of the continuous function: $\vec{u}=(u_0,u_1,...,u_N)^{T}$, to provide its discrete derivative: $\vec{u}_x=\tilde{D}\,\vec{u}$

You can see that $\tilde{D}$ is not a square matrix, and the discrete version of $(P)$, i.e. the equation $\tilde{D}\,\vec{u}=\vec{f}$ does not make any sense. This also occurs to $(1)$ where the equation $D\,u=f$ suffers from the same. The problem is that due to the fact that we have more variables than equations, the operator $\tilde{D}$ cannot be inverted and therefore the solution exists up to a constant. It can be shown that if $\vec{u}_0$ a the solution of the problem $(P)$ in its discrete form, any vector $\vec{v}$ such as $\vec{v} =\lambda (1,...,1)^T$ is also a solution. This means that the solution of this problem would be: $$\vec{u} = \vec{u}_0 + \lambda (1,...,1)^T$$ where $\lambda$ is an arbitrary constant. What is more, in the continuous case it is the same!! I mean, this results in continuum reads: $$u(x) = u_0(x) + \lambda$$ where $\lambda$ is again an arbitrary constant.

Now guess where do we set the value of $\lambda$ from!!

To set $\lambda$ only ONE boundary condition is needed to be defined, i.e. $u(0) = 0$ or $u(L) = 0$ for $(1)$ or $u_0=0$ or $u_N=0$ for $(2)$.

This, is theory. The operator $\tilde{D}$, in practice, is defined WITH the BC. For example, if the BC considered is $u_0 = 0$ the operator $\tilde{D}$ reads: $$ \tilde{D}=\frac{1}{\Delta x}\begin{bmatrix} 1 & 0 & 0 & \dots\\ -1 & 1 & 0 & \dots\\ \vdots & \vdots & \vdots & \ddots \end{bmatrix} $$ And it would be applied to the vector of unknowns: $\vec{u} = (u_1,...,u_n)^T$.

This reasoning, can be extended to any order and derivative. Try the operator that approximates the second derivative w.r.t. $x$. You will be convinced of the need for two BC's.

As a result, what you propose: applying the forward and backward discretisation on the first and end nodes respectively a proper BC on one of them: must lead to a non-solvable system.

Answer 4

You are doing explicit time stepping, so solvability is not much of a problem.

You can write for the last grid point $$ du_n/dt = (u_{n-1}-2u_n+u_{n+1})/h^2 $$ with the ghost value $u_{n+1}$ being $$ u_{n+1} = -u_{n-3} + 4 u_{n-2} - 6 u_{n-1} + 4 u_n $$ A Taylor expansion shows $$ u_{n+1} = u_n + h u'_n + (h^2/2) u''_n + (h^3/6) u'''_n + O(h^4) $$ So the ghost value is just obtained by a fourth order extrapolation of the solution.

If the Taylor expansion were of the form $$ u_{n+1} = h u'_n + (h^2/2) u''_n + (h^3/6) u'''_n + O(h^4) $$ then implicitly you would have applied zero Dirichlet bc. If it was of the form $$ u_{n+1} = u_n + (h^2/2) u''_n + (h^3/6) u'''_n + O(h^4) $$ then implicitly you would have applied zero Neumann bc.

But we see all terms upto $O(h^3)$, so no boundary condition is being applied, explicitly or implicitly.

This does not correspond to anything physical or mathematical. The solutions in fact seem to be growing unboundedly with time when I ran the Python code. So I would not attach any significance to the solution you are getting out of this.