NDSolve Extrapolating boundary conditions for FiniteElements

Question

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

My NDSolve call look basically like this:

NDSolveValue[{ 
  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?

Answer 1

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}]

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?

Answer 2

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.