# FindMaximum not working with ParametricNDSolve

I have the following code (ignore the manipulate parameters, assume them constant in this context).

For background, what I'm trying to do is take the results of a 2D differential equation, find where it intersects y=0, and maximize the x where that happens for different parameter values.

If I replace FindMaximum with Plot (and change {R, 0.05} to {R, 0.005, 0.05}) I get a perfect plot of the x intercept as a function of R, with a visible maximum around R=0.02. However, the below code throws a lot of errors, including

FindRoot::nlnum: The function value {ParametricFunction[................][R]} is not a list of numbers with dimensions {1} at {t} = {40.}.
ReplaceAll::reps: {FindRoot[(y[t]/. sol$40848)[R],{t,40,0.05,50}]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing. FindMaximum::nrnum: The function value -<<1>> is not a real number at {R} = {0.05}.  In some iterations of this, when I was using NDSolve, I was actually getting in an infinite loop where the kernel would hang on spewing repetitive error messages (mostly about R not being replaced in NDSolve). What's going on? ρm = 1600 ;(*Density of rock*) ρa = 1.225;(*Density of air*) Cd = 0.5;(*Typical drag coefficient for a somewhat rough sphere*) g = 9.8;(*Gravity*) Manipulate[ Module[{sol, v0}, v0 = Sqrt[2 A/(ρm 4/3 π r^3 + B)]; sol = ParametricNDSolve[{ρm 4/ 3 π r^3 y''[t] == - ρm 4/3 π r^3 g - 1/2 ρa π r^2 Sqrt[x'[t]^2 + y'[t]^2] Cd y'[t], ρm 4/ 3 π r^3 x''[t] == -(1/2) ρa π r^2 Sqrt[ x'[t]^2 + y'[t]^2] Cd x'[t], x[0] == 0, y[0] == 0, x'[0] == v0 Cos[θ], y'[0] == v0 Sin[θ]}, {x[t], y[t]}, {t, 0, 100}, { r}]; FindMaximum[(((x[t] /. sol)[ R]) /. {FindRoot[(y[t] /. sol) [R], { t, 40, 0.05, 50}]})[[1]], {R, 0.05}] ], {{θ, π/4}, 0, π/2}, {{A, 149, "A [Throwing Ability] (J)"}, 0.0001, 200, Appearance -> {"Open"}}, {{B, 0.5, "B [~1/3 arm mass] (kg)"}, 0.0001, 2, Appearance -> {"Open"}}, Item[StringForm[ "\!$$\*SubscriptBox[\(ρ$$, $$m$$]\) = \ kg/\!$$\*SuperscriptBox[\(m$$, $$3$$]\), \ \!$$\*SubscriptBox[\(ρ$$, $$a$$]\) = \ kg/\!$$\*SuperscriptBox[\(m$$, $$3$$]\), \!$$\*SubscriptBox[\(C$$, \ $$d$$]\) = ", ρm, ρa // StandardForm, Cd // StandardForm]]]  • Parameter v0 isn't defined! – Ulrich Neumann Oct 13 '20 at 8:18 • Define a function for the root search outside of FindMaximum e.g.: funy[R_] := (((x[t] /. sol)[ R]) /. {FindRoot[{Evaluate[(y[t] /. sol)[R]] == 0}, {t, 40, 0.05, 50}]})[[1]]; Then you can print funy[R] and see what is does. Then you will see that the root search failes because it can not find a root in the interval given. – Daniel Huber Oct 13 '20 at 9:41 • Use WhenEvent  and Reap  and Sow[x[t]]  to stop integration of NDSolve at y[t]==0 and to get the x[t] at that t as the last value in the reap list. Plot the x values you then get as a function of r with the other parameters varied as you suggested with Manipulate. – Akku14 Oct 13 '20 at 10:36 • @Daniel I did, that works fine. Like I said it Just Works when I plot it – Manishearth Oct 13 '20 at 14:27 • Hi, if I do as you said, replacing FindMaximium by Plot and changing the range, I get the same erro message: "FindRoot::nlnum: The function value {ParametricFunction[1,InternalBag[<1>],2,1,{{r$6738},SystemUtilitiesHashTable[<2>],{},{},{1},{Automatic,0,0}},{NDSolvebase\$6747,NDSolveNDSolveParametricFunction[0,{ParametricNDSolve,InternalBag[<2>],None,ParametricNDSolve},<<6>>,{},All]}][R]} is not a list of numbers with dimensions {1} at {t} = {40.}. " – Daniel Huber Oct 13 '20 at 15:11

There seems to be two problem. FindMaximum has the attribute "HoldAll" and this prevents the replacements. Plot does not have this attribute, what eliminates the problem. What exactly is going on, I was not able to find out. But when you declare a function with the replacements outside of FindMaximum it works. The second problem is with the start value of R. It is too far off and the algorithm goes to the wrong side.

I post the changed code below.

ρm = 1600;(*Density of rock*)ρa = 1.225;(*Density of \
air*)Cd = 0.5;(*Typical drag coefficient for a somewhat rough \
sphere*)g = 9.8;(*Gravity*)Manipulate[
Module[{sol, v0}, v0 = Sqrt[2 A/(ρm 4/3 π r^3 + B)];
sol = ParametricNDSolve[{ρm 4/3 π r^3 y''[
t] == -ρm 4/3 π r^3 g -
1/2 ρa π r^2 Sqrt[x'[t]^2 + y'[t]^2] Cd y'[
t], ρm 4/3 π r^3 x''[
t] == -(1/2) ρa π r^2 Sqrt[x'[t]^2 + y'[t]^2] Cd x'[
t], x[0] == 0, y[0] == 0, x'[0] == v0 Cos[θ],
y'[0] == v0 Sin[θ]}, {x[t], y[t]}, {t, 0, 100}, {r}];
fun[R_ /;
NumericQ[
R]] := (((x[t] /. sol)[
R]) /. {FindRoot[(y[t] /. sol)[R], {t, 40, 0.05,
50}]})[[1]] /. {θ -> π/4, A -> 149, B -> 0.5};
FindMaximum[fun[R], {R, 0.01}]
], {{θ, π/4},
0, π/2}, {{A, 149, "A [Throwing Ability] (J)"}, 0.0001, 200,
Appearance -> {"Open"}}, {{B, 0.5, "B [~1/3 arm mass] (kg)"},
0.0001, 2, Appearance -> {"Open"}},
Item[StringForm[
"\!$$\*SubscriptBox[\(ρ$$, $$m$$]\) = \
kg/\!$$\*SuperscriptBox[\(m$$, $$3$$]\), \
\!$$\*SubscriptBox[\(ρ$$, $$a$$]\) = \
kg/\!$$\*SuperscriptBox[\(m$$, $$3$$]\), \!$$\*SubscriptBox[\(C$$, \
$$d$$]\) = ", ρm, ρa // StandardForm,
Cd // StandardForm]], TrackedSymbols -> {θ, A, B}]