Re: Advice on Compile
- To: mathgroup at smc.vnet.net
- Subject: [mg112503] Re: Advice on Compile
- From: "Niels R. Walet" <niels.walet at manchester.ac.uk>
- Date: Fri, 17 Sep 2010 06:44:02 -0400 (EDT)
Carl K. Woll wrote:
> On 9/15/2010 7:02 PM, Niels R. Walet wrote:
>> I am trying to write the following code (simplified, Mathematica 7.0.2)
>> setZu[subs_] :=
>> Block[{Zmint}, Zmint = Evaluate[EFRp[q] /. subs]; Clear[Zmn];
>> Zmn = Compile[{{k, _Real}},
>> NIntegrate[Evaluate[Zmint], {q, 0, k}]];]
>> I call
>> setZu[EFRp -> Function[{q}, q^2/(1 + q^2)^2]]
>> and than
>> Zmn[1.]
>> gives an error (because Zmint has not been expanded).
>> Of course I can move the assigment on Zmint within the compile, but it
>> leads to ugly code, and I am not sure how efficient the compile is (does
>> it do the substitution every time, or is it "compiled out"?). In real
>> life the integrals I calculate are extremely expensive, and give me the
>> coefficients for an ODE I try to solve, so extreme efficiency would be
>> appreciated: I want to compile the expression with the substitution
>> done, not the substitution process itself.
>> Any light on this issue would be appreciated.
>> Niels
>>
>>
>
> You need to use something like With to get Zmint evaluated:
>
> setZu[subs_] := With[{Zmint = EFRp[q] /. subs},
> Zmn = Compile[{{k, _Real}},
> NIntegrate[Evaluate[Zmint], {q, 0, k}]];]
>
> However, for this particular integral, it would be worth considering
> doing something like:
>
> nint = f /.
> First@NDSolve[{f'[q] == q^2/(1 + q^2)^2, f[0] == 0}, f, {q, 0, 10}]
>
> nint[2.3]//Timing
> Zmn[2.3]//Timing
>
> Also, care to share what the real integrals you are trying to do
> numerically?
>
> Carl Woll
> Wolfram Research
Interesting. I always struggle with the difference between With, Block
and Module...
I like the NDSolve trick; Unfortunately, in my real integrals the
parameter is in the integrand, not the boundary!
As far as the integrals go, a simpler taster (of a set of very long
expressions) is
s=0.5;
tsc[k_, q_] = (Erf[(-q + k)/(k s)] + Erf[(q + k)/(k s)])/(2 Erf[1/s]);
subs=EFRp -> Function[{q, k}, q^2/(2) + 1/(2 ) k^2 tssc[k, q] - 1/2;
NIntegrate[(q^2/(6 q EFRp[q,k]^2) (2 D[EFRp[q, k],q] + q D[EFRp[q,
k],q,q] )/.subs,{q,0,\[Infinity\]}]
At this point I no longer integrate to infinity; it is easier to
integrate to something like {0, k(1+20s)}.
Niels
--
Prof. Niels R. Walet Phone: +44(0)1613063693
School of Physics and Astronomy Fax: +44(0)1613064303
The University of Manchester Mobile: +44(0)7905438934
Manchester, M13 9PL, UK room 7.7, Schuster Building
email: Niels.Walet at manchester.ac.uk
web: http://www.theory.physics.manchester.ac.uk/~mccsnrw