4
$\begingroup$

NDSolve Results

On the course of addressing question 181974, I encountered the following problem.

k2 = Rationalize[ 0.000194519, 0];
ϵ = 10^-4;
R = Rationalize[1.5472, 0];
s = NDSolveValue[{y''[r] + 2 y'[r]/r == k2 Sinh[y[r]], y'[R] == 0, y'[ϵ] == 0, 
    WhenEvent[r == 1, y'[r] -> y'[r] + 1]}, {y, y'}, {r, ϵ, R}, 
    Method -> {"Shooting", "StartingInitialConditions" -> {y[ϵ] == 
    Rationalize[-8.9941, 0], y'[ϵ] == 0}, Method -> "StiffnessSwitching"}, 
    WorkingPrecision -> 30];
Plot[First[s][r], {r, ϵ, R}, AxesLabel -> {r, y}, 
    ImageSize -> Large, LabelStyle -> {Black, Bold, Medium}]

enter image description here

which visibly is incorrect, because y'[R] == 0 is far from being satisfied. Moreover, this error is robust in that NDSolve yields the same solution for all credible initial guesses for y[ϵ] that I have tried. Interestingly, y[ϵ] == 0, as in the plot, would be correct, if the jump condition in y'[r] at r == 1 were not included. So, my question is, have I made an error? And, if not, is this a bug? For completeness, the correct solution is

ps = ParametricNDSolveValue[{y''[r] + 2 y'[r]/r == k2 Sinh[y[r]], 
    y[ϵ] == y0, y'[ϵ] == 0, WhenEvent[r == 1, y'[r] -> y'[r] + 1]}, {y, y'}, 
    {r, ϵ, R}, {y0}, Method -> "StiffnessSwitching", WorkingPrecision -> 30];
sol = FindRoot[Last[ps[y0]][R], {y0, -9, -8.8}, Evaluated -> False][[1, 2]]
(* -8.99411 *)
Plot[First[ps[sol]][r], {r, ϵ, R}, AxesLabel -> {r, y}, 
    ImageSize -> Large, LabelStyle -> {Black, Bold, Medium}]

enter image description here

Addendum: DSolve Results

In comments below, Alex Trounov observed that NDSolve also gives an incorrect result for the linearized version of the ODE, y''[r] + 2 y'[r]/r == k2 y[r]. (The plotted result is nearly identical to the first figure above.) For completeness, the correct result is

enter image description here

The linearized ODE also can be solved symbolically, and I have done so.

DSolveValue[{y''[r] + 2 y'[r]/r == k2 y[r], y'[R] == 0, y'[0] == 0, 
    WhenEvent[r == 1, y'[r] -> y'[r] + 1]}, y[r], {r, 0, R}]
(* Piecewise[{{0, 0 <= r <= 1}, {(5000*Sqrt[10/1151]*(-E^((13*Sqrt[1151/10])/5000)
   + E^((13*Sqrt[1151/10]*r)/5000)))/(13*E^((13*Sqrt[1151/10]*(1 + r))/10000)*r), 
   Inequality[1, Less, r, LessEqual, 967/625]}}, Indeterminate] *)

Plotting this expression again yields a plot nearly identical to the first one in this question. In other words, DSolve also yields the same wrong answer! To obtain the correct symbolic answer, begin by noting that the solution of the ODE without boundary conditions is of the form, c1 Sinh[Sqrt[k] r]/r + c2 Cosh[Sqrt[k] r]/r, as can be verified by,

Simplify[(y''[r] + 2 y'[r]/r == k y[r]) /. y -> 
    Function[r, c1 Sinh[Sqrt[k] r]/r + c2 Cosh[Sqrt[k] r]/r], k > 0]
(* True *)

The solution for 0 < r < 1 must be

fm[r_] = c1 Sinh[Sqrt[k] r]/r

so that it is not singular at the origin. The solution for 1 < r < rmax is

fp[r_] = c3 Sinh[Sqrt[k] r]/r + c4 Cosh[Sqrt[k] r]/r

Applying matching conditions at r == 1 and the y'[rmax] == 0 outer boundary condition determines the three coefficients.

sc = Solve[{fm[1] == fp[1], fm'[1] + 1 == fp'[1], fp'[rmax] == 0}, 
    {c1, c3, c4}] // Simplify // Flatten
(* {c1 -> (-Sqrt[k] rmax Cosh[Sqrt[k] (-1 + rmax)] + Sinh[Sqrt[k] (-1 + rmax)])
       /(k rmax Cosh[Sqrt[k] rmax] - Sqrt[k] Sinh[Sqrt[k] rmax]), 
    c3 -> (Sinh[Sqrt[k]] (Cosh[Sqrt[k] rmax] - Sqrt[k] rmax Sinh[Sqrt[k] rmax]))
       /(-k rmax Cosh[Sqrt[k] rmax] + Sqrt[k] Sinh[Sqrt[k] rmax]), 
    c4 -> -(Sinh[Sqrt[k]]/Sqrt[k])} *)

Consequently, the complete symbolic solution is

sa = Piecewise[{{fm[r], r < 1}, {fp[r], r > 1}}] /. sc /. {rmax -> R, k -> k2}

which is a bit long to reproduce here. Plotting it yields a figure identical to that immediately above. In conclusion, neither DSolve nor NDSolve is able to solve this relatively simple problem, in both cases returning essentially the same wrong answer. Certainly, this is a bug. I would welcome the perspectives of others.

$\endgroup$
  • 1
    $\begingroup$ In your case, ODE is non-linear, but I checked that the same error remains for the linear equation `s = NDSolve[{y''[r] + 2 y'[r]/r == k2*y[r], y'[R] == 0, y'[[Epsilon]] == 0, WhenEvent[{r == 1, y'[r] < 0}, y'[r] -> y'[r] + 1]}, y, {r, [Epsilon], R}, Method -> {"Shooting", "StartingInitialConditions" -> {y[[Epsilon]] == Rationalize[-8.9941, 0], y'[[Epsilon]] == 0}, Method -> "StiffnessSwitching"}, WorkingPrecision -> 30]' $\endgroup$ – Alex Trounev Sep 18 '18 at 5:00
  • 2
    $\begingroup$ You can remove all methods and solve a linear equation y''[r] + 2 y'[r]/r == k2 y[r], there remains a similar error. It seems that the boundary condition y'[R] == 0 is ignored in combination with WhenEvent, $\endgroup$ – Alex Trounev Sep 18 '18 at 5:58

Your Answer

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

Browse other questions tagged or ask your own question.