Skip to main content
deleted 2 characters in body
Source Link
xzczd
  • 68.4k
  • 9
  • 174
  • 489

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

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

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

added 1 character in body
Source Link
xzczd
  • 68.4k
  • 9
  • 174
  • 489

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

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

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

Source Link
xzczd
  • 68.4k
  • 9
  • 174
  • 489

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 I 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[_] term 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?