Re: Problem with NIntegrating a compiled function
- To: mathgroup at smc.vnet.net
- Subject: [mg121094] Re: Problem with NIntegrating a compiled function
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Sat, 27 Aug 2011 08:16:35 -0400 (EDT)
- Delivered-to: l-mathgroup@mail-archive0.wolfram.com
- References: <201108260924.FAA04553@smc.vnet.net>
On 08/26/2011 04:24 AM, Eric Michielssen wrote:
> Hi,
>
> I am running into unexpected problems using NIntegrate on a compiled
> function. In the code below, I NIntegrate the functions "myint" and
> "cmyint"; the latter is a compiled version of the former -- both routines
> are listed at the bottom of this email.
>
> ClearAll[dt,r,c,i];
> dt=0.1;
> c=1.;
> r=0.95;
> i=11;
> Print["Integrate using noncompiled modules"];
> NIntegrate[myint[r,c,dt,i],{r,0.95,1.05},EvaluationMonitor:>Print[r,"
> ",myint[r,c,dt,i]]]
> Print["Integrate using compiled modules"];
> NIntegrate[cmyint[r,c,dt,i],{r,0.95,1.05},EvaluationMonitor:>Print[r "
> ",cmyint[r,c,dt,i]," ",myint[r,c,dt,i] ," ",myint[r,c,dt,i]
> -cmyint[r,c,dt,i]]]
>
> The noncompiled NIntegrate yields a nonzero answer. The compiled one yields
> zero, even though EvaluationMonitor indicates the integrands evaluate to the
> same numerical value for all values of the integration variable. Changing
> the NIntegrate options (PrecisionGoal, MinRecursion etc) does not change the
> behavior.
>
> During evaluation of In[2288]:= Integrate using noncompiled modules [[I
> print integration variable and myint]]
>
> During evaluation of In[2288]:= 0.950796 0.000159148
> During evaluation of In[2288]:= 0.954691 0.000169447
> During evaluation of In[2288]:= 0.962292 0.000191239
> During evaluation of In[2288]:= 0.973077 0.000225252
> During evaluation of In[2288]:= 0.986018 0.000269762
> During evaluation of In[2288]:= 1. 0.00032485
> During evaluation of In[2288]:= 1.01398 0.000381631
> During evaluation of In[2288]:= 1.02692 0.000407476
> During evaluation of In[2288]:= 1.03771 0.000412319
> During evaluation of In[2288]:= 1.04531 0.00040786
> During evaluation of In[2288]:= 1.0492 0.000403319
> During evaluation of In[2288]:= 0.950398 0.000158134
> During evaluation of In[2288]:= 0.952346 0.000163168
> During evaluation of In[2288]:= 0.956146 0.000173454
> During evaluation of In[2288]:= 0.961538 0.000188989
> During evaluation of In[2288]:= 0.968009 0.000208878
> During evaluation of In[2288]:= 0.975 0.000231631
> During evaluation of In[2288]:= 0.981991 0.000255507
> During evaluation of In[2288]:= 0.988462 0.000278613
> During evaluation of In[2288]:= 0.993854 0.000298889
> During evaluation of In[2288]:= 0.997654 0.000314236
> During evaluation of In[2288]:= 0.999602 0.000322896
> During evaluation of In[2288]:= 1.0004 0.000326926
> During evaluation of In[2288]:= 1.00235 0.000336681
> During evaluation of In[2288]:= 1.00615 0.000353817
> During evaluation of In[2288]:= 1.01154 0.000374006
> During evaluation of In[2288]:= 1.01801 0.000392226
> During evaluation of In[2288]:= 1.025 0.000405096
> During evaluation of In[2288]:= 1.03199 0.00041149
> During evaluation of In[2288]:= 1.03846 0.000412151
> During evaluation of In[2288]:= 1.04385 0.000409176
> During evaluation of In[2288]:= 1.04765 0.000405299
> During evaluation of In[2288]:= 1.0496 0.000402774
> During evaluation of In[2288]:= 0.950199 0.000157629
> During evaluation of In[2288]:= 0.951173 0.000160117
> During evaluation of In[2288]:= 0.953073 0.000165091
> During evaluation of In[2288]:= 0.955769 0.000172408
> During evaluation of In[2288]:= 0.959005 0.000181561
> During evaluation of In[2288]:= 0.9625 0.000191865
> During evaluation of In[2288]:= 0.965995 0.000202557
> During evaluation of In[2288]:= 0.969231 0.000212766
> During evaluation of In[2288]:= 0.971927 0.000221482
> During evaluation of In[2288]:= 0.973827 0.000227732
> During evaluation of In[2288]:= 0.974801 0.000230967
> During evaluation of In[2288]:= 0.975199 0.000232295
> During evaluation of In[2288]:= 0.976173 0.000235561
> During evaluation of In[2288]:= 0.978073 0.000241993
> During evaluation of In[2288]:= 0.980769 0.000251257
> During evaluation of In[2288]:= 0.984005 0.000262586
> During evaluation of In[2288]:= 0.9875 0.000275109
> During evaluation of In[2288]:= 0.990995 0.00028799
> During evaluation of In[2288]:= 0.994231 0.000300358
> During evaluation of In[2288]:= 0.996927 0.000311187
> During evaluation of In[2288]:= 0.998827 0.000319334
> During evaluation of In[2288]:= 0.999801 0.000323855
> During evaluation of In[2288]:= 0.975099 0.000231963
> During evaluation of In[2288]:= 0.975586 0.000233592
> During evaluation of In[2288]:= 0.976536 0.000236786
> During evaluation of In[2288]:= 0.977885 0.000241352
> During evaluation of In[2288]:= 0.979502 0.000246885
> During evaluation of In[2288]:= 0.98125 0.000252926
> During evaluation of In[2288]:= 0.982998 0.000259035
> During evaluation of In[2288]:= 0.984615 0.000264752
> During evaluation of In[2288]:= 0.985964 0.000269565
> During evaluation of In[2288]:= 0.986914 0.000272986
> During evaluation of In[2288]:= 0.987401 0.000274748
> During evaluation of In[2288]:= 0.987599 0.00027547
> During evaluation of In[2288]:= 0.988086 0.000277242
> During evaluation of In[2288]:= 0.989036 0.000280721
> During evaluation of In[2288]:= 0.990385 0.000285708
> During evaluation of In[2288]:= 0.992002 0.000291784
> During evaluation of In[2288]:= 0.99375 0.000298485
> During evaluation of In[2288]:= 0.995498 0.000305369
> During evaluation of In[2288]:= 0.997115 0.00031197
> During evaluation of In[2288]:= 0.998464 0.000317724
> During evaluation of In[2288]:= 0.999414 0.000322009
> During evaluation of In[2288]:= 0.999901 0.000324346
> Out[2294]= 0.0000312806
>
> During evaluation of In[2288]:= Integrate using compiled modules [[I print
> integration variable, myint, cmyint, and the difference between the latter
> two]]
> During evaluation of In[2288]:= CompiledFunction::cfsa: Argument r at
> position 1 should be a machine-size real number.>>
> During evaluation of In[2288]:= 0.950796 0.000159148 0.000159148
> -1.29508*10^-16
> During evaluation of In[2288]:= 0.954691 0.000169447 0.000169447
> 2.46385*10^-17
> During evaluation of In[2288]:= 0.962292 0.000191239 0.000191239
> -3.69713*10^-17
> During evaluation of In[2288]:= 0.973077 0.000225252 0.000225252
> -2.46927*10^-17
> During evaluation of In[2288]:= 0.986018 0.000269762 0.000269762
> 6.17995*10^-18
> During evaluation of In[2288]:= 1. 0.00032485 0.00032485 6.17995*10^-18
> During evaluation of In[2288]:= 1.01398 0.000381631 0.000381631
> 6.17995*10^-18
> During evaluation of In[2288]:= 1.02692 0.000407476 0.000407476
> 3.08456*10^-17
> During evaluation of In[2288]:= 1.03771 0.000412319 0.000412319
> -1.11022*10^-16
> During evaluation of In[2288]:= 1.04531 0.00040786 0.00040786
> -7.39968*10^-17
> During evaluation of In[2288]:= 1.0492 0.000403319 0.000403319
> 4.93312*10^-17
> During evaluation of In[2288]:= 0.950398 0.000158134 0.000158134
> -4.31784*10^-17
> During evaluation of In[2288]:= 0.952346 0.000163168 0.000163168
> -6.1664*10^-17
> During evaluation of In[2288]:= 0.956146 0.000173454 0.000173454
> 6.17995*10^-18
> During evaluation of In[2288]:= 0.961538 0.000188989 0.000188989
> 2.46927*10^-17
> During evaluation of In[2288]:= 0.968009 0.000208878 0.000208878
> -1.23599*10^-17
> During evaluation of In[2288]:= 0.975 0.000231631 0.000231631
> 6.20706*10^-18
> During evaluation of In[2288]:= 0.981991 0.000255507 0.000255507
> -8.63567*10^-17
> During evaluation of In[2288]:= 0.988462 0.000278613 0.000278613
> -6.16911*10^-17
> During evaluation of In[2288]:= 0.993854 0.000298889 0.000298889
> 5.55654*10^-17
> During evaluation of In[2288]:= 0.997654 0.000314236 0.000314236
> -7.4051*10^-17
> During evaluation of In[2288]:= 0.999602 0.000322896 0.000322896
> -8.64109*10^-17
> During evaluation of In[2288]:= 1.0004 0.000326926 0.000326926
> -8.63025*10^-17
> During evaluation of In[2288]:= 1.00235 0.000336681 0.000336681
> -7.39968*10^-17
> During evaluation of In[2288]:= 1.00615 0.000353817 0.000353817
> -1.84856*10^-17
> During evaluation of In[2288]:= 1.01154 0.000374006 0.000374006
> -8.01225*10^-17
> During evaluation of In[2288]:= 1.01801 0.000392226 0.000392226
> 6.17995*10^-18
> During evaluation of In[2288]:= 1.025 0.000405096 0.000405096
> 1.04951*10^-16
> During evaluation of In[2288]:= 1.03199 0.00041149 0.00041149
> -3.70255*10^-17
> During evaluation of In[2288]:= 1.03846 0.000412151 0.000412151
> -1.84856*10^-17
> During evaluation of In[2288]:= 1.04385 0.000409176 0.000409176
> -1.23382*10^-16
> During evaluation of In[2288]:= 1.04765 0.000405299 0.000405299
> 6.16369*10^-17
> During evaluation of In[2288]:= 1.0496 0.000402774 0.000402774
> 7.39968*10^-17
>
> During evaluation of In[2288]:= NIntegrate::izero: Integral and error
> estimates are 0 on all integration subregions. Try increasing the value of
> the MinRecursion option. If value of integral may be 0, specify a finite
> value for the AccuracyGoal option.>>
>
> Out[2296]= 0.
>
> The routines for myint and cmyint are below -- they both call other
> noncompiled and compiled routines. myintdet and cmyintdet.
>
> ClearAll[myintdet,cmyintdet,myint,cmyint];
>
> myintdet[t1_,t2_,a_,dt_,to_,j_?IntegerQ]:=Module[
> {at2,at1,lat},
> at2=Sqrt[-1. a^2+t2^2] ; at1= Sqrt[-1. a^2+t1^2];lat=
> Log[(at2+t2)/(at1+t1)];
> Switch[j,
> 1,0.05555555555555555` (2.` a^2 (at1-at2)+at1 t1^2-at2 t2^2+9.` (a^2 lat-at1
> t1+at2 t2) (dt+0.5` to)+3.` lat (dt+to) (2.` dt+to) (3.` dt+to)+3.`
> (at1-at2) (11.` dt^2+12.` dt to+3.` to^2)),
> 2,0.16666666666666666` (2.` a^2 (-at1+at2)-at1 t1^2+at2 t2^2+3.` lat (dt-to)
> (dt+to) (2.` dt+to)+1.5` (-1.` a^2 lat+at1 t1-at2 t2) (2.` dt+3.` to)+3.`
> (at1-at2) (dt^2-4.` dt to-3.` to^2)),
> 3,0.16666666666666666` (2.` a^2 (at1-at2)+at1 t1^2-at2 t2^2+3.` lat (dt-to)
> (2.` dt-to) (dt+to)+3.` (a^2 lat-at1 t1+at2 t2) (-1.` dt+1.5` to)+3.`
> (-at1+at2) (dt^2+4.` dt to-3.` to^2)),
> 4,0.05555555555555555` (2.` a^2 (-at1+at2)-at1 t1^2+at2 t2^2+3.` lat (dt-to)
> (2.` dt-to) (3.` dt-to)+9.` (a^2 lat-at1 t1+at2 t2) (dt-0.5` to)-3.`
> (at1-at2) (11.` dt^2-12.` dt to+3.` to^2)) ,
> _,
> 0.
> ]];
>
> cmyintdet=Compile[
> {{t1, _Real},{t2, _Real}, {a, _Real},{dt, _Real},{to,_Real},{j, _Integer}},
> Module[
> {at2,at1,lat},
> at2=Sqrt[-1. a^2+t2^2] ; at1= Sqrt[-1. a^2+t1^2];lat=
> Log[(at2+t2)/(at1+t1)];
> Switch[j,
> 1,0.05555555555555555` (2.` a^2 (at1-at2)+at1 t1^2-at2 t2^2+9.` (a^2 lat-at1
> t1+at2 t2) (dt+0.5` to)+3.` lat (dt+to) (2.` dt+to) (3.` dt+to)+3.`
> (at1-at2) (11.` dt^2+12.` dt to+3.` to^2)),
> 2,0.16666666666666666` (2.` a^2 (-at1+at2)-at1 t1^2+at2 t2^2+3.` lat (dt-to)
> (dt+to) (2.` dt+to)+1.5` (-1.` a^2 lat+at1 t1-at2 t2) (2.` dt+3.` to)+3.`
> (at1-at2) (dt^2-4.` dt to-3.` to^2)),
> 3,0.16666666666666666` (2.` a^2 (at1-at2)+at1 t1^2-at2 t2^2+3.` lat (dt-to)
> (2.` dt-to) (dt+to)+3.` (a^2 lat-at1 t1+at2 t2) (-1.` dt+1.5` to)+3.`
> (-at1+at2) (dt^2+4.` dt to-3.` to^2)),
> 4,0.05555555555555555` (2.` a^2 (-at1+at2)-at1 t1^2+at2 t2^2+3.` lat (dt-to)
> (2.` dt-to) (3.` dt-to)+9.` (a^2 lat-at1 t1+at2 t2) (dt-0.5` to)-3.`
> (at1-at2) (11.` dt^2-12.` dt to+3.` to^2)),
> _,
> 0.
> ]]
> ];
>
> myint[r_?NumericQ,c_?NumericQ,dt_?NumericQ,i_?NumericQ]:=Module[
> {a,to,limits,tab,t1,t2,j},
> a=r/c;
> to=i dt;
> limits={{to,dt+to},{-dt+to,to},{-2 dt+to,-dt+to},{-3 dt+to,-2 dt+to}};
> tab=Table[
> If[a>limits[[j,2]],
> 0.,
> {t1,t2}={Max[limits[[j,1]],a],limits[[j,2]]};
> myintdet[t1,t2,a,dt,to,j]
> ],
> {j,1,4}];
> Total[tab]
> ];
>
> cmyint=Compile[
> {{r,_Real},{c, _Real},{dt, _Real},{i, _Integer}},
> Module[
> {a,to,limits,tab,t1,t2,j,u=0.},
> a=r/c;
> to=i dt;
> limits={{to,dt+to},{-dt+to,to},{-2 dt+to,-dt+to},{-3 dt+to,-2 dt+to}};
> tab=Table[
> If[a>limits[[j,2]],
> u=0.,
> {t1,t2}={Max[limits[[j,1]],a],limits[[j,2]]};
> u=cmyintdet[t1,t2,a,dt,to,j] (*** FIX to cmyintdet ***)
> ];
> u,
> {j,1,4}];
> Return[Total[tab]]
> ]];
>
> Eric Michielssen
> Radiation Laboratory
> Department of Electrical Engineering and Computer Science
> University of Michigan
> 1301 Beal Avenue
> Ann Arbor, MI 48109
> Ph: (734) 647 1793 -- Fax: (734) 647 2106
This one had me puzzling a bit.
The problem is that the "function" created by NIntegrate is not what you
think it is. At first, either 'r is unitialized, or has a numeric value.
In the former case we get
In[74]:= r =.
In[75]:= cmyint[r, c, dt, i]
During evaluation of In[75]:= CompiledFunction::cfsa: Argument r at
position 1 should be a machine-size real number. >>
Out[75]= 0.
In the latter case we'll of course get a number back. Either way, that
result is exactly what NIntegrate preprocessing constructs as the
numeric function it must evaluate.
So your evaluation monitor is printing cmyint[r,c,dy,i] with 'r being
the current value of the variable of integration. But the actual
function that NIntegrate is evaluating is just some constant value.
Hence the NIntegrate result is zero.
Here is a viable workaround. Inhibit the integrand so that it cannot
evaluate unless r is numeric.
cmyintWrapper[r_?NumericQ, c_?NumericQ, dt_?NumericQ, i_?NumericQ] :=
cmyint[r, c, dt, i]
Now you can do
NIntegrate[cmyintWrapper[r, c, dt, i], {r, 0.95, 1.05}]
and get a result that is plausible.
Daniel Lichtblau
Wolfram Research
- References:
- Problem with NIntegrating a compiled function
- From: "Eric Michielssen" <emichiel@eecs.umich.edu>
- Problem with NIntegrating a compiled function