2

I am trying to use ParametricNDSolve to find a desired value for a parameter in an ODE, and the output seems to misbehave.

First, consider the following example given in the documentation:

sol1 = ParametricNDSolve[{x''[t] - x'[t] + x[t] == 0, x[0] == 1, x'[0] == s}, x, {t, 0, 30}, {s}]
(* ::Output:: *) (*
{x -> ParametricFunction[ <> ]}
*)
root1 = FindRoot[Evaluate[x[s][10] /. sol1], {s, 6}]
(* ::Output:: *) (*
{s -> 1.40296}
*)

Secondly, my equation with the Automatically chosen method:

sol2 = ParametricNDSolve[{X'[t] == t + (Sqrt[3] / Pi) Log[a/(1 - a)] Abs[t X[t]] - X[t], X[0] == 1}, {X}, {t, 0, 2}, a]
(* ::Output:: *) (*
{x -> ParametricFunction[ <> ]}
*)

root2 = FindRoot[Evaluate[X[a][2] /. sol2] == 3, {a, .5}]
(* ::Output:: ) (
[!] InterpolatingFunction: Input {2} lies outside the range of data... (InterpolatingFunction::dmval)
{a -> 0.812787}
*)

In which I don't understand how the range is exceeded. Furthermore, specifying a method in the call yields:

sol3 = = ParametricNDSolve[{X'[t] == t + (Sqrt[3] / Pi) Log[a/(1 - a)] Abs[t X[t]] - X[t], X[0] == 1}, {X}, {t, 0, 2}, a, Method -> {"TimeIntegration" -> "ExplicitEuler"}]
(* ::Output:: *) (*
{x -> ParametricFunction[ <> ]}
*)

X[.7][2] /. sol3
(* ::Output:: ) (
2.01701
*)


root3 = FindRoot[Evaluate[X[a][2] /. sol3] == 3, {a, .5}]
(* ::Output:: ) (
[!] ParametricNDSolve: Encountered invalid NDSolveSensitivityMethod method data object at point t=0. (ParametricNDSolve::mdata)
[!] FindRoot: At {a}={0.5}, function value {-3+ParametricFunction[1,InternalBag[&lt;1&gt;],1,1,False,{{a$71966},&lt;&lt;5&gt;&gt;,{0}},{NDSolvebase$71973,NDSolveNDSolveParametricFunction[0,{ParametricNDSolve,InternalBag[<2>],None,ParametricNDSolve},{{{<<9>>},{<<9>>}},{0,{<<3>>},{<<2>>},{<<3>>}},None,{{<<7>>},{<<1>>},None,{}}},{X},<<4>>,{Cache->True,CacheTableLength->19,CacheTableWidth->7,CacheKeyMaxBytes->1000000,CacheResultMaxBytes->1000000,KeyComparison->None,ResultComparison->LessEqual},{},<<1>>]}][0.5][2]} is not a {1} dimensional list of numbers  (FindRoot::nlnum)
{a -> 0.5}
*)


X[.7][2] /. sol3
(* ::Output:: ) (
ParametricFunction[ <> ][0.7][2]
*)

Note that t=2 is indeed in the interpolation range, and the call to the ParametricFunction instance fails after the unsuccessful FindRoot run.

How to understand this behavior, and is there a current workaround if it's a bug?

Related:

  1. ParametricFunction from ParametricNDSolveValue changes when evaluated?
  2. Vector ParametricNDSolve and FindRoot interaction
  3. Issue in ParallelTable after evaluating another function using NDSolve and FindRoot
Gravifer
  • 862
  • 5
  • 18
  • In the first link, @Carl Woll's answer suggests that there is indeed something wrong with ParametricFunction, however I am not sure if the behavior in the current post is of the same origin. Root solving is a large numeric topic, and I haven't finished all the tech notes and SE QA's just yet; do be kind enough and point me to anything that I am supposed to read. – Gravifer Dec 14 '22 at 11:42
  • 1
    I think "Findroot" makes a first guess using something like the tangent. This can lead to a bad first guess. Therefore, try a starting value closer to the searched root. E.g.: FindRoot[Evaluate[X[a][2] /. sol1] == 3, {a, .7}] – Daniel Huber Dec 14 '22 at 12:17
  • @DanielHuber This indeed yields a good root and prevents the ParametricFunction from being spoiled – Gravifer Dec 14 '22 at 13:06
  • 1
    To get a good starting value use Plot, i.e., Plot[(X[a][2] - 3) /. sol3, {a, 0, 1}] indicates the solution is in the vicinity of 0.8 Alternatively, avoid derivatives by using "FindRoot[lhs == rhs, {x, Subscript[x, 0], Subscript[x, 1]}] searches for a solution using Subscript[x, 0] and Subscript[x, 1] as the first two values of x, avoiding the use of derivatives." – Bob Hanlon Dec 14 '22 at 15:27

3 Answers3

4

Evaluated -> False

Clear[sol2, root2];
sol2 = ParametricNDSolve[{X'[t] == 
    t + (Sqrt[3]/Pi) Log[a/(1 - a)] Abs[t X[t]] - X[t], 
   X[0] == 1}, {X}, {t, 0, 2}, a]
root2 = FindRoot[(X[a][2] /. sol2) == 3, {a, .5}, Evaluated -> False]

{a -> 0.812787}

cvgmt
  • 72,231
  • 4
  • 75
  • 133
  • I find this indeed works. Would you explain what the option means and does? Thanks in advance – Gravifer Dec 14 '22 at 15:27
  • 1
    Clear[sol]; sol = ParametricNDSolveValue[{X'[t] == t + (Sqrt[3]/Pi) Log[a/(1 - a)] Abs[t X[t]] - X[t], X[0] == 1}, {X}, {t, 0, 2}, {a}]; FindRoot[Indexed[sol[a], 1][2] == 3, {a, .5}] – cvgmt Mar 12 '23 at 06:23
2

and is there a current workaround

Some strange evaluation order? Here is a workaround

sol2 = ParametricNDSolve[{X'[t] == t + (Sqrt[3]/Pi) Log[a/(1 - a)] Abs[t X[t]] - X[t], 
    X[0] == 1}, {X}, {t, 0, 2}, a];
f[a_] := X[a] /. sol2
root2 = FindRoot[Evaluate[f[a][2] /. sol2] == 3, {a, .5}]

Mathematica graphics

No beep and no warnings and no range exceeded messages.

V 13.1 on windows 10

Screen shot. No messages. Mathematica graphics

Nasser
  • 143,286
  • 11
  • 154
  • 359
  • v13.1 here too. Weirdly enough this doesn't seems to fix it for me. I'll clear the kernel and try again. – Gravifer Dec 14 '22 at 12:57
  • @Gravifer no messages/warnings for me on V 13.1 on windows 10. Yes, make sure to Quit the kernel. Added screen shot – Nasser Dec 14 '22 at 13:21
  • I can now replicate your result. However, what dictates the result for me is /.sol2 in the FindRoot; with this replacement no exception is thrown. I don't fully understand why this is the case yet. – Gravifer Dec 14 '22 at 15:24
2

Without Method -> {"TimeIntegration" -> "ExplicitEuler"} try more direct solution using ParametricNDSolveValue(only return X[2])

X3 = ParametricNDSolveValue[{X'[t] == t + (Sqrt[3]/Pi) Log[a/(1 - a)] Abs[t X[t]]- X[t], X[0] == 1}, 
X[2] , {t, 0, 2}, a  ]
root3 = FindRoot[ X3[a]   == 3, {a, 1/2 }]
(* {a -> 0.812787} *)

NMinimize evaluates without message " Input value {2} lies outside..."

NMinimize[ (X3[a]  - 3)^2, a]
(*{6.10426*10^-16, {a -> 0.812787}}  *)
Ulrich Neumann
  • 53,729
  • 2
  • 23
  • 55