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've got some problems with my code, and I try to make it faster. Some of you suggested me to split my problem, and I'm here...

I post the same function with and without Compile.

This is a support function and does not involve my question

r=8.314472
fug1 = Compile[{v, p, t, a, b},
Module[{y, z, vbv, vb, f1, f2, f3, f4, f},
y = b/(4 v);
z = (p v)/(r t);
vbv = Log[(v + b)/v];
vb = v + b;
f1 = (4.*y - 3.*y^2.)/(1 - y)^2.;
f2 = (4.*y - 2.*y^2.)/(1 - y)^3.;
f3 = (2.*vbv)/(r t*b)*a;
f4 = (vbv/b - 1./vb)/(r t)*a;
f = f1 + f2 - f3 + f4 - Log[z];
Exp[f]]]

g2 doesen't work and I don't know why

g2 = Compile[{p, t, a0, a1, a2, b0, b1, b2},
Module[{a, b, csd, vol, sol, vliquid, vvapor, fl, fv},
a = a0*Exp[a1*t + a2*t^2];
b = b0 + b1*t + b2*t^2;
csd = a/(r*t*(b + v)) - (-(b^3/(64.*v^3)) + b^2/(16.*v^2.0) + 
    b/(4.*v) + 1.)/(1 - b/(4*v))^3 + (p*v)/(r*t);
vol = NSolve[csd == 0. && v > 0., v, Reals] // Quiet;
sol = v /. vol;
vliquid = Min[sol];
vvapor = Max[sol];
fl = fug1[vliquid, p, t, a, b];
fv = fug1[vvapor, p, t, a, b];
Print[{t, p, vol, Abs[fl - fv]}];
Abs[fl - fv]],
RuntimeAttributes -> {Listable}]

This works without Compile!

g[p_?NumericQ, t_?NumericQ, a0_?NumericQ, a1_?NumericQ, a2_?NumericQ, 
  b0_?NumericQ, b1_?NumericQ, b2_?NumericQ] := 
Module[{a, b, csd, vol, sol, vliquid, vvapor, fl, fv},
a = a0*Exp[a1*t + a2*t^2];
b = b0 + b1*t + b2*t^2;
csd = a/(r*t*(b + v)) - (-(b^3/(64.*v^3)) + b^2/(16.*v^2.0) + 
    b/(4.*v) + 1.)/(1 - b/(4*v))^3 + (p*v)/(r*t);
vol = NSolve[csd == 0. && v > 0., v, Reals];
sol = v /. vol;
vliquid = Min[sol];
vvapor = Max[sol];
fl = fug1[vliquid, p, t, a, b];
fv = fug1[vvapor, p, t, a, b];
Print[{t, p, vol, Abs[fl - fv]}];
Abs[fl - fv]]; 

g works very well and is the same of g2!

FindRoot[g[p, 100, 500., -4.4627562855*10^-3, -2.7625748*10^-6, 
7.30402014*10^-2, -2.2222592*10^-4, 9.42486*10^-8], {p, 
34.376}] // Timing

g2 doesen't work

FindRoot[g2[p, 100, 500., -4.4627562855*10^-3, -2.7625748*10^-6, 
7.30402014*10^-2, -2.2222592*10^-4, 9.42486*10^-8], {p, 
34.376}] // Timing
share|improve this question
5  
There are many issues with with your functions. Variables are not localized, r is not defined in fug1, NSolve is not compilable and you probably mean to use g2 but are using function g1 in the FindRoot. –  asim Jan 27 '13 at 16:00
    
You're right. I've corrected it, I've localized my variables and I've inserted all your corrections, but it doesn't work again –  meriens Jan 27 '13 at 18:16
add comment

1 Answer 1

Let's start with solving your equation with Solve:

rv = r -> 8.314472;
csd = a/(r*t*(b + v)) - (-(b^3/(64*v^3)) + b^2/(16*v^2) + b/(4*v) + 
     1)/(1 - b/(4*v))^3 + (p*v)/(r*t);
rootsols = Solve[csd == 0, v] /. rv;
v /. Last[rootsols]

Mathematica graphics

I change fug1 slightly:

fug1 = Compile[{{v, _Real, 1}, p, t, a, b}, 
   Module[{y, z, vbv, vb, f1, f2, f3, f4, f, r},
    r = 8.314472;
    y = b/(4 v);
    z = (p v)/(r t);
    vbv = Log[(v + b)/v];
    vb = v + b;
    f1 = (4.*y - 3.*y^2.)/(1 - y)^2.;
    f2 = (4.*y - 2.*y^2.)/(1 - y)^3.;
    f3 = (2.*vbv)/(r t*b)*a;
    f4 = (vbv/b - 1./vb)/(r t)*a;
    f = f1 + f2 - f3 + f4 - Log[z];
    Exp[f]]];

Now, g2, using Root instead of NSolve:

    g2 = With[{fug1 = fug1}, 
   Compile[{p, t, a0, a1, a2, b0, b1, b2}, 
    Module[{a, b, csd, vol, vliquid, vvapor, fl, fv}, 
     a = a0*Exp[a1*t + a2*t^2];
     b = b0 + b1*t + b2*t^2;
     sol = {};
     i = 0;
     While[i < 5, i++;
      soli =
       Root[-a b^3 + 
          8.314472` b^4 t + (12 a b^2 - b^4 p - 
             24.943416` b^3 t) #1 + (-48 a b + 11 b^3 p - 
             166.28944` b^2 t) #1^2 + (64 a - 36 b^2 p - 
             665.15776` b t) #1^3 + (16 b p - 532.126208` t) #1^4 + 
          64 p #1^5 &, i];
      If[Im[soli] == 0, sol = Append[sol, soli]]
      ];

     vliquid = Min[sol];
     vvapor = Max[sol];
     {fl, fv} = fug1[{vliquid, vvapor}, p, t, a, b];

     Abs[fl - fv]], {{soli, _Complex}, {sol, _Real, 1}}, 
    CompilationOptions -> {"InlineCompiledFunctions" -> True},
    RuntimeAttributes -> {Listable}
    ]];

Testing with g:

g[p_?NumericQ, t_?NumericQ, a0_?NumericQ, a1_?NumericQ, a2_?NumericQ, 
   b0_?NumericQ, b1_?NumericQ, b2_?NumericQ] := 
  Module[{a, b, csd, vol, sol, vliquid, vvapor, fl, fv}, 
   a = a0*Exp[a1*t + a2*t^2];
   b = b0 + b1*t + b2*t^2;
   csd = a/(r*t*(b + v)) - (-(b^3/(64.*v^3)) + b^2/(16.*v^2.0) + 
        b/(4.*v) + 1.)/(1 - b/(4*v))^3 + (p*v)/(r*t);
   vol = NSolve[csd == 0. && v > 0., v, Reals];
   sol = v /. vol;
   vliquid = Min[sol];
   vvapor = Max[sol];
   {fl, fv} = fug1[{vliquid, vvapor}, p, t, a, b];
   (*Print[{t,p,vol,Abs[fl-fv]}];*)
   Abs[fl - fv]];
r = 8.314472; 
FindRoot[g[p, 100, 500., -4.4627562855*10^-3, -2.7625748*10^-6, 
    7.30402014*10^-2, -2.2222592*10^-4, 9.42486*10^-8], {p, 34.376}] //
   Quiet // Timing

gives

(* 
       {0.218401, {p -> 170.441}}
*)

and with the compile function (need a simple extra function g3 here):

g3[args__?NumberQ] := g2[args];
FindRoot[g3[p, 100, 500., -4.4627562855*10^-3, -2.7625748*10^-6, 
    7.30402014*10^-2, -2.2222592*10^-4, 9.42486*10^-8], {p, 34.376}] //
   Quiet // Timing

gives

(*
    {0.093601, {p -> 170.441}}
*)

So a speed-up of a factor of 2. Why not more? Because even though g2 compiles there is still a MainEvaluate (call-back to main Mathematica) to evaluate the Root object. Maybe things will get even faster if you use a C-compiler.

share|improve this answer
    
I think the performance improvement here is due to your choice to express the problem slightly differently rather than because of Compile. Calls out of the VM are very expensive and there are just too many of them here to make compilation of g2 worthwhile. (Although you forgot to localize i, it doesn't seem to make much difference overall.) Replacing Compile with Function actually improves the timings for me. –  Oleksandr R. Jan 28 '13 at 8:44
    
Thank's Oleksandr. Basically I just wanted to show how to get Compile to work at all. Did you test CompilationTarget -> "C" ? I have no time to investigate this further. Feel free to edit and extend my answer. –  Rolf Mertig Jan 28 '13 at 11:55
add comment

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.