How to solve Neumann boundary ODE with shooting method?

Question

$\rho(x)$ is the distribution function.

I'm trying to solve the following equation:

$D_t\frac{\partial^2\rho}{\partial x^2}-\frac{\partial (v_0\rho)}{\partial x}=0$

$D_t\frac{\partial\rho}{\partial x}-v_0\rho=0,\ x=1\ \text{or}\ x=-1$

The equation has two solution

(1) $\rho(x)=0$

(2) $\rho(x)=C_1\exp(v_0(x-1))$

I want to get the (2) solution. I try with the shooting method.

rWall = 1;
v0[x_] := 1;
Subscript[D, t] = 1
With[{rho = rho[x]}, 
 pde = D[rho v0[x] - Subscript[D, t] D[rho, x], x] == 
    NeumannValue[Sign[x] v0[x] rho, Abs[x] == rWall];] ; usolPoisson =
  NDSolveValue[{pde}, rho, {x, -rWall, rWall}, 
  Method -> {"Shooting", 
    "StartingInitialConditions" -> {x[0] == 1, x'[0] == 1}}];

But the solution is (1) $\rho(x)=0$.

How shall I modify it?

Thanks!

Answer 1

First of all, please notice your statement "The equation has two solution" is improper, because the problem has infinite many solutions, and your solution (1) is just a particular case of solution (2) when $C_1=0$.

Then, your trial doesn't have any effect because, the option setting has actually been ignored in this case. Notice you've already mistakenly write x instead of rho! To make it more obvious:

NDSolveValue[{pde}, rho, {x, -rWall, rWall},
 Method -> {"Shooting", 
   "StartingInitialConditions" -> {rho[rWall] == 酋长!, rho'[rWall] == 我们回部落吧!}}]

I think this should be considered as a bug of error checking of NDSolve.

Remark

Inadequate error checking is not rare in NDSolve, if you're in v12.3 or earlier, try the following for fun:

NDSolve[{x'[t] == 1, x[0] == 2}, x, {t, 0, 1}, 
 Method -> {"ExplicitEuler", "兽人!" -> "永不为奴!"}]

This is fixed at least since v13.1, though.

Then why is the method setting be ignored? Because it's not valid. If you read the document carefully, you'll notice the built-in Shooting method solves BVP by searching for the corresponding IVP that satisfies the b.c.s, and IVP solver of NDSolve is utilized during the searching process, but FiniteElement method is not a IVP solver! Related:

How to use Finite elements to solve an initial value ODE with NDSolve?

"OK, then how to circumvent this?" One choice is, if you still want the built-in Shooting method, forget about FiniteElement:

With[{rho = rho[x]},
  bc = rho - D[rho, x] == 0 /. {{x -> -rWall}, {x -> rWall}}];
NDSolveValue[{pde, bc} /. _NeumannValue -> 0, rho, {x, -rWall, rWall}, 
 Method -> {"Shooting", "StartingInitialConditions" -> {rho[0] == 1, rho'[0] == 1}}]
ListPlot@%

NDSolveValue::bvluc

Remark

I've used an undocumented syntax of ListPlot here. See this post for more info.

Another choice is, forget about Shooting, and add a Dirichlet b.c. into NDSolve:

With[{rho = rho[x]}, 
  pde = D[rho v0[x] - Subscript[D, t] D[rho, x], x] == 
    NeumannValue[-v0[x] rho, True]];
NDSolveValue[{pde, rho[rWall] == 3}, rho, {x, -rWall, rWall}]
Plot[%[x], {x, -rWall, rWall}]

Remark

I've removed the Sign[x] because NDSolve doesn't handle it correctly. I believe it's a bug of FiniteElement.

Why can we do this? Because one of the Neumann b.c. is actually redundant. By imposing either of the Neumann b.c., the other Neumann b.c. is automatically satisfied:

With[{rho = rho[x]}, eq = D[rho v0[x] - Subscript[D, t] D[rho, x], x] == 0;
 dbc = rho v0[x] - Subscript[D, t] D[rho, x] == 0 /. {{x -> -rWall}, {x -> rWall}};]

DSolve[{eq, dbc[[1]]}, rho, {x, -rWall, rWall}]
(* {{rho -> Function[{x}, E^x C[1]]}} *)
DSolve[{eq, dbc[[2]]}, rho, {x, -rWall, rWall}]
(* {{rho -> Function[{x}, E^x C[1]]}} *)

So the simplest way to determine a non-trivial solution is, add a b.c. that's not redundant.