# DSolve, NDSolve with WhenEvent Give Incorrect Solution for Simple ODE

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


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


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

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.

• 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]' – Alex Trounev Sep 18 '18 at 5:00
• 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`, – Alex Trounev Sep 18 '18 at 5:58