Issue in ParallelTable after evaluating another function using NDSolve and FindRoot. What is wrong with this inverse?

Question

I am trying to find the inverse of a function which is defined through NDSolve and NIntegrate.

The question is pretty similar to How to invert an integral equation. Just that in my case $g_1(x)=1$ and $g_2(x)$ is the solution of an NDSolve.

Following the first answer https://mathematica.stackexchange.com/a/4604/6138, I tried to implement this using NDSolve and FindRoot, but I get an error message and I don't understand where it comes from.

The bug appears in the fact that if I don't evaluate the inverse function, there is no problem with ParallelTable or Table. After evaluating the inverse function, ParallelTable (but not Table) shows me the same error message. I would also like to know why FindRoot doesn't work in this case.

Below is my minimal working code and some comments:

Needs["NumericalCalculus`"]
kvalues = {0.01,0.1,0.5,1.0,2.0};
tini = -Log[100]; 
tfin = 0;
timeGrid = Range[tini,tfin,(tfin-tini)/10];
fSpace[min_, max_, steps_, f_: Log] :=
  InverseFunction[f] /@ Range[f @ min, f @ max, (f @ max - f @ min)/(steps -1)]
kin = 0.00072427;
kfin = 2.159;
klogGrid = fSpace[kin, kfin, 10];

funcA = 
  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, kin, kfin}]

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

funcD[td_, kd_] := ND[Log[funcB[tt, kd]], tt, td]

Now a simple evaluation with ParallelTable:

BTable = ParallelTable[{{tt, kk}, funcB[tt, kk]}, {tt, timeGrid}, {kk,
klogGrid}]

The output is (shortened):

{{{-Log[100], 0.00072427}, 1.}, {{-Log[100], 0.00176173}, 1.}, {{-Log[100], 0.00428528}, 1.}, {{-Log[100], 0.0104236}, 1.}, {{-Log[100], 0.0253546}, 1.}, {{-Log[100], 0.0616732}, 1.}, {{-Log[100], 0.150015}, 1.}, {{-Log[100], 0.3649}, 1.}, {{-Log[100], 0.887592}, 1.}, {{-Log[100], 2.159}, 1.}...

Now if I try:

inverseN = FindRoot[funcD[tt, kk] == funcA[tt, kk], {tt, -1}]

I get:

NIntegrate::inumr: The integrand InterpolatingFunction[{{-4.60517,0.},{0.00072427,2.159}},{4,5,1,{54,25},{4,5},0,0,0,0,Automatic,{},{},False},{<<1>>},{Developer`PackedArrayForm,{<<1>>},{1.,-5.10517,1.,-5.09284,1.,-5.05626,1.,-4.99541,1.,-4.9103,1.,-4.80093,1.,-4.6673,1.,-4.5094,1.,-4.32725,1.,-4.12083,1.,-3.89015,1.,<<5>>,1.,-2.72483,1.,-2.37285,1.,-1.99661,1.,-1.5961,1.,-1.17133,1.,-0.722307,1.,-0.249018,1.,0.248533,1.,0.770344,1.,1.31642,1.,1.88675,<<2650>>}},{Automatic,Automatic}][<<1>>] has evaluated to non-numerical values for all sampling points in the region with boundaries {{-Log[100],-4.60186}}. >>

The same thing happens if I use the inverse function method suggested in the Mathematica documentation of FindRoot:

inv[f_, s_] := Function[{t}, s /. FindRoot[f - t, {s, -2}]]
fInverse = inv[funcB[tt, kk], tt]

Now the bug: If I know evaluate again ParallelTable, I get the same above error message:

BTable = ParallelTable[{{tt, kk}, funcB[tt, kk]}, {tt, timeGrid}, {kk,
klogGrid}]

...

NIntegrate::inumr: The integrand InterpolatingFunction[{{-4.60517,0.},{0.00072427,2.159}},....Automatic}][<<1>>] has evaluated to non-numerical values for all sampling points in the region with boundaries {{-Log[100],-4.60186}}. >>

But after the error messages I get the correct output from above.

Nevertheless, with Table nothing bad happens:

Table[{{tt, kk}, funcB[tt, kk]}, {tt, timeGrid}, {kk, klogGrid}]

{{{{-Log[100], 0.00072427}, 1.}, {{-Log[100], 0.00176173}, 1.}, {{-Log[100], 0.00428528}, 1.}, {{-Log[100], 0.0104236}, 1.}, {{-Log[100], 0.0253546}, 1.}, {{-Log[100], 0.0616732}, 1.},....

I hope you can reproduce this bug and by the way help me with finding the inverse function.

I tried this in Mathematica 10.0 and 9.0

Answer 1

The reasons for the change in the behavior of ParallelTable are subtle. The main source of the problem is that in funcB, the argument k_ is not protected with ?NumericQ like this:

funcB[t_?NumericQ, k_?NumericQ] := (* a solution *)
 funcB[t, k] = Exp[NIntegrate[funcA[et, k], {et, tini, t}]]

But more on that later. The problem does not appear in the first ParallelTable, but the trap is set by executing

inverseN = FindRoot[funcD[tt, kk] == funcA[tt, kk], {tt, -1}]

I will explain how the trap is set later. First observe that this command naturally produces error messages because tt is numeric but kk is not. Thus the NIntegrate is ultimately invoked with a nonnumeric integrand when funcD calls funcB. If kd_ or k_ were protected with ?NumericQ, only FindRoot would complain and evaluation would never reach NIntegrate.

The key to the trap is that evaluation reaches funcB, which memoizes its evaluation. Since the argument k_ is the symbol kk, NIntegrate complains and returns unevaluated. This unevaluated NIntegrate is saved in funcB[tt, kk].

?funcB

Global`funcB

(*
  funcB[-1.,kk] = E^NIntegrate[funcA[et, kk], {et, tini, -1.}]
  ...
  funcB[t_?NumericQ, k_] := funcB[t, k] = Exp[NIntegrate[...
*)

On the face of it, this should be no problem since ParallelTable substitutes numeric, not symbolic, values for kk. Ah, but the surprise is that is not where the problem arises. The problem arises from distributing the definitions! ParallelTable distributes definitions before it begins tabulating.

If the definitions above are evaluated on the parallel kernels, the NIntegrate will be evaluated, which will give errors since kk is a symbol. The following will show it.

(* Evaluate the OP's definitions *)
inverseN = FindRoot[funcD[tt, kk] == funcA[tt, kk], {tt, -1}]
(* abort the above after the NIntegrate::inumr messages *)

DistributeDefinitions[funcB]

NIntegrate::inumr: The integrand <<1>> has evaluated to non-numerical values for all sampling points in the region with boundaries {{-Log[100], -4.60186}}. >>
...

For the sake of curiosity, the following definition works, which shows the difference between Set and SetDelayed:

funcB[t_?NumericQ, k_] := 
  funcB[t, k] := Evaluate@Exp[NIntegrate[funcA[et, k], {et, tini, t}]]

The definition given at the beginning is the best one to use. It would be good to do something similar with funcD and any other function which does not make sense on nonnumeric inputs.

funcD[td_?NumericQ, kd_?NumericQ] := ND[Log[funcB[tt, kd]], tt, td]