differential equations - Limitations of ParametricNDSolve family w.r.t objective functions

Observation:

I can see even for very simple modification in case of an scalar objective involving an definite integral in time ParametricNDSolve fails. Here is an example!

eqns = {y''[t] + y[t] == 3 a Sin[y[t]], y[0] == y'[0] == 1};

pfun =ParametricNDSolveValue[eqns,Integrate[y[s] - a s, {s, 0, 5}], {t, 0, 5}, {a},
Method -> {"ParametricSensitivity" -> "ForwardSensitivity"}];
pfun[1.5]

Meaningless output!

Same kind of output for pfun'[1.5] but from pfun''[1.5] onwards for higher derivatives we get numerical values which I guess are totally wrong.

However everything will be fine if one uses Integrate[y[s], {s, 0, 5}]! So I tried {"ParametricSensitivity" -> "AdjointSensitivity"} which is most suitable for integrated objective functions. Again failure but this time for both the cases. We get the following error

ParametricNDSolveValue::adjsens: The adjoint sensitivity method cannot be used for the output function {t,0,5}. It can only be used for output functions that are at a particular time or are a definite integral over time. >>

I feel this is a major inconsistency of implementation internal. Using Trace I found some esoteric IntegrateImproperDump` andInternalDependsOnQ.

What should be pfun[1.5]:

We know

Distribute[Integrate[y[s] - a s, {s, 0, 5}], Plus] === 
Integrate[-a s, {s, 0, 5}] + Integrate[y[s], {s, 0, 5}]

True

So we first can find pfun[1.5] using

eqns = {y''[t] + y[t] == 3 a Sin[y[t]], y[0] == y'[0] == 1};
pfun = ParametricNDSolveValue[eqns, 
Integrate[y[s], {s, 0, 5}], {t, 0, 5}, {a}, 
Method -> {"ParametricSensitivity" -> "ForwardSensitivity"}];
(pfun[1.5] + Integrate[-a s, {s, 0, 5}]) /. a -> 1.5

-7.86673

and the first order sensitivity will be

(pfun'[1.5] + D[Integrate[-a s, {s, 0, 5}], a]) /. a -> 1.5

-7.87591

Crosschecking the first order sensitivity below!

fun1[aval_?NumericQ] :=NIntegrate[Block[{a = aval}, 
Evaluate@(y /. First@NDSolve[Evaluate@eqns, y, {t, 0, 5}])[t] -a t] , {t, 0, 5}];

Needs["NumericalCalculus`"];
ND[fun1[x], x, 1.5]

-7.87604

Pretty much as expected.

Question:

• It will be great to know if we can use ParametricNDSolve family to find parameter dependency of integrated objective like the following: $$ G(p)=\int_a^{b} g(y(s),s,p) \,ds$$ where $g$ is a function of the dependent variable $y(s)$ of the underlying differential equation system and $p$ represents the parameter with respect to which the sensitivity $\frac{dG}{dp}$ is sought (i.e a in the above example).



• Also why {"ParametricSensitivity" -> "AdjointSensitivity"} fails in the above example?

For some math reference check here.

BR

Answer

OK,

eqns = {y''[t] + y[t] == 3 a Sin[y[t]], y[0] == y'[0] == 1};
pfun = ParametricNDSolveValue[eqns, 
   Integrate[y[s] - a s, {s, 0, 5}], {t, 0, 5}, {a}
   (*,Method\[Rule]{"ParametricSensitivity"\[Rule]\"ForwardSensitivity"}*)];

pfun[1.5]

does not return 'meaningless' stuff but a symbolic integral. And in fact you have requested it to return an Integrate. So to evaluate it just call N.

N[pfun[1.5]]
(* -7.86673 *)

And the same holds for the derivatives. Now, you can not request an NIntegrate in NDSolve, since that would give messages during the function call since then the input to NIntegrate were then symbolic. And similar for the derivative. As to why the second derivative self evaluates, I do not know. If you feel passionate about it, you might report a bug to WRI.

For a sanity check we can use:

Block[{a = 1.5}, 
 eqns = {y''[t] + y[t] == 3 a Sin[y[t]], y[0] == y'[0] == 1}; 

 nds = NDSolveValue[eqns, y[t], {t, 0, 5}];
 NIntegrate[nds - a t, {t, 0, 5}]
 ]
(* -7.86673 *)

Using the

, Method -> {"ParametricSensitivity" -> "AdjointSensitivity"}

gives a warning message:

(*

ParametricNDSolveValue::adjsens: "The adjoint sensitivity method cannot be used for the output function {t,0,5}. It can only be used for output functions that are at a particular time or are a definite integral over time."
*)

So only something like this is possible:

eqns = {y''[t] + y[t] == 3 a Sin[y[t]], y[0] == y'[0] == 1};
pfun = ParametricNDSolveValue[eqns, 
   Integrate[y[1.5], {s, 0, 5}], {t, 0, 5}, {a}
   , Method -> {"ParametricSensitivity" -> "AdjointSensitivity"}];
pfun[1.5]

I think ParametricNDSolve is quite a useful function and produces everything else then meaningless output.