4
$\begingroup$

I need to solve a ODE for a function $f_\sigma(x)$, where $\sigma$ is a free parameter to tune. All the solutions are supposed to be singular for some value of $x$, except one for a given value $\sigma=\sigma^*$. $f_\sigma(x)$ behaves nicely for every $x$, and as I get close to $\sigma^*$, the value at which $f(x)$ diverges moves further away (i.e. to larger $x$).

I need to find $\sigma^*$. My idea would be to do a iterative procedure in which I compute $f_\sigma(x)$ for different values $\sigma$, and then focus on the range of values for which the singularity happens the latest (the largest $x$). However I do not know how to extract this information from NDSolve.

kd = (2^(d - 1) π^(d/2) Gamma[d/2]);
d = 3;
ext = 100;
s2[σ_] := 
  NDSolve[{f''[p] -  2 f[p] f'[p] + (1 + d/2) f[p] + (1 - d/2) p f'[p] == 0,  
           f'[0] == σ, 
           f[0] == 0
          }, f[p], {p, 0, ext}
  ]

Doing for example s2[-.3], Mathematica tells me

NDSolve::ndsz: At p == 4.160343901296758`, step size is effectively zero; singularity or stiff system suspected. >>

From this I would like to extract the value p=4.1603....., but I do not know how to do this.

Any tips?

$\endgroup$
4
$\begingroup$

Although the overall conclusions of this answer are unchanged, details have been edited substantially.

I surmise from the question and comments that the OP is seeking a solution to the ODE that does not grow exponentially at large p. Two solutions already identified are f[p] = 0 and f[p] = p. In addition to these trivial solutions, another exists near σ = -0.228601, and other answers have focused on obtaining more accurate values for this σ. Unfortunately, even with highly accurate values for σ, f[p] has been computed only for modest values of p. This answer takes a different approach, finding the desired f[p] directly and determining the corresponding σ only as a byproduct.

If f[p] is not to explode in magnitude for large p, f''[p] must be small compared to other terms in the ODE. So, let us drop this term from the ODE to obtain an asymptotic equation, which can be solved for f'[p].

Solve[- 2 f[p] f'[p] + (1 + d/2) f[p] + (1 - d/2) p f'[p] == 0, f'[p]][[1, 1]]
(* f'[p] -> (5 f[p])/(p + 4 f[p]) *)

As shall be seen, f[p] is approximately equal to p for large p, so, f'[p] approaches 1 asymptotically. This is sufficient information to obtain the desired solution, computed out to pmax = 10000.

pmax = 10000; 
{a1, a2} = NDSolveValue[{f''[p] - 2 f[p]  f'[p] + (1 + d/2) f[p] + (1 - d/2) p f'[p] == 0, 
    f[0] == 0, - 2 f[pmax]  f'[pmax] + (1 + d/2) f[pmax] + (1 - d/2) pmax f'[pmax] == 0},
    {f, f'}, {p, 0, pmax}, WorkingPrecision -> 30, MaxSteps -> 10^6, 
    Method -> {"Shooting", "StartingInitialConditions" -> {f[pmax] == pmax - 1754/100,
    - 2 f[pmax] f'[pmax] + (1 + d/2) f[pmax] + (1 - d/2) pmax f'[pmax] == 0}}];

Plot[{a1[p], a2[p]}, {p, 0, 10.4}, AxesLabel -> {p, "f, f'"}, PlotRange -> {All, {-1, 7}}]

enter image description here

which agrees with the plot in my earlier answer for as far as the former is valid. However, the current answer is valid to p -> 10000. As an accuracy test (as suggested by Michael E2),

a1[0]
(* -1.43919293846595796615*10^-11 *)
σ = a2[0]
(* -0.22859820245788122180105821189 *)

The accuracy of σ, although good, agrees with Michael E2's value only to ten significant figures. However, that is beside the point. By integrating from large p to small p rather than the reverse, an excellent value for σ is not needed.

To plot the solution over the entire range, {p, 0, 10000}, it is convenient to do so in terms of a new dependent variable.

g[p] = f[p] - p
LogLogPlot[{-b1[p] + p, -b2[p] + 1, b2'[p]}, {p, .01, pmax}, 
    AxesLabel -> {p, "-g, -g', g''"}]

enter image description here

where the blue and orange curves are -g[p] and -g'[p], and the green curve is g''[p]. As expected, f''[p] = g''[p] is very small asymptotically, and g[p] = f[p] - p is much less than p.

As an aside, the asymptotic ODE can be solved to obtain the asymptotic solution. First, substitute g for f,

Simplify[Unevaluated[D[f[p], p, p] - 2 f[p] D[f[p], p] + (1 + d/2) f[p] + 
    (1 - d/2) p D[f[p], p] ] /. f[p] -> p + g[p]]
(* 1/2 (g[p] - 5 p g'[p] - 4 g[p] g'[p] + 2 g''[p] *)

drop the second derivative term, and solve using DSolve.

dsol = DSolveValue[g[p] - 5 p g'[p] - 4 g[p] g'[p] == 0, g[p], p] /. C[1] -> Log[c]
(* Root[-c p - c #1 + #1^5 &, 1] *)

The constant c can be determined by fitting the asymptotic solution to the full solution above.

FindRoot[(dsol /. p -> pmax) == a1[pmax] - pmax, {c, -167}][[1]]
(* c -> -166.625 *)

The same calculation for pmax = 1000 yields c -> -166.457, indicating that the asymptotic state has indeed been reached.

$\endgroup$
  • $\begingroup$ I'm not sure what to think about this method (other than it's clever). The value obtained for σ = a1'[0] does not seem to be very good....Or do you have an argument why it is? $\endgroup$ – Michael E2 Oct 17 '15 at 21:02
  • $\begingroup$ @MichaelE2 For this answer I took as a goal not to obtain a better value for σ but to obtain a solution for f that satisfies a1[0] =0 reasonably well ( -6.8455349099784265381*10^-12 for WorkingPrecision->30) and is valid to p = 100. I shall update my second answer now with this higher precision result. Thanks for the question. $\endgroup$ – bbgodfrey Oct 17 '15 at 21:18
  • $\begingroup$ Fair enough, but it seems clear the question asks how to estimate the critical value σ*. One could push this method to Infinity (I think), but it's not clear one gets a more accurate estimate of the critical value σ* the further one goes. Integrating this vector field (ode) backwards seems stable, but I don't know how to estimate the error. OTOH, it's not clear that shooting as far as one can from p == 0, as in the other answers, is superior either....If I get a chance, I'll think about it. This method might actually reduce error. $\endgroup$ – Michael E2 Oct 17 '15 at 21:32
  • $\begingroup$ The Q says, "I need to find σ*." That's what I was thinking of. Maybe someone added that. In any case, getting it to 11 digits is probably much more than anyone needs in practice. Assuming they're correct. :) $\endgroup$ – Michael E2 Oct 17 '15 at 21:43
  • $\begingroup$ @bbgodfrey thanks, the answer looks great. I will probably need this for the future. A couple of questions. How did you see that asymptotically f[p]=p? and how did you other term in the initial shooting condition, f[pmax] == pmax - 1754/100? $\endgroup$ – bnado Oct 19 '15 at 17:32
2
$\begingroup$

I modified your NDSolve a bit for convenience (NDSolveValue to get rid of the rule, and f instead of f[p] to get a pure function):

s2[σ_] := 
 NDSolveValue[{f''[p] - 2 f[p] f'[p] + (1 + d/2) f[p] + (1 - d/2) p f'[p] == 0, 
   f'[0] == σ, f[0] == 0}, f, {p, 0, ext}]

s2[σ]["Domain"] can be used to inspect the domain. Now examining a range of sigmas:

Table[{σ, s2[σ]["Domain"][[1, 2]] // Quiet}, {σ, -5, 5, .05}] // ListPlot

Mathematica graphics

Looks like there are three candidate positions. Examining their neighborhood:

ListPlot[
  Table[{σ, s2[σ]["Domain"][[1, 2]] //Quiet}, {σ, -0.229, -0.2284, .0000001}], 
  PlotRange -> All]

Mathematica graphics

ListPlot[
  Table[{σ,s2[σ]["Domain"][[1, 2]] // Quiet}, {σ, -0.001, 0.001, .00002}], 
  PlotRange -> All]

Mathematica graphics

ListPlot[
  Table[{σ, s2[σ]["Domain"][[1, 2]] // Quiet}, {σ, 0.99999, 1.00001, .0000001}], 
  PlotRange -> All]

Mathematica graphics

Looks like the one at σ = 0 is the winner. X runs up to its full range at 100.

$\endgroup$
  • $\begingroup$ The problem with σ = 0 is that it yields f identically zero. $\endgroup$ – bbgodfrey Oct 17 '15 at 14:34
  • 1
    $\begingroup$ @bbgodfrey Where in the question did you read that that should not be happening? ;-) $\endgroup$ – Sjoerd C. de Vries Oct 17 '15 at 14:40
  • $\begingroup$ @SjoerdC.deVries Indeed f idetically zero is not the interesting solution. The actual value should be something around -0.228601... which is between the other two candidates in your theory! Thank you very much for the answer $\endgroup$ – bnado Oct 17 '15 at 15:26
  • $\begingroup$ @bnado It appears that 1 works well, as I shall show in a mod to my answer soon. $\endgroup$ – bbgodfrey Oct 17 '15 at 15:29
  • $\begingroup$ @bbgodfrey yes that indeed works. Basically 1 and 0 are the "boring" solutions that I'm not interested into, while -0.22... is a physically interesting one. So for these three points I should get no divergence $\endgroup$ – bnado Oct 17 '15 at 15:49
2
$\begingroup$

The point at which the integration halts can be determined by

pf[σ_?NumericQ] := NDSolveValue[{f''[p] - 2 f[p] f'[p] + (1 + d/2) f[p] + (1 - d/2) p f'[p]
  == 0, f'[0] == Rationalize[σ, 0], f[0] == 0}, f, {p, 0, ext}, WorkingPrecision -> 30]
  ["Domain"][[1, 2]] 

pf[-.3]
(* 4.83068220440801267749372370940 *)

Addendum

Adding WorkingPrecision -> 30 and Rationalize[σ, 0] to the definition of pf proves helpful. Then,

pf[1]
(* 100.000000000000000000000000000 *)
pf[-0.228598202]
(* 10.0150679186698881369975149637 *)

So, there may be another solution very near -0.228598202, and

Quiet@FindRoot[pf[σ == 100, {σ, -0.228598202}, WorkingPrecision -> 30][[1, 2]]
(* -0.228598202437027556607595293218 *)

should find it (if it exists), but the computation is painfully slow. Moreover, the improvement is modest in this case

pf[%]
(* 10.3369620229137045868752090085 *)

So, at this level of precision the solution just found is not very good. For completeness, here is a plot of the corresponding function and its derivative

s = %%; plt = Quiet@s2[s]; Quiet@
Plot[{plt[p], plt'[p]}, {p, 0, pf[s]}, AxesLabel -> {p, "f, f'"}, PlotRange -> {-1, 7}]

enter image description here

$\endgroup$
1
$\begingroup$

One can find the value the last value in the domain of an interpolating function ifn with any of the following:

ifn["Domain"][[-1, -1]]
ifn["Grid"][[-1, -1]]
ifn["Coordinates"][[-1, -1]]

To maximize, now that we have an expression for it, we can use FindMaximum. A slight alteration of the OP's s2, similar to Sjoerd's, to make it easier to get at the InterpolatingFunction.

ClearAll[s2];
With[{kd = (2^(d - 1) π^(d/2) Gamma[d/2]), d = 3},
 s2[σ_?NumericQ] := Quiet[
    NDSolveValue[
     {f''[p] - 2 f[p] f'[p] + (1 + d/2) f[p] + (1 - d/2) p f'[p] == 0, 
      f'[0] == σ, f[0] == 0},
     f, {p, 0, Infinity}, 
     WorkingPrecision -> Precision[σ]],
    NDSolveValue::ndsz]
 ]

Machine precision will be much faster:

obj[s_?NumericQ] := s2[s]["Domain"][[-1, -1]];
{resmp, solmp} = FindMaximum[obj[σ], {σ, -0.225}]

FindMaximum::sdprec: Line search unable to find a sufficient increase in the function value with MachinePrecision digit precision. >>

(*  {10.903, {σ -> -0.228598}}  *)

We can use the machine precision result as a start for the higher precision search:

{res, sol} = FindMaximum[obj[σ], {σ, σ /. solmp}, WorkingPrecision -> 20]

FindMaximum::sdprec: Line search unable to find a sufficient increase in the function value with 20.` digit precision. >>

(*  {12.182688412696092583, {σ -> -0.22859820215629144528}}  *)

More precision, slightly longer interval, but much, much more time to compute:

{res, sol} = FindMaximum[obj[σ], {σ, σ /. solmp}, WorkingPrecision -> 30]

FindMaximum::sdprec: Line search unable to find a sufficient increase in the function value with 30.` digit precision. >>

(*
  {13.5616831717857365714489087304, 
   {σ -> -0.228598202437027556607597968677}}
*)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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