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?

-
I believe this is the related quesiton I was thinking of and it's not a duplicate, but perhaps interesting nevertheless: (11772) –  Mr.Wizard Mar 19 '13 at 12:39
Possible duplicate: 14645 –  Mr.Wizard Oct 20 '13 at 5:50

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}]

-

Mathematica 10 introduces something like listPart, with additional functionality, in Indexed:

Indexed can be used to indicate components of symbolic vectors, matrices, tensors, etc.

When expr is a list, Indexed[expr,i] gives expr[[i]].

When expr is a list, Indexed[expr,{i,j,...}] gives Indexed[expr[[i]],{j,...}].

Indexed can be used instead of Part for single-element extraction (at any level).

Applied to the case at hand:

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

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

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

-