Skip to main content
3 of 9
added 135 characters in body; edited title; edited tags
m_goldberg
  • 108.1k
  • 16
  • 104
  • 259

Problem finding inverse of numerical function using NDSolve and FindRoot

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.

My first attempt of inverting was by using interpolation of the "inversed data", as in this answer: http://mathematica.stackexchange.com/a/7195/6138 but since the grid becomes non-regular, I get many problems of precision and I don't know if I can evaluate this function safely, I get weird plots.

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, but I don't understand where it comes from. At the end there seems to be a bug when I repeat the process with previous functions.

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?NumericQ, k] = Exp[NIntegrate[funcA[et,k], {et, tini,t}]]

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

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]

fInverse[1.0]

Now the bug: If I know evaluate something that should be no problem at all I get the same 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}}. >>

This shouldn't happen since the function funcB is clearly well defined there.

If I restart Mathematica (or quit the Kernel) and don't evaluate these inverse functions, I can without problem perform this last evaluation of ParallelTable.

It seems that this is a bug, because it happens only after this error of the inverse functions. With (serial) Table, nothing bad seems to happen.

Edit: After checking the corresponding output of ParallelTable and Table, they are equal, but with ParallelTable I get the strange error message.

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

Santiago
  • 1.2k
  • 8
  • 16