6

So I have a set of differential equations (which I know are analytically integrable) but I want to integrate them numerically to get the form of the solution (the solution is not explicit and hence I want to numerically integrate it).

g=7/5

v = Function[r, 1/(n[r] r^2)]

eulerEq = v[r]*v'[r] + n'[r]*n[r]^(g - 2) + 2/r^2 == 0

initcond = n[1] == 1 (*These are my equation and initial condition *)

Now to numerically solve it.

S = NDSolve[{eulerEq, initcond}, n[r], {r, 0.5, 10}]

When I hit shift+Enter, it shows me the following error:

Power::infy: Infinite expression 1/0. encountered. >>

Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered. >>

NDSolve::ndnum: Encountered non-numerical value for a derivative at r == 1.`. >>

I tried using the debugger and I realized that it is computing $n'[1]$ in the process and it is getting a 'division-by-zero'. I can integrate the equation easily and get the constant of integration using the initial condition (it is getting $n$ explicitly in terms of $r$ that is troublesome). Is there a way to do this still using NDSolve or is there an alternative? I want to do it numerically only because this equation is just a special case which is analytically integrable and I need to get the other solutions also.

bbgodfrey
  • 61,439
  • 17
  • 89
  • 156
Sunil S.
  • 135
  • 7

3 Answers3

9

Determining n'[r] from

Solve[eulerEq, Derivative[1][n][r]][[1, 1, 2]]
(* -((2 n[r] (-1 + r^3 n[r]^2))/(r (-1 + r^4 n[r]^(12/5)))) *)

we see that it is indeterminate at r = 1, and it is for this reason that NDSolve fails there. Nonetheless, the limiting value can be obtained as follows,

Limit[% /. n[r] -> 1 + a (r - 1), r -> 1]
(* -((5 (3 + 2 a))/(10 + 6 a)) *)
sl = Solve[% == a]
(* {{a -> 1/6 (-10 - Sqrt[10])}, {a -> 1/6 (-10 + Sqrt[10])}} *)

and used to define boundary conditions very near r = 1.

dr = .0001;
s1 = NDSolveValue[{eulerEq, n[1 + dr] == 1 + dr a /. First[sl]}, n, {r, 1 + dr, 10}];
s2 = NDSolveValue[{eulerEq, n[1 + dr] == 1 + dr a /. Last[sl]}, n, {r, 1 + dr, 10}];
s3 = NDSolveValue[{eulerEq, n[1 - dr] == 1 - dr a /. First[sl]}, n, {r, .5, 1 - dr}];
s4 = NDSolveValue[{eulerEq, n[1 - dr] == 1 - dr a /. Last[sl]}, n, {r, .5, 1 - dr}];

p1 = Plot[{s1[r], s2[r]}, {r, 1 + dr, 10}, PlotRange -> {{0, 10}, {0, 5}}];
p3 = Plot[{s3[r], s4[r]}, {r, .5, 1 - dr}, PlotRange -> {{0, 10}, {0, 5}}];
Show[p1, p3, AxesLabel -> {r, n}]

enter image description here

bbgodfrey
  • 61,439
  • 17
  • 89
  • 156
  • good solution you managed to get it. Maple does actually solve this analytically for same initial condition, and I get the same plot as you show. The solution n(r) it gives is in terms of a complicated non-linear function in n(r) – Nasser Dec 18 '15 at 03:51
  • @Nasser so Maple uses a different algorithm to compute the numerical integral? Is it possible to implement it here? Because it would be tedious to compute the limit every time such a situation pops up. – Sunil S. Dec 18 '15 at 04:17
  • I think you can explain a bit more about why you choose n[r] -> 1 + a (r - 1) but not n[r] -> 1 + a (r - 1)^2, n[r] -> 1 + a (r - 1)^3, etc. :) – xzczd Dec 18 '15 at 04:23
  • @SunilS. I do not know how Maple solved it, but here is the solution it gives. There is no numerical anything. it solved the analytical ode. Not numerical. The plots seems similar to what in the above answer: Mathematica graphics here is plain text of the solution of the ODE itself. 1+5*n(r)^(12/5)*r^4+(-2*r^4-4*r^3)*n(r)^2 = 0 so it is non-linear in n(r) – Nasser Dec 18 '15 at 04:32
  • @Nasser I see. @bbgodfrey Solve[eulerEq, Derivative[1][n][r]][[1, 1, 2]] seems to give me the error Part::partw: "Part 1 of {} does not exist." edit: Oops, my bad! Sorry about that. Used wrong variable names – Sunil S. Dec 18 '15 at 04:40
  • 1
    @xzczd Linear was my first guess, and it worked. Once assumed, it is easy to show that the guess is self-consistent, as is apparent from the plot.. – bbgodfrey Dec 18 '15 at 05:06
  • Can you explain what First[s1] and Last[s1] do here? I'm unable to understand their function. Thanks a lot guys – Sunil S. Dec 18 '15 at 05:06
  • @SunilS. sl has two values, and First and Last pick out the first and second values, respectively. – bbgodfrey Dec 18 '15 at 05:08
  • @bbgodfrey. But s1 seems to have a lot more than two values. s1[[3]] and s1[[4]] return some tables. Also, can you explain how it is recursively using First[s1] to find s1 – Sunil S. Dec 18 '15 at 05:12
  • @SunilS. No. sl has just two values, as shown in my answer above. {a -> 1/6 (-10 - Sqrt[10])} and {a -> 1/6 (-10 + Sqrt[10])}. – bbgodfrey Dec 18 '15 at 05:17
  • Let us continue this discussion in chat. – Sunil S. Dec 18 '15 at 05:19
  • @SunilS. Too late at night. Sorry. – bbgodfrey Dec 18 '15 at 05:19
7

bbgodfrey has already given a good answer, here's just a not that rigorous but cheap and accurate enough way to resolve the problem:

S = NDSolveValue[{eulerEq, n@1 == 1 + #}, n[r], {r, 0.5, 10}, 
    Method -> "StiffnessSwitching"] & /@ ({1, -1} $MachineEpsilon)

Plot[S, {r, 0.5, 10}, PlotRange -> All]

enter image description here

Community
  • 1
xzczd
  • 65,995
  • 9
  • 163
  • 468
  • Why does {1.7, -1.7} work as the perturbation, when {1.8, -1.8} does not? Perhaps, the boundary of the basins of attraction. In any case, +1. – bbgodfrey Dec 18 '15 at 05:02
  • @bbgodfrey Interesting observation. I guess it might be related to how "StiffnessSwitching" works (try Manipulate[With[{S = NDSolveValue[{eulerEq, n@1 == 1 + #}, n[r], {r, 0.5, 10} , Method -> "StiffnessSwitching"] & /@ ({1.7, -a} $MachineEpsilon)}, Plot[S, {r, 0.5, 10}, PlotRange -> {0, 5}, Epilog -> {PointSize@Medium, Point[{1, 1}]}]], {a, 1, 10}]), but currently I can't find a better Method to handle the singularity. – xzczd Dec 18 '15 at 05:58
2

Using the limiting values from bbgodfrey's answer, we can take a more direct approach with NDSolve using "EquationSimplification" -> "Residual":

sol = Join @@@
       (NDSolve[{eulerEq, initcond, n'[1] == #}, 
         n[r], {r, 0.5, 10},
         Method -> {"EquationSimplification" -> "Residual"}] & /@
          {1/6 (-10 + Sqrt[10]), 1/6 (-10 - Sqrt[10])});

Plot[n[r] /. sol // Evaluate, {r, 0.5, 10}, PlotRange -> All, 
 AxesOrigin -> {0, 0}]

Mathematica graphics

Community
  • 1
Michael E2
  • 235,386
  • 17
  • 334
  • 747