8
$\begingroup$

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

thirdFuncK

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

thirdInterp

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.

|improve this question
$\endgroup$
  • 1
    $\begingroup$ You should include the definitions of tini, tfin in your question -- I guessed below. $\endgroup$ – Michael E2 Aug 8 '14 at 16:49
  • $\begingroup$ you're right, done! $\endgroup$ – Santi Aug 10 '14 at 10:29
9
$\begingroup$

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

Mathematica graphics

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.

|improve this answer
$\endgroup$
  • $\begingroup$ Sorry for the late reply. Thanks for all your comments and suggestions. So, if I understand right, using Block should give me the right result or otherwise using First before NDSolve. But this shouldn't be necessary since NIntegrate should localize all variables and this is the bug right? $\endgroup$ – Santi Aug 10 '14 at 10:34
  • $\begingroup$ @Sanit That's right. (I may reorganize the answer when I have time. The two major edits have resulted in a somewhat disorganized answer and there are some ungrammatical or unclear bits.) $\endgroup$ – Michael E2 Aug 10 '14 at 12:29
  • $\begingroup$ Great, thanks! What I still don't understand from your minimal example is the difference between Table[fn[x], {x, 4}] and fn /@ Range@4. Why do they give different results? $\endgroup$ – Santi Aug 10 '14 at 13:24
  • $\begingroup$ Table effectively uses Block to set x equal to 1, 2, 3, 4 in turn; the bug allows the symbol x in the integrand of NIntegrate to be replaced by the numerical value before the integration. So NIntegrate integrates a constant; e.g. when x = 4, it calculates NIntegrate[{4^2}, {x, 1, 4}]. fn /@ Range@4 calculates fn[1], fn[2], etc., without changing x. fn[4] calculates NIntegrate[{x^2}, {x, 1, 4}]. The Table behavior can be imitated by Block[{x = #}, fn[x]] & /@ Range@4. $\endgroup$ – Michael E2 Aug 10 '14 at 14:34
4
$\begingroup$

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.

|improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.