How to optimize an reduce the computation time

while Fg >= 0.001
    if Fg > 0.15
        Psat = Psat + 1;
    elseif Fg > 0.02
        Psat = Psat + 0.1;
    elseif Fg > 0.009
        Psat = Psat + 0.01;
    else
        Psat = Psat + 0.001;
    end
  Pri = Psat./PcMPa;
  k = exp(5.37.*(1+w).*(1-Tri.^(-1)))./Pri;
  syms fg
  fun = sum(x.*(k-1)./(1+fg*(k-1)));
  Fgmass = double(solve(fun==0,fg));
  Fg = Fgmass(Fgmass<0.99 & Fgmass>0);
end

Hello. I am a complete beginner. How can i reduce computation time for this process? Now it takes around a minute. I need to do 2-20 thousand of such computations. Thank you in advance.

[Yash]

Hi @Ruslan

Please share the values of each of these variables to reproduce the issue: Fg, Psat, PcMPa, Tri.

[Ruslan]

x = [0.0035, 0.0314, 0.5426, 0.0857, 0.0572, 0.0076, 0.0245, 0.0075, 0.012, 0.0153, 0.026, 0.0302, 0.021, 0.0174, 0.0136, 0.011, 0.0111, 0.0095, 0.0086, 0.0068, 0.006, 0.0056, 0.0051, 0.0408];
PcMPa = [3.399, 7.382, 4.604, 4.880, 4.249, 3.648, 3.797, 3.381, 3.369, 3.012, 2.736, 2.486, 2.289, 2.096, 1.979, 1.824, 1.7, 1.524, 1.486, 1.4, 1.325, 1.294, 1.173, 1.08];
TcK= [126.3, 304.2, 190.6, 305.4, 369.8, 408.2, 425.2, 460.4, 469.7, 507.4, 540.3, 568.8, 594.6, 617.7,639, 658, 675, 693, 708, 723, 734, 748, 755, 768];
w =[0.045, 0.231, 0.0115, 0.0908, 0.1454, 0.1756, 0.1928, 0.2273, 0.251, 0.2957, 0.3506, 0.3978, 0.4437, 0.4902, 0.535, 0.575, 0.619, 0.659, 0.706, 0.742, 0.77, 0.79, 0.827, 0.907];
TK = 200;
R = 8.314;
Psat = 0.01;
Tri = TK./TcK;
Fg = 100;
while Fg >= 0.001
    if Fg > 0.15
        Psat = Psat + 1;
    elseif Fg > 0.02
        Psat = Psat + 0.1;
    elseif Fg > 0.009
        Psat = Psat + 0.01;
    else
        Psat = Psat + 0.001;
    end
  Pri = Psat./PcMPa;
  k = exp(5.37.*(1+w).*(1-Tri.^(-1)))./Pri;
  syms fg
  fun = sum(x.*(k-1)./(1+fg*(k-1)));
  Fgmass = double(solve(fun==0,fg));
  Fg = Fgmass(Fgmass<0.99 & Fgmass>0);
end

[Walter Roberson]

  k = exp(5.37.*(1+w).*(1-Tri.^(-1)))./Pri;

In that calculation, w and Tri are calculated before the while loop. The entire expression exp(5.37.*(1+w).*(1-Tri.^(-1))) is effectively constant . So pre-calculate it,

  wexp = exp(5.37.*(1+w).*(1-Tri.^(-1)));
  

before the loop, and then inside the loop

  k = wexp ./ Pri;

[Walter Roberson]

  syms fg
  fun = sum(x.*(k-1)./(1+fg*(k-1)));
  Fgmass = double(solve(fun==0,fg));
  Fg = Fgmass(Fgmass<0.99 & Fgmass>0);

Are you expecting multiple solutions for the solve() ? If not then it would be better to use vpasolve() instead of solve()

Answer 1

[Torsten]

I don't understand the sense of your while-loop.

If you want to define Psat as a function of fg, the problem becomes more complicated.

Try to use "fzero" or "fsolve" instead of symbolic computations like done below if you want to make the code perform faster.

x = [0.0035, 0.0314, 0.5426, 0.0857, 0.0572, 0.0076, 0.0245, 0.0075, 0.012, 0.0153, 0.026, 0.0302, 0.021, 0.0174, 0.0136, 0.011, 0.0111, 0.0095, 0.0086, 0.0068, 0.006, 0.0056, 0.0051, 0.0408];
PcMPa = [3.399, 7.382, 4.604, 4.880, 4.249, 3.648, 3.797, 3.381, 3.369, 3.012, 2.736, 2.486, 2.289, 2.096, 1.979, 1.824, 1.7, 1.524, 1.486, 1.4, 1.325, 1.294, 1.173, 1.08];
TcK= [126.3, 304.2, 190.6, 305.4, 369.8, 408.2, 425.2, 460.4, 469.7, 507.4, 540.3, 568.8, 594.6, 617.7,639, 658, 675, 693, 708, 723, 734, 748, 755, 768];
w =[0.045, 0.231, 0.0115, 0.0908, 0.1454, 0.1756, 0.1928, 0.2273, 0.251, 0.2957, 0.3506, 0.3978, 0.4437, 0.4902, 0.535, 0.575, 0.619, 0.659, 0.706, 0.742, 0.77, 0.79, 0.827, 0.907];
TK = 200;
R = 8.314;
PPsat = 0.01:0.01:3.38;
Tri = TK./TcK;
fg0 = 0.7;
for i = 1:numel(PPsat)
  Psat = PPsat(i);
  Pri = Psat./PcMPa;
  k = exp(5.37.*(1+w).*(1-Tri.^(-1)))./Pri;
  %syms fg
  %sol = vpa(solve(fun==0,fg));
  %sol = sol(abs(imag(sol))<1e-6);
  %sol = sol(sol>0 & sol<1);
  %Fg(i) = double(sol);
  fun = @(fg)sum(x.*(k-1)./(1+fg*(k-1)));
  Fg(i) = fsolve(fun,fg0,optimset('Display','none'));
  fg0 = Fg(i);
end
plot(PPsat,Fg)
xlabel('Psat')
ylabel('fg')

[Ruslan]

My while loop is trying to find value of Psat at which fg as close to zero as possible. I don't know the end value of Psat even remotely. It can be anything from 2 to 100.

[Stephen23]

Do not reinvent the wheel.

As Torsten wrote, just use FZERO or FSOLVE or similar.

Avoid symbolic mathematics (if computation time is important).

[Torsten]

My while loop is trying to find value of Psat at which fg as close to zero as possible

Then insert fg = 0 (or something close to 0) and solve for Psat:

x = [0.0035, 0.0314, 0.5426, 0.0857, 0.0572, 0.0076, 0.0245, 0.0075, 0.012, 0.0153, 0.026, 0.0302, 0.021, 0.0174, 0.0136, 0.011, 0.0111, 0.0095, 0.0086, 0.0068, 0.006, 0.0056, 0.0051, 0.0408];
PcMPa = [3.399, 7.382, 4.604, 4.880, 4.249, 3.648, 3.797, 3.381, 3.369, 3.012, 2.736, 2.486, 2.289, 2.096, 1.979, 1.824, 1.7, 1.524, 1.486, 1.4, 1.325, 1.294, 1.173, 1.08];
TcK= [126.3, 304.2, 190.6, 305.4, 369.8, 408.2, 425.2, 460.4, 469.7, 507.4, 540.3, 568.8, 594.6, 617.7,639, 658, 675, 693, 708, 723, 734, 748, 755, 768];
w =[0.045, 0.231, 0.0115, 0.0908, 0.1454, 0.1756, 0.1928, 0.2273, 0.251, 0.2957, 0.3506, 0.3978, 0.4437, 0.4902, 0.535, 0.575, 0.619, 0.659, 0.706, 0.742, 0.77, 0.79, 0.827, 0.907];
TK = 200;
R = 8.314;
Tri = TK./TcK;
Pri = @(Psat)Psat./PcMPa;
k = @(Psat)exp(5.37.*(1+w).*(1-Tri.^(-1)))./Pri(Psat);
fg = 0;
fun = @(Psat)sum(x.*(k(Psat)-1)./(1+fg*(k(Psat)-1)));
Psat = fsolve(fun,5)

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

Psat = 3.3467