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

### 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.

• Well, your (3) is valid iff ne != 0 ... Nov 23, 2015 at 19:54
• Fore sure. That is a question how to find a suitable analog... Nov 23, 2015 at 19:56

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


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.

• Thank you so much, Dr. Wolfgang! Let me also ask about nonsymmetrical case. Suppose, that we consider 3D cartesian $\varphi=\varphi(x,y,z)$ $$\nabla^2 \equiv \partial_x + \partial_y + \partial_z. \tag{\ast}$$ How should we solve BVP (or IVP) of $$\nabla^2 \varphi = -8 (1 - \exp(\varphi))$$ to obtain ''hat'' like solution? Nov 25, 2015 at 12:17
• @olekravchenko: some remarks: 1) what are the boundary conditions for your partial differental equation in (x,y,z)? 2) I suggest to consider the linear approximation first, as it is a classical case, where all solution methods are well known. 3) In the extremely nonlinear case where you drop the 1, for 1D you get the solution Log[E^A*Sec[(Sqrt[E^A]*x)/Sqrt[2]]^2];which shows a periodic divergence in a the points xn = pi Sqrt(2) Exp(-A) (n+1/2), n integer. Nov 25, 2015 at 15:55
• 1) well, that is the question I've asked in the very beginning, even in 1D case I still have no idea about coressponding BVP with numerical solution like in IVP (1), (2). For example, in 2D I expect somewhat like "hat" function which one could obtain by rotation of $n_3(x^2+y^2)$ solution. In 3D I also expect somewhat like a "semi-hat" solution. So, BC in 2D to be like that: $$\left.\varphi(x,y)\right|_{x=0} = n_3(y),\, \left.\varphi(x,y)\right|_{y=0} = n_3(x),\, \left.\varphi(x,y)\right|_{x=L} = n_3(L),\, \left.\varphi(x,y)\right|_{y=L} = n_3(L).$$ Nov 25, 2015 at 20:01
• @olekravchenko: the answer to your preliminary question is "no", as belisarius already said. You can't have n(0) = n0 with n !=0 and n(x0) = 0 for any finite x0. It is exactly the region of small phi for which the linearization is valid, and you can study the approach to 0. I don't understand what you mean by "hat" function. I have shown how to solve the 3D symmetric case. You can do the 2D symmetric case similarly by yourself. Also, as I said, in the linear approximation you can easily get the results for any reasonable boundary conditions using the standard procedures (separation of vars.) Nov 25, 2015 at 20:21
• $\text{@Dr.WolfgangHintze}:$ thank you for your answer! I mean "hat", cause the considered are is a quoter of a square with center point in $(0;0)$. Dec 3, 2015 at 19:42

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


• Thank's a lot for the fruitful comment. Definitely yes, it's an answer. But, if fact we solve IVP for $n_e$, don't we? What about corresponding BVP? Is it possible to find suitable BVP which gives the similar numerical solution as IVP, isn't it? Nov 24, 2015 at 6:34
• What about multidimensional case? For example, let's consider $$\nabla^2 \varphi = -8 (1 - \exp(\varphi)). \tag{\ast}$$ Should we consider IVP problem for $(\ast)$ instead of BVP, shouldn't we? Nov 24, 2015 at 6:54
• @olekravchenko To solve the 2nd order ODE requires an appropriate boundary condition at x = 6. ne[6] == 0 is not the right choice, because it leads to ne = 0 everywhere, which is inconsistent with ne[0] == .999. You should pose your PDE as a separate question, because it is well beyond the scope of the present question. Nov 24, 2015 at 14:02