4
$\begingroup$

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: Numerical solution of (1), (2)

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:

enter image description here

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.

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

2 Answers 2

3
$\begingroup$

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

enter image description here

Additional remarks

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

enter image description here

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]

enter image description here

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.

$\endgroup$
5
  • $\begingroup$ 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? $\endgroup$ Nov 25, 2015 at 12:17
  • 1
    $\begingroup$ @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. $\endgroup$ Nov 25, 2015 at 15:55
  • $\begingroup$ 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). $$ $\endgroup$ Nov 25, 2015 at 20:01
  • 1
    $\begingroup$ @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.) $\endgroup$ Nov 25, 2015 at 20:21
  • $\begingroup$ $\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)$. $\endgroup$ Dec 3, 2015 at 19:42
3
$\begingroup$

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

enter image description here

gives the desired answer.

$\endgroup$
3
  • $\begingroup$ 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? $\endgroup$ Nov 24, 2015 at 6:34
  • $\begingroup$ 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? $\endgroup$ Nov 24, 2015 at 6:54
  • 1
    $\begingroup$ @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. $\endgroup$
    – bbgodfrey
    Nov 24, 2015 at 14:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.