I'm trying to solve a 1+1 (1d time + 1d space) partial differential diffusion equation.

My NDSolve call look basically like this:

  D[u[x, t], t] == (* RHS *), 0.1, 0.1, 2], 
  u[x, 0] == (* something *), u[0, t] == 0, (* ? *)}, 
  u[x, t], {x, 0, xmax}, {t, 0, tmax}, 
  Method -> {"MethodOfLines", "TemporalVariable" -> t, 
   "SpatialDiscretization" -> {"FiniteElement"}}

I have an initial and a left boundary condition for u. The problem I'm facing is, that I want to have linear extrapolated right boundary conditions.

I know, how this is done, when I discretize the spatial domain myself, I just do

$u_{i+1} = 2u_i - u_{i-1}$,

where $i$ is the spatial index.

Maybe another question I wound like to get answered is, if the mesh for the spatial domain gets adapted during the time integration?

  • $\begingroup$ @andre314 I'm trying to do a linear extrapolation $\endgroup$
    – Bomel
    Dec 22, 2022 at 21:25
  • 2
    $\begingroup$ Your question is not very clear, at least to me. For example what is a 'linear extrapolated right boundary condition' ? You should provide an actual example in stead of having us guess what you mean. $\endgroup$
    – user21
    Dec 24, 2022 at 5:42
  • $\begingroup$ It appears that you wish to have u''[x] vanish at the x = xmax boundary. This is not permitted, if the highest derivative in x in the PDE also is second order. However, you could try applying u''[x]->0 to the PDE at x = xmax, which might give a simpler boundary condition, perhaps that u[x] is constant in time there. $\endgroup$
    – bbgodfrey
    Dec 28, 2022 at 15:11

2 Answers 2


If I understand what you mean by "linear extrapolated right boundary conditions" correctly: you don't need such stuff for NDSolve, no matter what spatial discretization method you've chosen. (Currently there're only TensorProductGrid and FiniteElement, of course. )

Let me explain with a toy example:

$$u'(x)=\sin(x),\ u(0)=0,\ x\in[0,2\pi]$$

Suppose we want to solve it via finite difference method, and decide to discretize the equation based on the central difference formula

$$u'(x)\approx\frac{u(x+\Delta x)-u(x-\Delta x)}{2\Delta x}$$

This can be easily done in Mathematica as follows:

{eq, ic} = {D[u[x], x] == Sin[x], u[0] == 0};

domain = {xL, xR} = {0, 2 Pi};
points = 25;
dx = (xR - xL)/(points - 1);
grid = Array[# &, points, domain];

deq = Table[(u[x + dx] - u[x - dx])/(2 dx) == Sin[x], {x, grid}];

The last step is to solve deq together with ic… wait!:

Cases[{deq, ic}, u[_], Infinity] // Union
(* {u[0], u[-(π/12)], u[π/12], u[π/6], u[π/4], u[π/3], 
 u[(5 π)/12], u[π/2], u[(7 π)/12], u[(2 π)/3], u[(3 π)/4], 
 u[(5 π)/6], u[(11 π)/12], u[π], u[(13 π)/12], u[(7 π)/6], 
 u[(5 π)/4], u[(4 π)/3], u[(17 π)/12], u[(3 π)/2], u[(19 π)/12], 
 u[(5 π)/3], u[(7 π)/4], u[(11 π)/6], u[(23 π)/12], u[2 π], 
 u[(25 π)/12]} *)
% // Length
(* 27 *)
{deq, ic} // Flatten // Length
(* 26 *)

There're 27 different u[_] terms in the system, but we only have 26 equations at hand, the system is not closed!

That's the reason why the "linear extrapolated right boundary conditions" (here I simply follow your terminology) is utilized. Based on the truth

$$\lim_{x\rightarrow 2\pi^+}u'(x)=\lim_{x\rightarrow 2\pi^-}u'(x)$$

and one-sided first order difference formula

$$u'(x)\approx\frac{u(x+\Delta x)-u(x)}{\Delta x}$$

we deduce

$$u(x+\Delta x)=2u(x)-u(x-\Delta x)$$

deqextra = u[x] == (u[x - dx] + u[x + dx])/2 /. x -> xR;

assuming the grid is uniform. Now the system is closed and can be solved without difficulty:

sollst = NSolveValues[{deq, ic, deqextra} // Flatten, 
   Table[u[x], {x, xL - dx, xR + dx, dx}]][[1]]

It's clear the "linear extrapolated right boundary conditions" is not a boundary condition in usual sense. It's merely a technique for supplying the equation system when discretizing differential equations with difference formula whose order is higher than the differential order of the equation. This is not needed for FiniteElement method:

solfem = NDSolveValue[{eq, ic}, u, {x, xL, xR}, Method -> "FiniteElement"]

You'll see NDSolveValue::femcscd warning, but it's not a problem in this case:

ListPlot[sollst[[2 ;; -2]], DataRange -> domain]~Show~Plot[solfem[x], {x, xL, xR}]

enter image description here

Elaborating on how TensorProductGrid and FiniteElement closes the system is too complicated not quite related to your problem, so I'd like to stop here.

As to the question

Is the mesh for the spatial domain gets adapted during the time integration?

It is already answered here:

Is it possible to solve a differential equation with a user-defined variable mesh in NDSolve?

  • $\begingroup$ But as I understand your answer, you look at a 1d IVP (initial value problem). Which I can integrate from x=xL to x=xR. Your ode doesn’t need a second boundary condition, since you eq and initial condition determines your solution. For a 1+1 PDE, I need to impose a boundary to calculate the next time slice. $\endgroup$
    – Bomel
    Dec 23, 2022 at 7:25
  • $\begingroup$ @Bomel …The issue is exactly the same. Notice I've intentionally used central difference formula to solve the problem as a BVP to illustrate the problem. Think carefully about what the discretized PDE is. $\endgroup$
    – xzczd
    Dec 23, 2022 at 7:42
  • $\begingroup$ @xzcxd Okay, I think I understand, what you mean. But what is Mathematica doing then at the undefined boundaries? Since the solution of Mathematica clearly doesn't match with the reference solution I have, when I discretize the spatial domain by hand and solve the coupled ode. $\endgroup$
    – Bomel
    Dec 23, 2022 at 8:37
  • $\begingroup$ @Bomel That should not happen as long as the reference solution itself is correct and NDSolve is adjusted properly. Please show us your specific problem. Also, please notice that, if by "diffusion equation" you mean something similar to $u_t=u_{xx}$, then you cannot close this type of system with "linear extrapolated right boundary condition" because it's not a valid approximation for boundary at infinity, and it doesn't even have a clear physical meaning in this case. Please read the following post carefully (don't miss the link there): mathematica.stackexchange.com/q/73961/1871 $\endgroup$
    – xzczd
    Dec 23, 2022 at 8:56
  • $\begingroup$ I'm not trying to solve a pde with a nonlinear diffusion flux $u_t = (F(u_x, t))_x$ and u becomes "non-analytic" over time (this is why I want to use the Inactive feature of the $F(u_x, t)$ derivative). I try to create an example problem, since I cannot post the original problem because I'm not so sure about the copyright stuff. Furthermore, I'm not trying to model boundaries at Infinity at $x=x_{max}$, mathematically $u(x_{max},t)$ has to be linear. $\endgroup$
    – Bomel
    Dec 23, 2022 at 10:36

Automatic temporal adaptivity:yes ( adaptive stepsize control), spatial adaptivity: no. For extrapolation see the documentation for NDSolve FEM options. You are looking for the 'ExtrapolationHandler', This is an option for 'ElementMeshInterpolation'. If you need something else then you need to re-create an Interpolation function.

  • $\begingroup$ I don't quiet get, how I create an extrapolating boundary condition for my PDE using options for ElementMeshInterpolation? $\endgroup$
    – Bomel
    Dec 22, 2022 at 20:48
  • $\begingroup$ @Bomel, see update. $\endgroup$
    – user21
    Dec 22, 2022 at 21:06
  • $\begingroup$ But as I understand this, it only defines the Extrapolation behaviour for the ElementMeshInterpolation object. If I could configure the ExtrapolationHandler in a way to return a linear extrapolation. Will NDSolve make use of that and extrapolate at xmax? $\endgroup$
    – Bomel
    Dec 22, 2022 at 21:17
  • $\begingroup$ If you want linear extrapolation use InterpolationOrder->1 $\endgroup$
    – user21
    Dec 24, 2022 at 6:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.