Comparison between numerical solution of nonlinear ode and nonlinear ode of second order

Question

Background

Let's consider the following initial value problem for nonlinear system $$ \begin{cases} E' &=& 1 - n_e, \\ n_e' &=& -8\,n_e\,E, \end{cases} \tag{1} $$ with the following initial conditions $$ n_e(0) = 0.999,\quad E(0) = 4\cdot10^{-4}. \tag{2} $$ WM code is as follows:

L = 6;
a = -8;
{ne, Ef} = NDSolveValue[{Ef'[x] == 1 - ne[x],ne'[x] == a ne[x] Ef[x], 
Ef[0] == 4 10^-4, ne[0] == 0.9999},
{ne, Ef}, {x, 0, L}, Method -> "StiffnessSwitching"];
Plot[{ne[x], Ef[x]}, {x, 0, L}, PlotRange -> All, ImageSize -> 400]

Plot of numerical solution of (1), (2) to be the following:

Also, one could express $E$ from (2): $$ E = -\frac{n'_e}{8n_e} \tag{3} $$ and substitute (3) in (1): $$ -\biggl(\frac{n'_e}{8n_e}\biggr)' = 1 - n_e. \tag{4} $$ Then, nonlinear equation is arised: $$ n''_e - \dfrac{(n'_e)^2}{n_e} = -8(1 - n_e)n_e. \tag{5} $$ Now, boundary conditions are the following $$ n_e(0) = 0.999,\quad n_e(6) = 0. \tag{6} $$ WM code:

sol01 = NDSolveValue[{ne''[x] - 1/ne[x] (ne'[x])^2 == -8 (1 - ne[x]) ne[x], 
ne[0] == 0.999, ne[6] == 0}, ne, {x, 0, 6}]

Output is as follows:

Question

So, my question wether it is possible to find similar numerical solution of (5),(6) to numerical solution of (1),(2) or not?

Notice

One could also find the following substitution $$ \varphi = \ln{n_e},~n_e = \exp{\varphi}. \tag{7} $$ So, (5) transforms into (8) $$ \varphi'' = -8(1 - \exp{\varphi}). \tag{8} $$ The question expand on (8) as well.

UPD1. Misprints in (5) fixed.

Answer 1

We solve here the 3D equation

$$\nabla ^2\phi =-8 (1-\exp (\phi ))$$

in the special case of spherical symmetry, and compare the numerical solution with the exact soluton of the linearized problem. We leave a possible derivation of this equation in the 3D case to the author of the OP.

The equation to be solved is in Mathematica notation

eq3 = u''[r] + 2/r u'[r] == 8 (Exp[u[r]] - 1);

To simplify the typesetting in Mathematica we formulate the equation in terms of u[r] instead of \[Phi][r].

In analogy to the 1D case we shall also consider the quantity

n[r_] = Exp[u[r]];

Numerical solution

The numerical solution must be accompanied by two boundary conditions. It seems natural to let u'[0] == 0 and u[0] == u30, a small negative number taken from the 1D case. In the numerical treatment we have to replace 0 by some small number, say r0 = 0.001.

u30 = Log[0.999]

(* Out[86]= -0.0010005 *)

uu[r_] = u[r] /. 
    NDSolve[{eq3, u[0.001] == u30, u'[0.001] == 0}, u[r], {r, 0, 8}] // Quiet;

n3 = Exp[uu[r]];

Exact solution of the linearized equation

The linearized equation is

eq3L = u3L''[r] + 2/r u3L'[r] == 8 u3L[r];

Its solution is

uu3L = u3L[r] /. DSolve[eq3L && u3L[0] == u3L0, u3L[r], r][[1]]

(*
Out[92]= (E^(-2 Sqrt[2] r) (-1 + E^(4 Sqrt[2] r)) u3L0)/(4 Sqrt[2] r)
*)

It is interesting that the finiteness at r = 0 alone already takes care for the condition u'[0]==0.

Indeed

Series[((E^(-2 Sqrt[2] r)) (-1 + E^(4 Sqrt[2] r)) )/(4 Sqrt[2] r), {r, 0, 2}]//Normal

(*
Out[108]= 1 + (4 r^2)/3
*)

Now, for comparison, we let

u3L0 = u30;

n3L = Exp[uu3L];

and plot both solutions together

Plot[{n3, n3L}, {r, 0, 6}, AxesLabel -> {"r", "n3, n3L"}, 
 PlotLabel -> 
  "Solution of 3D eqution in spherical symmetry\nnumerical solution -> n = \
blue\nlinear approximation -> n = yellow"]

Additional remarks

The same comparison can be made for the 1D case with the result

In the 1D case you get a good overview over the solution using StreamPlot[]:

StreamPlot[{- 8 fE n, 1 - n}, {n, -1, 3}, {fE, -1, 1}, 
 PlotLabel -> "ODEs E' = 1 - n, n' = - 8 E n\nStreamplot of dE/dn", 
 AxesLabel -> {"n", "E"}, Frame -> False, Axes -> True]

The explicit for of the streamlines (orbits) are easily obtained exactly by solving 8 dE/dn = (1-n)/(E n) to give

fE = {Sqrt[c + n - Log[n]], -Sqrt[c + n - Log[n]]};

with an integration constant c.

The separatrix is given by c = -1.

Answer 2

Equation (4) should read ne''[x] - ne'[x]^2/ne[x] == -8 (1 - ne[x]) ne[x]. Additionally, boundary conditions derived from the original equations are ne[0] == 0.999, ne'[0] == -32*^-4 ne[0] With these changes,

sol01 = NDSolveValue[{ne''[x] - ne'[x]^2/ne[x] == -8 (1 - ne[x]) ne[x], ne[0] == 0.999, 
    ne'[0] == -32*^-4 ne[0]}, ne, {x, 0, 6}, Method -> "StiffnessSwitching"];
Plot[sol01[x], {x, 0, 6}, AxesLabel -> {x, ne}]

gives the desired answer.