Addendum (orginial question below). Thank you for the responses! After thinking some more about this problem, and thanks to @bbgodfreys helpful comment, I realised that the problem as posted below is incomplete. What I need instead is to find $R$ and $f$ such that

  1. $f(R) = g(\lambda R)$,
  2. $f'(R) = g'(\lambda R)$,
  3. $\lim_{r \to 0} \frac{r f'(r)}{f(r)} = -\gamma$, and

$f$ also satisfies the ODE below. This leads to two questions.

Q1. Which two conditions I should pass to NDSolve, and which one should be the focus of the FindRoot operation? Further, what is the correct way to pass condition 3 to NDSolve?

Q2. In either case, the small-$r$ limit and corresponding singularity of $f$ in condition 3 makes life difficult - what is the best way to deal with this?

Perhaps Q1 is not so critical - Q2 seems to me to be the bottleneck. My current attempt, using the methods in the answers below, is to pass conditions 1 and 2 to NDSolve, and use FindRoot for condition 3. I have however been unsuccessul thus far. To wit, if I set $r = 0$ as the left boundary of NDSolve, I first of all get the expected warning

NDSolveValue::ndsz: At r == 1.2839295947759316`*^-19, step size is effectively zero; singularity or stiff system suspected.

This sets a lower bound on $r$ for FindRoot, since I assume from the warning above that I am not to trust the solution for $r < 10^{-19}$. I hence have to make an ad hoc choice for FindRoot for condition 3, for example,

trig = R/.FindRoot[
  -\[Gamma] - (r Derivative[1][solver[R]][r])/solver[R][r] /. CONS /. 
   r -> 10^-15,
  {R, (1/100/1000)/5}

(* 4.28145*10^-15 *)

I am unconvinced by FindRoot's solution, since (i) it seems to depend on my ad hoc choice for "$r \to 0$", and (ii) doesn't seem to perform well anyway:

sol = solver[trig]

-\[Gamma] - (r Derivative[1][sol][r])/sol[r] /. CONS /. r -> 10^-15
(* 0.0230768 *)

I'm aware that this addendum significantly alters the question - my apologies! Would be grateful for any suggestions on how to proceed from here. Thanks again for your help thus far.

Original question. I have the following ODE for $f(r)$:

pdef = 2 (\[Beta] f[r] + 
r (-\[Mu] + \[Rho]) Derivative[1][f][
      r] + ((-1)^(\[Gamma]/(-1 + \[Gamma])) r^((
     2 \[Gamma])/(-1 + \[Gamma])) (-1 + \[Gamma]) \
\[Zeta]^(-(\[Gamma]/(-1 + \[Gamma]))) \[Rho]^(\[Gamma]/(-1 + \
      Derivative[1][f][r]^(\[Gamma]/(-1 + \[Gamma])))/\[Gamma]) == 
 r^2 \[Sigma]^2 (f^\[Prime]\[Prime])[r]

I am looking for an $R > 0$ such that for the following known function $g$,

g[r_] := (2^(1 - \[Gamma]) (1/
  r)^\[Gamma] \[Eta]^\[Gamma] \[Rho]^-\[Gamma] ((-2 \[Beta] + \
\[Gamma] (-2 \[Mu] + 
      2 \[Rho] + (1 + \[Gamma]) \[Sigma]^2))/(-1 + \[Gamma]))^(-1 + \

we have $f(R) = g(\lambda R)$, and $f'(R) = g'(\lambda R)$ for $\lambda$ a constant. Further, for $r < R$, we must have $f(r) > g(r)$, as well as $\lim_{r \to\ 0} f(r) = \infty$

Since I do not know $R$, my idea is to wrap NDSolve in FindRoot as follows. First, for some guess $R$, define

solver[R_?NumericQ] := NDSolveValue[{
   pdef /. CONS,
   DirichletCondition[f[r] == g[\[Lambda] r] /. CONS, r == R]
  {r, R/10, R}


CONS = {\[Rho] -> 0.02, 
  J -> 5, \[Gamma] -> 0.5, \[Mu] -> 0.03, \[Sigma] -> 0.2, \[Zeta] -> 
   0.5, \[Beta] -> 0.012, k -> 10, \[Eta] -> 0.8, \[Lambda] -> 1.03}

is a set of constants. Then use FindRoot to ensure that the condition on the derivative is specified:

  Derivative[1][solver[R]][R] - D[g[\[Lambda] R], R] /. CONS,
  {R, 1/100000}

However, I run into trouble with my function solver, which gives me the error message

DiscretizePDE::femdpop: The FEMLoadElements operator failed.

A plot shows that the solution, instead of "shooting" up and above the orange curve, collapses to the $r$-axis.

enter image description here

What is wrong with my NDSolve setup? I would be grateful for any suggestions on how to improve my approach or attack the problem with an alternative method.

PS. I have resorted to the above since I simply could not understand how to implement NeumannValue for this problem. I was also unsuccessful in implementing Method -> "Shooting" with the smooth fit conditions: the solution "collapsed" just as above, leading me to believe that the derivative condition was simply being ignored when shooting.

  • 1
    $\begingroup$ Because pdef is a second order ODE, NDSolve requires that two boundary conditions be specified. One of them clearly is f[R]. Since f'[R] is to be used to determine R. it cannot be used for NDSolve. This suggests that a boundary condition very near r = 0 be specified. This, in turn, requires that an asymptotic solution be determined there, perhaps i the form of a power series. $\endgroup$
    – bbgodfrey
    Commented Jul 13, 2023 at 2:20
  • 1
    $\begingroup$ Perhaps, r*f'[r]/f[r] == -1/2 for very small r could be the inner boundary condition. $\endgroup$
    – bbgodfrey
    Commented Jul 13, 2023 at 2:32
  • $\begingroup$ @bbgodfrey You're absolutely correct that a third boundary condition is necessary; in fact, this is what makes the free-boundary problem overdetermined. But I wasn't sure how to deal numerically with $\lim_{r \to 0} f(r) = \infty$; how did you derive the boundary condition that you suggest? $\endgroup$
    – Anthony
    Commented Jul 13, 2023 at 7:46
  • 1
    $\begingroup$ Sorry for the slow reply. I shall not have much time to work on this until late today. I obtained the small-r boundary conditin by expanding the ODE there and keeping the leading term, which is proportional to 1/Sqrt[r]/ In p rinciple, you can simply insert the equation in my last comment into NDSolve. However, I encountered overflow errors after doing so, which I need to investigate. $\endgroup$
    – bbgodfrey
    Commented Jul 13, 2023 at 15:50
  • 1
    $\begingroup$ Does my answer below indicate why I have abandoned `-1/2`` as the proper limit? $\endgroup$
    – bbgodfrey
    Commented Jul 21, 2023 at 0:54

3 Answers 3


The question's ode with the given constants, rationalized, is

pdfc = 2 ((3 f[r])/250 + 25/(r^2 f'[r]) - 1/100 r f'[r]) == 1/25 r^2 f''[r]

and the function g is 200/Sqrt[r]. Answering the question is difficult, because the ode has an essential singularity at r = 0, and Limit[r f'[r]/f[r], r->0] can have multiple values.

Approximate symbolic solution of pdfc

Some insight can be gained, however, by obtaining an approximate solution to the ODE, obtained by dropping the term 25/(r^2 f'[r]).

sa = DSolveValue[2 ((3 f[r])/250 - 1/100 r f'[r]) == 1/25 r^2 f''[r], f[r], r]
(* r^(1/2 (1/2 - Sqrt[53/5]/2)) C[1] + r^(1/2 (1/2 + Sqrt[53/5]/2)) C[2] *)

The two exponents have numerical values {-0.563941, 1.06394}. If, then, f[r] is set equal to the first exponential, it is easy to show that the term 25/(r^2 f'[r]) indeed can be ignored for very small r. The same is not true of the second exponential, which should be dropped (e.g., C[2] = 0) from the approximate solution. Consequently, this solution satisfies the limit,

Limit[r f'[r]/f[r], r->0] = 1/4 (1 - Sqrt[53/5])

Interesting, an exact solution also exists, c/Sqrt[r], with {c -> -50 Sqrt[10], c -> 50 Sqrt[10]}.

Simplify[pdfc /. f -> Function[{r}, -50 Sqrt[10]/Sqrt[r]]]
(* True *)

and similarly for the second value of c. As will become apparent below, this solution is not of much value. For completeness, this solution has the small r limit,

Limit[r f'[r]/f[r], r->0] = -1/2

There may be other symbolic solutions, but I have been unable to find any.

Numerical solutions

Obtaining numerical solutions may seem straightforward: With initial conditions, {f[rm] == g[1.03 rm], f'[rm] == g'[1.03 rm]}, integrate toward r = 0. However, integrating toward an essential singularity invites numerical errors. So, I first integrated to a very small value of r0 = 10^-30 and then integrated back again to rm to ascertain whether the same values of `{f[rm], f'[rm]} resulted. It turned out that very high precision was required. For instance,

r0 = 10^-30; rm = 1; wp = 150;
s1 = NDSolveValue[{pdfc, f[rm] == 200/Sqrt[103/100 rm], 
    f'[rm] == -100/Sqrt[103/100 rm]^3}, {f[r], f[r0], f'[r0]}, {r, r0, rm}, 
    WorkingPrecision -> wp, Method -> "StiffnessSwitching", 
    InterpolationOrder -> All];
s2 = NDSolveValue[{pdfc, f[r0] == SetPrecision[s1[[2]], wp], 
    f'[r0] == SetPrecision[s1[[3]], wp]}, {f[r], f[rm], f'[rm]}, {r, r0, rm}, 
    WorkingPrecision -> wp, Method -> "StiffnessSwitching", 
    InterpolationOrder -> All];

A measure of the accuracy of the calculation is

{s2[[2]]/g[103/100 rm]) - 1, s2[[3]]/g'[103/100 rm] - 1} // N
(* {-6.2683*10^-29, 1.27023*10^-28} *)

(Smaller values of wp often perform much less well.) The curves themselves are shown by

LogLogPlot[{First[s1], g[r], First[s2]}, {r, r0, rm}, 
    AxesLabel -> {r, "f,g"}, LabelStyle -> {12, Bold, Black}]

enter image description here

where the blue and green curves (superimposed) are the inward and outward numerical integrations, and the orange curve is g[r]. More informative is a plot of r f'[r]/f[r].

LogLinearPlot[Evaluate[r D[s2[[1]], r]/s2[[1]]], {r, r0, rm}, 
    AxesLabel -> {r, "r f'/f"}, LabelStyle -> {12, Bold, Black}]

enter image description here

which converges to -0.564 for very small r, as might be expected from the symbolic analysis above. Note that the numerical calculation takes several minutes on my pc. Thus, if this is accepted as the desired r = 0 limit, then rm = 1 satisfies the question. In fact, every value of rm I have tried does so, and in each case Limit[r f'[r]/f[r], r->0] converges to the same value. Curves for rm = {10^-5, 1, 10^5} are

enter image description here

In no case did r f'[r]/f[r] converge to -1/2.

  • $\begingroup$ Thanks so much for this comprehensive and clear answer. Do I understand correctly that your conclusion is the same as Alex's above, namely, that the problem as posed does not admit a unique solution? $\endgroup$
    – Anthony
    Commented Jul 25, 2023 at 7:35
  • 1
    $\begingroup$ @Anthony Yes, that is correct. $\endgroup$
    – bbgodfrey
    Commented Jul 25, 2023 at 15:08
  • $\begingroup$ Thanks, this whole discussion has been enlightening. I'll have to think hard about my problem in order to understand what condition I need to impose in order to obtain a unique solution. $\endgroup$
    – Anthony
    Commented Jul 25, 2023 at 17:30

I'm more accustomed to working in a slightly different format, but this should be equivalent to what you're doing. Define functions to build the pdes and g:

pdef[\[Beta]_, \[Gamma]_, \[Sigma]_, \[Mu]_, \[Rho]_, \[Zeta]_] := 
  2 \[Beta] f[r] + 
    2 r (-\[Mu] + \[Rho]) f'[
      r] + (2 (-1)^(\[Gamma]/(-1 + \[Gamma])) r^((2 \[Gamma])/(-1 + \
\[Gamma])) (-1 + \[Gamma]) \[Zeta]^(-(\[Gamma]/(-1 + \[Gamma]))) \
\[Rho]^(\[Gamma]/(-1 + \[Gamma])) f'[
         r]^(\[Gamma]/(-1 + \[Gamma])))/\[Gamma] - 
    r^2 \[Sigma]^2 f''[r] == 0;

g[\[Beta]_, \[Gamma]_, \[Mu]_, \[Eta]_, \[Rho]_, \[Sigma]_, \
   r_] := (2^(1 - \[Gamma]) (1/
        r)^\[Gamma] \[Eta]^\[Gamma] \[Rho]^-\[Gamma] ((-2 \[Beta] + \
\[Gamma] (-2 \[Mu] + 
             2 \[Rho] + (1 + \[Gamma]) \[Sigma]^2))/(-1 + \
\[Gamma]))^(-1 + \[Gamma]))/\[Gamma];

The solver with most of the constants plugged in. (You don't seem to use J or k for anything as far as I can tell.)

solver[rSol_?NumericQ] := 
  NDSolve[{pdef[0.012, 1/2, 0.2, 0.03, 0.02, 1/2], 
    f[rSol] == g[0.012, 1/2, 0.03, 0.8, 0.02, 0.2, 1.03][rSol], 
    Derivative[1][f][rSol] == 
     Derivative[1][g[0.012, 1/2, 0.03, 0.8, 0.02, 0.2, 1.03]][rSol]}, 
   f, {r, rSol/10, rSol}];

The main difference is that I removed the Dirichlet code and just used a simple f[rSol] = g[..][rSol] and the corresponding derivatives directly, using the conditions you wrote in the text.

FindRoot complained about the solution being non-numeric, so I created an auxillary function to keep it happy:

rsf[rs_?NumericQ] := 
  Derivative[1][f][rs] /. Flatten@solver[rs];

I changed the constant to a value different from yours, to verify that it did indeed search for the root:

 rsf[rSol] - 
  Derivative[1][g[0.012, 1/2, 0.03, 0.8, 0.02, 0.2, 1.03]][
   rSol], {rSol, .0002(*1/100000*)}]

Out[39]= {rSol -> 0.00001}

I reran once last time to generate plots and verify other conditions. The output was of the form {f->InterpolatingFunction[..]}, which I won't copy here since it's long and ugly.

nds = Flatten@solver[.00001]

Here are f,g, and f-g:

Plot of f,g,f-g

Since the scales are quite different, here's f-g alone:

Plot of f-g alone

Here's f'[r]


And here are {{f[rSol],g[rSol]}{f'[rSol],g'[rSol]}}. f,g and f',g' are equal at rSol, as required.

{f[r] /. nds, 
   g[0.012, 1/2, 0.03, 0.8, 0.02, 0.2, 1.03][r]}, {f'[r] /. nds, 
   D[g[0.012, 1/2, 0.03, 0.8, 0.02, 0.2, 1.03][r], r]}} /. 
 r -> 9.999999999999999`*^-6

Out[53]= {{63245.6, 63245.6}, {-3.16228*10^9, -3.16228*10^9}}
  • $\begingroup$ Thanks for this helpful answer! I need to play with it for a little bit. Your method of simply using Derivative[1][f]... in NDSolve already helps a great deal. One comment: in your definition of rsf, I assume that the 10^-5 should be replaced with rs? $\endgroup$
    – Anthony
    Commented Jul 13, 2023 at 11:28
  • $\begingroup$ Thanks again for your response. I've updated the question, which was incompletely posed; would be great if you had any further ideas. $\endgroup$
    – Anthony
    Commented Jul 14, 2023 at 12:43
  • $\begingroup$ When I run your code, FindRoot returns as an answer whatever initial guess I give it, as it must for the definition of solver. in your code. $\endgroup$
    – bbgodfrey
    Commented Jul 14, 2023 at 16:55

We can map $(0, R)$ on $(0,1)$ using substitution $r\rightarrow R r$, and use ParametricNDSolveValue as follows

eq = 2 (\[Beta] f[r] + 
      r (-\[Mu] + \[Rho]) f'[
        r] + ((-1)^(\[Gamma]/(-1 + \[Gamma])) (R r)^((2 \[Gamma])/(-1 \
+ \[Gamma])) (-1 + \[Gamma]) \[Zeta]^(-(\[Gamma]/(-1 + \[Gamma]))) \
\[Rho]^(\[Gamma]/(-1 + \[Gamma])) (f'[r]/
            R)^(\[Gamma]/(-1 + \[Gamma])))/\[Gamma]) - 
   r^2 \[Sigma]^2 f''[r];

CONS = {\[Rho] -> 0.02, 
  J -> 5, \[Gamma] -> 0.5, \[Mu] -> 0.03, \[Sigma] -> 0.2, \[Zeta] -> 
   0.5, \[Beta] -> 0.012, k -> 10, \[Eta] -> 0.8, \[Lambda] -> 1.03}; 
g[r_] := (2^(1 - \[Gamma]) (1/R/
        r)^\[Gamma] \[Eta]^\[Gamma] \[Rho]^-\[Gamma] ((-2 \[Beta] + \
\[Gamma] (-2 \[Mu] + 
             2 \[Rho] + (1 + \[Gamma]) \[Sigma]^2))/(-1 + \
\[Gamma]))^(-1 + \[Gamma]))/\[Gamma] /. CONS;

sol = ParametricNDSolveValue[{eq == 0, f[1] == g[\[Lambda]], 
    f'[1] == \[Lambda] g'[\[Lambda]]} /. CONS, f, {r, 10^-3, 1}, {R}];

Visualization numerical solution (solid line) and function $g(\lambda r)$ (dashed red line) for R=0.5,1,1.5,2

  LogLogPlot[sol[R][r], {r, 10^-3, 1}, PlotStyle -> Blue, 
   PlotLabel -> R, Frame -> True, GridLines -> Automatic, 
   FrameTicks -> {{Automatic, None}, {{0.001, 1}, Automatic}}], 
  LogLogPlot[g[\[Lambda] r] /. CONS, {r, 10^-3, 1}, 
   PlotStyle -> {Red, Dashed}]], {R, .5, 2, .5}]

Figure 1

Update 1. To satisfier condition at $r\rightarrow 0$ we suppose that $f(r)=a/r^{\gamma}$, and using eq we can define a as follows

rule = {f[r] -> a/r^\[Gamma] , f'[r] -> -\[Gamma]  a/r^(\[Gamma] + 1),
    f''[r] -> \[Gamma]  (\[Gamma] + 1) a/r^(\[Gamma] + 2)};

eq /. rule /. \[Gamma] -> 1/2 // FullSimplify; sl = Solve[% == 0, a];

Then we define function

f0 = a/r^\[Gamma] /. sl[[2]] /. CONS

Out[]= 158.114/(r^0.5 Sqrt[R])

Analyzing picture above we can suggest that function f0+b, where b is some constantan, should correlate with solution sol. To verifier this suggestion we use NMinimize as follows

eqn = (f0 + b - sol1[R][r]) /. r -> 10^-3;

 NMinimize[{eqn^2, {R > 0}}, {R, b}]

(*Out[]= {2.42103*10^-14, {R -> 138.914, b -> 215.84}}*)
  • $\begingroup$ Thanks for your response! I've updated the question, which was incompletely posed; would be great if you had any further ideas. $\endgroup$
    – Anthony
    Commented Jul 14, 2023 at 12:44
  • $\begingroup$ @Anthony Please, see update to my answer. $\endgroup$ Commented Jul 15, 2023 at 3:38
  • $\begingroup$ Thanks for your effort, and sorry for taking a few days to get back to you. I'm just parsing through the logic of your answer, and have a rather simple question to begin with: why is the mapping $r \mapsto Rr$ useful? I am also wondering how exactly you obtained eq: did you do something like pdef /. r -> Rr? Because if so, I get a different expression for eq than you do. $\endgroup$
    – Anthony
    Commented Jul 17, 2023 at 14:27
  • $\begingroup$ @Anthony Could you show your expression eq? $\endgroup$ Commented Jul 18, 2023 at 3:16
  • 1
    $\begingroup$ @Anthony For my opinion it means that your problem has no unique solution. $\endgroup$ Commented Jul 18, 2023 at 12:42

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.