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,
4 个评论
Walter Roberson
2021-5-15
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
2021-5-15
Walter Roberson
2021-5-15
Are some of the values constants? If so what are their values?
hns
2021-5-15
采纳的回答
Walter Roberson
2021-5-16
编辑：Walter Roberson
2021-5-16
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
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Axis Labels 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
