Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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 http://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

share|improve this question
    
Try funcB[t_?NumericQ, k_] := funcB[t, k] = Exp[NIntegrate[funcA[et,k], {et, tini,t}]] –  m_goldberg Aug 15 at 16:14
    
That's what I already have, isn't it? –  Santi Aug 15 at 16:44
    
No, not in your post. You have funcB[t_?NumericQ, k_] := funcB[t?NumericQ, k] = Exp[NIntegrate[funcA[et, k], {et, tini, t}]]. The t?NumericQ is bad. –  m_goldberg Aug 15 at 16:53
    
Ok thanks! I'll try that. Maybe that causes the strange behaviour without warnings. –  Santi Aug 15 at 16:59
    
@m_goldberg I tried that and I still get the same result. Maybe I should rephrase the question to make it more accesible and understandable. The basic problem is that when I run this code without calling the inverse functions, I get no error messages at all. If I run the inverse function, I get the error message there and also when calling ParallelTable. Clearing the functions A and B doesn't help anymore then. –  Santi Aug 17 at 12:16

1 Answer 1

up vote 2 down vote accepted

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]
share|improve this answer
    
thanks for taking the time to explain carefully the problem, I would have never found this subtlety. Should I change the title of my answer and remove the word bug? –  Santi Aug 20 at 11:39
    
@Santi You're welcome. I suppose "bug" is misleading. Something like "problem" or "issue" might seem more appropriate. –  Michael E2 Aug 20 at 12:21

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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