How can I convert an implicit function to an explicit function?

I wish to convert an implicit function to an explicit one even if I sacrifice some accuracy.

What is the approximate expression for T(C)?

Best,

[Walter Roberson]

Is the left hand side a constant relative to T(c)? If so then the approximate solution is probably T(c) = 0 unless you are willing to narrow the ranges of consideration a fair bit.

[hns]

P(kPa) is proportional to T(c).

Will knowing that 50 <= T(c) <= 80 help?

[Walter Roberson]

Are some of the values constants? If so what are their values?

[hns]

P(kPa)= 70

y1= 0.6

y2= 0.4

A1=14.8950

A2=14.7513

B1=3413.10

B2=3331.7

C1=250.523

C2=227.6

T(c) ~= ? (Can you get it by explicit expression obtained from MATLAB?)

Answer 1

[Walter Roberson]

On my system at home, the following program takes less than 6 seconds. However, because I was iterating through versions during development, it might have "learned" solutions to some of the equations, so potentially it could take significantly longer the first time you run it. On the online facility, it takes more than 55 seconds, so I cannot directly show the result here.

The variable N controls the maximum taylor degree to use. If you change it to 5, then the code takes a fair number of minutes to run (about 20 minutes or so.) The best result for N = 5 is fairly close to the best result for N = 4, which executes much faster. It is possible that the best result for N = 5 is a a bit better than for N = 4, but the cost is so high that it does not seem to be worth doing.

For example it might be worth taking the best solution for N = 4 to use as the starting point for fsolve() to find exact values.

tic
Q = @(v) sym(v);

syms PkPa positive
y1 = Q(0.6);
y2 = Q(0.4);
A1 = Q(14.8950);
A2 = Q(14.7513);
B1 = Q(3413.10);
B2 = Q(3331.7);
C1 = Q(250.523);
C2 = Q(227.6);

pkpa = 70;

syms Tc

eqn = 1/PkPa == y1 / exp(A1 - B1/(Tc + C1)) + y2 / exp(A2 - B2/(Tc + C2))
E = simplify(eqn, 'steps', 50);
N = 4;
eqn_t = zeros(N, 1, 'sym');
sol_t = cell(N, 1);
sol_c = repmat({}, N, 1);
sol_rc = repmat({}, N, 1);
sol_selt = cell(N,1);

for K = 1 : N
     eqn_t(K) = (taylor(lhs(E), Tc, 65, 'order', K) == rhs(E));
     sol_t{K} = solve([eqn_t(K), imag(Tc)==0], Tc, 'returnconditions', true, 'maxdegree', 4);
     ncond = max(length(sol_t{K}.conditions),1);
     sol_c{K} = symfalse(ncond, 1);
     sol_rc{K} = symfalse(ncond, 1);
     sol_selt{K} = sym.empty;
     for J = 1 : length(sol_t{K}.conditions)
          cond = sol_t{K}.conditions(J);
          if isequal(cond, symtrue)
             sol_c{K}(J) = symtrue;
             sol_rc{K}(J) = symtrue;
             sol_selt{K}(end+1,1) = sol_t{K}.Tc(J);
          else
            temp = solve(sol_t{K}.conditions(J), 'returnconditions', true, 'maxdegree', 4);
            if isAlways(temp.conditions, 'unknown', 'false')
                sol_c{K}(J) = symtrue;
                sol_rc{J}(J) = symtrue;
                sol_selt{K}(end+1,1) = sol_t{K}.Tc(J);
            else
              sol_c{K}(J) = vpa(temp.conditions);
              sol_rc{K}(J) = simplify(subs(sol_c{K}(J), temp.parameters, pkpa));
              if sol_rc{K}(J)
                  sol_selt{K}(end+1,1) = sol_t{K}.Tc(J);
              end
            end
          end
     end
end

%%
PK = linspace(60,80,20);

orig = arrayfun(@(pk) vpasolve(subs(eqn, PkPa, pk), 65), PK);

plot(PK, orig, 'displayname', 'orig');
xlabel('P(kPa)')
ylabel('T(c)');

hold on

for K = 1 : N
    selt = sol_selt{K};
    for J = 1 : length(selt)
        this = double(subs(selt(J), PkPa, PK));
        plot(PK, this, '--', 'displayname', sprintf('O(%d) #%d', K, J));
    end
end
hold off
legend show
ylim auto

vpa(sol_selt{4}(1), 16)

toc

[hns]

Impressive!

Thanks,