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

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

• What do you do at $i=0$? Are you just setting $f(x_{i-2})$ and $f(x_{x-1})$ to zero? That's a kind of boundary condition as well! – Wolfgang Bangerth Oct 12 '17 at 12:39
• When you approximate a continuous operator by a discrete one, you must always preserve its form. You cannot use other scheme for the borders, this would imply that you are using a boundary condition, whose result is the backward difference formula. The natural discretisation of the first derivative, is a forward/backward scheme, due to the fact that it is unsymmetrical w.r.t. x. – HBR Oct 12 '17 at 14:53
• @hbr You mean, if central difference formula is already used for the discretization of internal area, the one-sided formula can only be used for discretization of boundary condition then? – xzczd Oct 12 '17 at 15:09
• Maybe I have not been very concise. I will extend my comment in an answer to this question. Let me a few minutes. – HBR Oct 12 '17 at 15:12
• @WolfgangBangerth The original example turns out to be improper, I've modified the problem a little, see my edit. – xzczd Oct 12 '17 at 17:12

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.

• Yeah you've fully understanded my question, I think. The main reason that makes me suspect there is a hidden b.c. is, no matter how I build $A$, the final solution will be the same, as long as I've used one-sided formula for the missing boundary. I think it's natural to expect something more interesting. (If the result varied randomly by number of grid points, difference order, etc., I wouldn't have asked this question. ) – xzczd Oct 22 '17 at 4:07
• @xzczd I would expect the numerical solution to be convergent to the same answer because the underlying difference equation dynamics does not really change as you refine or change the discretization schemes. – spektr Oct 22 '17 at 4:12
• Then is it possible to express this underlying dynamics with an equivalent boundary condition? :D – xzczd Oct 23 '17 at 3:16
• @xzczd No. These dynamics are already defined by the matrix $A$ and the initial condition. These dynamics don't depend on $A$ resulting from a well posed physics formulation. There's nothing you will learn from having a unique solution to this approximate solution. – spektr Oct 23 '17 at 3:19

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?

• I took the liberty to add a graphic for the grid to your answer, feel free to roll back if you don't like it :) . – xzczd Oct 16 '17 at 12:28
• Figure is great, as it should be. Note above is just loose attempt to answer your question. – likask Oct 16 '17 at 12:31
• – xzczd Oct 16 '17 at 13:31
• Though still far from the goal, this is the only answer intending to figure out the equivalent boundary condition. You deserve the bounty. :) – xzczd Oct 23 '17 at 3:14

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.

• Well… "As a result, what you propose (applying the forward discretisation on the first node instead a proper BC) leads to a hidden BC. " Then what's this hidden b.c.? Does it have a clear meaning? Can it be expressed in a simpler form? Why can't we use one-sided formula in this way (if using it in this way is just wrong)? These are what I'm interested in. – xzczd Oct 12 '17 at 16:28
• I had a mistake at the end. The system is non-solvable unless you provide a BC, there is no hidden BC. You always will have a non-square matrix. You can verify it if you write down the matrix with the discretisation you proposed (central in the middle and forward/backward at extrema). – HBR Oct 12 '17 at 16:31
• Well, as mentioned in the body of question, this does lead to a solvable system, and the result seems to be quite stable. I've verified it in Mathematica (I can add the source code if you like). Also, notice the first order difference formula seems to be an exception, in that case the matrix is not square, but when the difference order becomes higher, the system becomes solvable. – xzczd Oct 12 '17 at 16:33
• This is a first order hyperbolic equation. The BC is implicitly defined. The solution of this type of equation is well-known to be any function $u$ such as it can be written as a function of the variable $x\pm t$. Therefore the initial condition is just translated from left to the right or vice-versa, depending on the sign of the velocity. The BC is simply the value on the initial time at one of the extrema (depending on the sign of the velocity). You can verify it. But i will repeat myself, every derivative operator must come with a proper BC. – HBR Oct 12 '17 at 16:41
• OK… I admit the previous example chosen in the question isn't proper. (I chose it for simplicity, but forgot the first order hyperbolic equation is actually special. ) Nevertheless, the problem is general, it also appears in other type of PDEs. I've modified the example in the question a little, see my edit. – xzczd Oct 12 '17 at 16:50

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.

• Er… I suggest you to elaborate a bit more on the deduction of $u_{n+1} = -u_{n-3} + 4 u_{n-2} - 6 u_{n-1} + 4 u_n$, though I managed to figure it out, it's not immediately clear, I think. Then, "No boundary condition is being applied", I think it's better to say "No boundary condition is being explicitly applied"? And what I'm interested in is, whether this implicit b.c. can be expressed in a clearer way, any idea? – xzczd Oct 17 '17 at 4:43
• There is no implicit bc being applied. You are doing pure extrapolation to the ghost cell. – cfdlab Oct 17 '17 at 5:34
• Hmm… but as the title asked, why did we find a unique solution if no b.c. is applied on the right side? – xzczd Oct 17 '17 at 6:05
• If you do explicit time stepping, getting a solution is no big deal. You dont know what problem this solution is solving. So it is pointless to go along this path. – cfdlab Oct 17 '17 at 6:25
• Just because you got some answer instead of getting inf, nan or crashing, does not mean there is something real happening here. You should start with a properly posed problem and then construct a consistent and stable scheme for that. Instead you write down a random scheme and ask "Why am I getting something instead of garbage". This is wrong question to ask. – cfdlab Oct 18 '17 at 5:35