MathGroup Archive 2005

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: hints to speed up numerical integration?

  • To: mathgroup at smc.vnet.net
  • Subject: [mg60159] Re: [mg60114] hints to speed up numerical integration?
  • From: Bob Hanlon <hanlonr at cox.net>
  • Date: Mon, 5 Sep 2005 02:27:51 -0400 (EDT)
  • Reply-to: hanlonr at cox.net
  • Sender: owner-wri-mathgroup at wolfram.com

The second method below is only marginally faster (approx. 16%), but it may 
help.

datanw={{0.0400,36.38,0.19},{0.050,36.84,0.21},
      {0.0307,35.90,0.20},{0.0560,37.31,0.18},
      {0.0331,35.54,0.20},{0.0141,34.13,0.25},
      {0.0460,36.35,0.21},{0.0265,35.64,0.20},
      {0.101,38.73,0.20},{0.075,37.77,0.19}};

Clear[f,ff,rr,chi2f2,chi2mar];
ndat=10;
z=datanw[[#1,1]]&;
mi=datanw[[#1,2]]&;
smi=datanw[[#1,3]]&;
SetAttributes[smi,Listable];
f[x_,om_,w_]:=
    1/Sqrt[om*(1+x)^3+(1-om)*(1+x)^(3*(1+w))];
rr[1,om_,w_]:=
    rr[1,om,w]=NIntegrate[f[x,om,w],{x,0,z[1]},Compiled->False];
rr[i_?IntegerQ,om_,w_]:=rr[i,om,w]=
      rr[i-1,om,w]+NIntegrate[f[x,om,w],{x,z[i-1],z[i]},Compiled->False];
ff[i_?IntegerQ,om_,w_]:=5*Log[10,rr[i,om,w]*(1+z[i])];
ci=Sum[1/smi[i]^2,{i,1,ndat}];
chi2f2[(om_)?NumericQ,(w_)?NumericQ]:=
    Module[{vec=((mi[#1]-ff[#1,om,w])/smi[#1]&)/@Range[ndat]},
      vec.vec-Total[vec/smi[Range[ndat]]]^2/ci];
chi2mar[w_]:=
    -2*Log[
        NIntegrate[Exp[-chi2f2[0.23+omvar,w]/2],
          {omvar,-0.07,0.07},
          Compiled->False]];

m1={chi2mar[0.1]//Timing,chi2mar[0.3]//Timing};

Clear[rr,chi2f2,chi2mar];
ndatr=Range[Length[datanw]];
z=First/@datanw;
bnds=Partition[z,2,1];
mi=#[[2]]&/@datanw;
smi=Last/@datanw;
ci=Total[1/smi^2];
rr[om_,w_]:=FoldList[
      #1+NIntegrate[
            1/Sqrt[om*(1+x)^3+(1-om)*(1+x)^(3*(1+w))],
            {x,#2[[1]],#2[[2]]}]&,
      NIntegrate[
        1/Sqrt[om*(1+x)^3+(1-om)*(1+x)^(3*(1+w))],
        {x,0,z[[1]]}], bnds];
chi2f2[om_?NumericQ,w_?NumericQ]:=
    Module[
      {vec=(mi-(5*Log[10,rr[om,w]*(1+z)]))/smi},
      vec.vec-Total[vec/smi]^2/ci];
chi2mar[w_]:=
    -2*Log[
        NIntegrate[Exp[-chi2f2[0.23+omvar,w]/2],
          {omvar,-0.07,0.07}]];

{chi2mar[0.1]//Timing,chi2mar[0.3]//Timing}/m1;

{{0.8369460350513395, 1.}, {0.8332640830442495, 0.9999999999999999}}


Bob Hanlon

> 
> From: "Ruth Lazkoz" <ruth.lazkoz at ehu.es>
To: mathgroup at smc.vnet.net
> Date: 2005/09/03 Sat AM 02:06:01 EDT
> Subject: [mg60159] [mg60114] hints to speed up numerical integration?
> 
> 
> Dear colleagues,
> 
> I would like to speed up this code, because I have to minimize the 
> function chi2mar defined below. I have tried Compiled -> True in the 
> numerical integrals but it would only make things worse.
> 
> Thanks for your cooperation,
> 
> Ruth Lazkoz
> 
> In[1]:=
> datanw = {{0.0400, 36.38, 0.19}, {0.050, 36.84, 0.21}, {0.0307, 35.90, 
> 0.20},
> {0.0560, 37.31, 0.18}, {0.0331, 35.54, 0.20}, {0.0141, 34.13, 0.25}, 
> {0.0460,
> 36.35, 0.21}, {0.0265, 35.64, 0.20}, {0.101, 38.73, 0.20}, {0.075,
> 37.77, 0.19}};
> 
> In[2]:=
> ndat = 10;
> 
> In[3]:=
> z = datanw[[#1, 1]] & ; mi = datanw[[#1, 2]] & ; smi 
> datanw[[#1, 3]] & ; SetAttributes[smi, Listable];
> f[x_, om_, w_] := 1/Sqrt[om*(1 + x)^3 + ( 1 - om)*(1 + x)^(3*(1 + 
> w))];
> rr[1, om_, w_] :=
> rr[1, om, w] = NIntegrate[f[x, om, w], {x, 0, z [1]}, Compiled -> 
> False];
> rr[i_?IntegerQ, om_, w_] :=
> rr[i, om, w] =
> rr[i - 1, om, w] +
> NIntegrate[f[x, om, w], {x, z[i - 1], z[i]}, Compiled -> False];
> ff[i_?IntegerQ, om_, w_] := 5*Log[ 10, rr[i, om, w]*(1 + z[i])]; ci =
> Sum[1/smi[i]^2, {i, 1, ndat}];
> chi2f2[(om_)? NumericQ, (w_)?NumericQ] :=
> Module[{vec = ((mi[#1] - ff[#1, om, w])/ smi[#1] & ) /@ Range[ndat]},
> vec . vec - Total[vec/smi[Range[ndat]]]^2/ci]
> 
> In[4]:=
> chi2mar[w_] := -2*
> Log[NIntegrate[Exp[-chi2f2[0.23 + omvar, w]/2], {omvar, -0.07, 0.07},
> Compiled -> False]]
> 
> 
> In[5]:=
> chi2mar[0.1] // Timing
> Out[5]=
> {0.078 Second, 15.5573}
> 
> 
> In[6]:=
> chi2mar[0.3] // Timing
> Out[6]=
> {0.078 Second, 15.7998}
> 
> 


  • Prev by Date: Re: Re: Re: inconsistency with Inequality testing and Floor (long)
  • Next by Date: TableHeadings for 1x1 tables
  • Previous by thread: Re: Re: hints to speed up numerical integration?
  • Next by thread: Re: Does ContourPlot behave correctly?