FindRoot and NDSolve question with multiple parameters

I'm trying to reproduce the results from 1, which solve the equations of an elastic ring (a closed loop elastica) under various loadings. Here is the relevant part:

I need to solve the governing equations in (2), with m0 and p as unknown parameters for a range of prescribed values of f. They indicate in the text that they solve these equations using Mathematica with NDSolve and FindRoot.

I found the following question helpful as an example: How can I use FindRoot on an expression from NDSolve?, and so I implemented my code in a similar manner. However, I'm trying to use FindRoot to solve for two parameters, and I'm not sure how to use the results of NDSolve to get the necessary two equations. Here is my code:

sol[p_?NumericQ, m0_?NumericQ, f_?NumericQ] := {\[Theta], x, y, m} /.
First@NDSolve[{
x'[s] == Cos[\[Theta][s]],
y'[s] == Sin[\[Theta][s]],
\[Theta]'[s] == m[s] + 2 \[Pi],
m'[s] == f/2 Cos[\[Theta][s]] - p Sin[\[Theta][s]],
x[0] == y[0] == \[Theta][0] == 0,
m[0] == m0
},
\[Theta], {s, 0, 2 \[Pi]}]

FindRoot[{
sol[p, m0, 200][[1]][1/2] == \[Pi],
sol[p, m0, 200][[2]][1/2] == 0}, {p, 0}, {m0, -10}]


FindRoot gives me an error that states:

FindRoot: The function value {-3.14159+0.[0.5],(-10)[0.5]} is not a list of numbers with dimensions (2) at {p,m0}={0.,-10.}.


It seems to me that FindRoot is evaluating things in the wrong order. This feels like the type of Mathematica question that is answered with "Use ?NumericQ", for example, however I am doing that here, so I'm confused. I have tried wrapping Evaluate around the functions in that list, but so far I've had no luck.

I suspect the issue is my lack of knowledge of how to call sol and use it in FindRoot in the presence of ?NumericQ. For example, it works just fine, and I can visually see the roots when I plot:

Plot3D[{sol[p, m0, 200][[1]][1/2],
sol[p, m0, 200][[2]][1/2], \[Pi]}, {p, -1, 1}, {m0, -50, 50}]


So clearly, the information is there, I just can't access it correctly yet. Any help would be greatly appreciated.

1 L.N. Virgin et al., "Deformation and vibration of compressed, nested, elastic rings on rigid base", Thin-Walled Structures, 132, 167-175, (2018). Link (Paywall)

In your definition of sol you should list all functions {\[Theta], x, y,m}!

Change the definition of sol to

sol[p_?NumericQ, m0_?NumericQ, f_?NumericQ] :=NDSolveValue [{
x'[s] == Cos[\[Theta][s]],y'[s] == Sin[\[Theta][s]],
\[Theta]'[s] == m[s] + 2 \[Pi],
m'[s] == f/2 Cos[\[Theta][s]] - p Sin[\[Theta][s]],
x[0] == y[0] == \[Theta][0] == 0, m[0] == m0}, {\[Theta], x, y,m}, {s,0, 2 \[Pi]} ]


CounturPlot now shows possible solutions

ContourPlot[{sol[p, m0, 200][[1]][1/2] == \[Pi],sol[p, m0, 200][[2]][1/2] == 0}, {p, -1 , 1 }, {m0, -15 , -10},FrameLabel -> {p, m0}]


With explicitly defined equations,

gln[p_?NumericQ, m0_?NumericQ] := {sol[p, m0, 200][[1]][1/2] - \[Pi],sol[p, m0, 200][[2]][1/2] }


FindRoot is able to solve

FindRoot[  gln[p, m0] , {{p, 0, -1, 1}, {m0, -10, -15, -5}} ]
(*{p -> 3.80678*10^-6, m0 -> -11.9541}*)


addendum Option Evaluated->False also works

FindRoot[  {sol[p, m0, 200][[1]][1/2] - \[Pi],sol[p, m0, 200][[2]][1/2] - 0 } , {{p, 0, -1,1}, {m0, -10, -15, -5}} , Evaluated -> False]


without predefinition!

• Thanks. I'm also able to visualize the roots, however I'm unable to get FindRoot to determine them. This is important because I plan to iterate of f, instead of just finding the roots at f=200. But, FindRoot seems to call Part in a different order than Plot or Plot3D. I think this question has a similar issue, but I haven't been able to implement it correctly yet: mathematica.stackexchange.com/questions/34355/… Mar 3, 2020 at 19:49
• @dpholmes I completed my answer. With predefined equations FindRoot is able to solve. Mar 3, 2020 at 20:09
• ah-ha! the functions need to be explicitly defined before passing them to FindRoot. Thank you! Mar 3, 2020 at 20:17
• @dpholmes You're welcome. Predefinition seems to prevent FindRoot from "SymbolicProcessing". Mar 3, 2020 at 20:21