Problem when defining function through NIntegrate and NDSolve and Interpolation - Bug?

Question

More than a single question, I have some doubts about the output of certain functions when defined through the result of other calculations. I am an active user of Mathematica, but maybe I haven't read deeply enough about Attributes or related things. Here is the following minimal example:

{tini, tfin} = {-Log[100], 0};
firstFuncK=NDSolve[{D[f[t,k],t]+f[t,k]^2+(1-t)*f[t,k]==3/2*
(1+k^2),f[tini,k]==1},f,{t,tini,tfin},{k,0.001,10}]

secondFuncK[t_?NumberQ, k_] := 
Exp[NIntegrate[f[et, k] /. firstFuncK, {et, tini, t}]]

thirdFuncK[t_,k_]:=Log[N[secondFuncK[t,k]/secondFuncK[-Log[100],k]]][[1]]

Here comes the first question, why exactly do I need the [[1]] or First and the N. Without them I get:

thirdFuncK[-1,0.01]
{Log[4.11782/secondFuncK[-Log[100],0.01]]}

It is weird that Mathematica doesn't output a number, since secondFuncK is already numerical.

And I don't understand why I get a List in this case:

thirdFuncLin[et_]=thirdFuncK[et,0.01];

Output:

thirdFuncLin[-1]
{1.41532}

Plot[thirdFuncLin[tt],{tt,-4,0}]

But the real problem is when I use this function to interpolate:

thirdInterp=Interpolation[Table[{et,thirdFuncLin[et][[1]]},{et,-Log[100],0,0.01}]]

I get a completely different result:

Plot[thirdInterp[tt], {tt, -4, 0}]

I know there are related questions of simple mistakes people do, but I haven't found anything that helps me really to understand the core of this problem.

Thanks for any suggestions also on style or optimization.

Answer 1

Update

The problem is subtler than my first analysis revealed. There is indeed a problem with the variable et in NIntegrate not being properly blocked. Part of the problem has to do with the extra braces in firstFuncK which has the form

{{f -> InterpolatingFunction[<>]}}

Somehow that leads to an evaluation of et in the integrand f[et, k] /. firstFuncK in the definition of secondFuncK. This may be observed in the following minimal example, in which x plays the role of the OP's et:

g0 = {{g -> (#^2 &)}};
fn[x1_?NumericQ] := NIntegrate[g[x] /. g0, {x, 1, x1}]

Table[fn[x], {x, 4}]
fn /@ Range@4
(*
  {{0.}, {4.}, {18.}, {48.}}
  {{0.}, {2.33333}, {8.66667}, {21.}}
*)

The outputs are the same if we define g0 = First @ {{g -> (#^2 &)}}. This happens V8/9/10.

A consequence is that the bug could be avoided with using First on NDSolve:

{tini, tfin} = {-Log[100], 0};
firstFuncK = First @ NDSolve[{
   D[f[t, k], t] + f[t, k]^2 + (1 - t)*f[t, k] == 3/2*(1 + k^2), f[tini, k] == 1}, 
   f, {t, tini, tfin}, {k, 0.001, 10}]
secondFuncK[t_?NumberQ, k_] := Exp[NIntegrate[f[et, k] /. firstFuncK, {et, tini, t}]]
thirdFuncK[t_, k_] := Log[N[secondFuncK[t, k]/secondFuncK[-Log[100], k]]]
thirdFuncLin[et_] = thirdFuncK[et, 0.01]
thirdInterp = Interpolation[Table[{et, thirdFuncLin[et]}, {et, -Log[100], 0, 0.01}]]

The problem and a fix

The problem arises because of the use of et as a variable in NIntegrate and as the iterator symbol in Table. Table effectively uses Block to set the value of et. This interferes with NIntegrate. You can either use a different variable or protect et by using Block like this:

secondFuncK[t_?NumericQ, k_] := 
 Block[{et}, Exp[NIntegrate[f[et, k] /. firstFuncK, {et, tini, t}]]]

With this definition, we get

thirdInterp = (* needs reevaluation *)
 Interpolation[Table[{et, thirdFuncLin[et][[1]]}, {et, -Log[100], 0, 0.01}]];

Plot[thirdFuncLin[tt], {tt, -4, 0}]

Original analysis of the problem

Keep the OP's definitions for testing.

Below we see that blocking et makes a difference and yields the same result as Table:

Block[{et = -0.`}, thirdFuncLin[et]]
(* {3.34341} *)

Block[{x = -0.`}, thirdFuncLin[x]]
(* {2.03644} *)

Table[thirdFuncLin[et], {et, {0.}}]
(* {{3.34341}} *)

Where the 3.34341 comes from. The integrand f[et, 0.01] of NIntegrate evaluates to f[0, 0.01], which is then integrated. Since this is constant, we can check by hand as follows:

{tini, tfin} = {-Log[100], 0}; (* omitted by OP from Q *)
f[et, 0.01]*(0 - tini) /. firstFuncK /. et -> 0.
(* {3.34341} *)

Comment: Bug?

This appears to be a bug in V10/V9/V8. The docs forNIntegrate state

NIntegrate has attribute HoldAll and effectively uses Block to localize variables.

This does not happen here. Confirmation would be appreciated.

Edit - The minor problems

I was so focused on the potential bug that I forgot about the OP's other issues.

The solution returned by NDSolve has the form of a List of a solutions, each solution being itself a List of substitution rules, one Rule for each variable. In the OP's case, it has the form

{{f -> InterpolatingFunction[<>]}}

When this is used with ReplaceAll (./ -- see the last "Basic Example" in the documentation), you get a list:

f /. {{f -> InterpolatingFunction[<>]}}    
(* {InterpolatingFunction[<>]} *)

It would be nicer to get just the function without the list. To do that, knowing there is only one solution, one define firstFuncK like this with First,

firstFuncK = First @ NDSolve[<etc>]

or like this (V9+) withNDSolveValue,

firstFuncK = NDSolveValue[<etc>]

That removes the need for [[1]] in later definitions and makes thirdFuncLin evaluate to a number instead of a List; this is because using firstFuncK with ReplaceAll results in a function, not a list:

f /. firstFuncK
(* InterpolatingFunction[<>] *)

Next, N needs to be used because the OP used the PatternTest _?NumberQ instead of _?NumericQ. The difference is that -Log[100] is a numeric expression, but not technically a number, which is one of Integer, Rational, Real, or Complex.

Answer 2

Two small changes make this work in a much nicer way:

First, the use of NDSolveValue instead of NDSolve gets rid of this rule replacement monkey business.

{tini, tfin} = {-Log[100], 0};
firstFuncK = 
 NDSolveValue[{D[f[t, k], t] + f[t, k]^2 + (1 - t)*f[t, k] == 
    3/2*(1 + k^2), f[tini, k] == 1}, 
  f, {t, tini, tfin}, {k, 0.001, 10}]
(* InterpolatingFunction[....] *)

This then makes the second and subsequent functions much cleaner:

secondFuncK[t_?NumericQ, k_] := 
 Exp[NIntegrate[firstFuncK[et, k], {et, tini, t}]]

thirdFuncK[t_, k_] := 
 Log[secondFuncK[t, k]/secondFuncK[-Log[100], k]]

thirdFuncLin[et_] = thirdFuncK[et, 0.01];

thirdFuncLin[-1]
(*1.415323046633079`*)

Plot[thirdFuncLin[tt], {tt, -4, 0}]

Note, the additional change: I replaced the _?NumberQ with a _?NumericQ since Log[100] is not a number (it's a symbolic expression) but it is a numeric quantity. That's also the reason why the N was needed; it forces the Log[100] to be _?NumberQ.