numerical integration - Integrating a list of functions

For the purpose of this minimal example, let's say we have a list of functions, like this:

f[y_?NumericQ] := {NIntegrate[z*y, {z, 0, 1}], a y}

I want to integrate an expression involving f, say

NDSolve[{y'[t] == f[y[t]][[1]], y[0] == 1}, y[t], {t, 0, 5}]

Now, the problem is that this doesn't return the expected result, because f[y[t]][[1]] evaluates to y[t] inside the NDSolve.

How can this be done correctly?

Answer

You can include the part extraction as an argument of your function, perhaps as a SubValues definition:

ClearAll[f]


f[y_?NumericQ][part_] := {NIntegrate[z*y, {z, 0, 1}], a y}[[part]]

NDSolve[{y'[t] == f[y[t]][1], y[0] == 1}, y[t], {t, 0, 5}]

{{y[t]->InterpolatingFunction[{{0.,5.}},<>][t]}}

Or, inside the primary body as an optional argument:

ClearAll[f]


f[y_?NumericQ, part_: All] := {NIntegrate[z*y, {z, 0, 1}], a y}[[part]]

NDSolve[{y'[t] == f[y[t], 1], y[0] == 1}, y[t], {t, 0, 5}]

This second method returns both values by default:

f[3.6]

{1.8, 3.6 a}

---

An alternative that comes to mind is to use a custom Part function that won't trigger when it should not, e.g.:

ClearAll[f]
listPart[x_List, part__] := x[[part]]

f[y_?NumericQ] := {NIntegrate[z*y, {z, 0, 1}], a y}

NDSolve[{y'[t] == f[y[t]] ~listPart~ 1, y[0] == 1}, y[t], {t, 0, 5}]