Empty sym: 0-by-1

Question

Alessio Falcone 2021-10-19

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1567338-empty-sym-0-by-1

评论： Alessio Falcone 2021-10-19

采纳的回答： Walter Roberson

在 MATLAB Online 中打开

Hi everyone !

I have some problems with a system of equations that I can't manage to solve. Could you please help me ?

Tc=400;
Pc=35;
R=0.0821;
Z=0.34;
syms Vc a b c;
X=R*Tc*(1-c)/(Vc)-(a*b*(c^2)/(Vc^2))-(R*Tc/b)*log(1-((c*b)/Vc))-Pc;
eqn1=R*Tc*(1-c)/(Vc)-(a*b*(c^2)/(Vc^2))-(R*Tc/b)*log(1-((c*b)/Vc))-Pc==0
eqn2=diff(X,Vc)==0
eqn3=diff(X,Vc,2)==0
eqn4=(Pc*Vc)/(R*Tc)-Z==0
sol=vpasolve([eqn1,eqn2,eqn3,eqn4],[a,b,c,Vc]);
Vc=sol.Vc
a=sol.a
b=sol.b
c=sol.c

Could you please help me ?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Walter Roberson 2021-10-19

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1567338-empty-sym-0-by-1#answer_812243

在 MATLAB Online 中打开

Q = @(v) sym(v);

Tc = Q(400);

Pc = Q(35);

R = Q(0.0821);

Z = Q(0.34);

syms Vc a b c;

%assume([b, c], 'real')

assumeAlso(b ~= 0 & c ~= 0)

X = R*Tc*(1-c)/(Vc)-(a*b*(c^2)/(Vc^2))-(R*Tc/b)*log(1-((c*b)/Vc))-Pc;

eqn1 = R*Tc*(1-c)/(Vc)-(a*b*(c^2)/(Vc^2))-(R*Tc/b)*log(1-((c*b)/Vc))-Pc == 0

eqn1 =

eqn2 = diff(X,Vc) == 0

eqn2 =

eqn3 = diff(X,Vc,2) == 0

eqn3 =

eqn4 = (Pc*Vc)/(R*Tc)-Z == 0

eqn4 =

eqns = simplify([eqn1, eqn2, eqn3, eqn4])

eqns =

Vc_sol = solve(eqns(end), Vc)

Vc_sol =

eqns2 = simplify(subs(eqns(1:end-1), Vc, Vc_sol))

eqns2 =

partial_a = solve(eqns2(2), a)

partial_a =

eqns3 = simplify(subs(eqns2([1 3:end]), a, partial_a))

eqns3 =

partial_b = solve(eqns3(2), b, 'returnconditions', true)

partial_b = struct with fields:

b: [2×1 sym] parameters: [1×0 sym] conditions: [2×1 sym]

partial_b.b

ans =

partial_b.conditions

ans =

eqns4_1 = (subs(eqns3([1 3:end]), b, partial_b.b(1)))

eqns4_1 =

eqns4_2 = (subs(eqns3([1 3:end]), b, partial_b.b(2)))

eqns4_2 =

sol_c_1 = vpasolve(eqns4_1, c)

sol_c_1 = Empty sym: 0-by-1

sol_c_2 = vpasolve(eqns4_2, c)

sol_c_2 =

141.68652558079973320034185988983

full_c = sol_c_2

full_c =

141.68652558079973320034185988983

full_b = subs(partial_b.b(2), c, full_c)

full_b =

0.00017449682895104484527054335897185

full_a = subs(subs(partial_a, b, partial_b.b(2)), c, full_c)

full_a =

19.294808147837450107863804722783

full_Vc = Vc_sol

full_Vc =

sol = [full_a, full_b, full_c, full_Vc]

sol =

subs([eqn1, eqn2, eqn3, eqn4], [a, b, c, Vc], sol)

ans =

vpa(ans)

ans =

So the solution works to within round-off error.

3 个评论
显示 1更早的评论隐藏 1更早的评论

Alessio Falcone 2021-10-19

在 MATLAB Online 中打开

Q = @(v) sym(v);
Tc = Q(400);
Pc = Q(35);
R = Q(0.0821);
Z = Q(0.34);
syms Vc a b c;
%assume([b, c], 'real')
assumeAlso(b ~= 0 & c ~= 0)
X = R*Tc*(1-c)/(Vc)-(a*b*(c^2)/(Vc^2))-(R*Tc/b)*log(1-((c*b)/Vc))-Pc;
eqn1 = R*Tc*(1-c)/(Vc)-(a*b*(c^2)/(Vc^2))-(R*Tc/b)*log(1-((c*b)/Vc))-Pc == 0
eqn2 = diff(X,Vc) == 0
eqn3 = diff(X,Vc,2) == 0
eqn4 = (Pc*Vc)/(R*Tc)-Z == 0
eqns = simplify([eqn1, eqn2, eqn3, eqn4])
Vc_sol = solve(eqns(end), Vc)
eqns2 = simplify(subs(eqns(1:end-1), Vc, Vc_sol))
partial_a = solve(eqns2(2), a)
eqns3 = simplify(subs(eqns2([1 3:end]), a, partial_a))
partial_b = solve(eqns3(2), b, 'returnconditions', true)
partial_b.b
partial_b.conditions
eqns4_1 = (subs(eqns3([1 3:end]), b, partial_b.b(1)))
eqns4_2 = (subs(eqns3([1 3:end]), b, partial_b.b(2)))
sol_c_1 = vpasolve(eqns4_1, c)
sol_c_2 = vpasolve(eqns4_2, c)
full_c = sol_c_2
full_b = subs(partial_b.b(2), c, full_c)
full_a = subs(subs(partial_a, b, partial_b.b(2)), c, full_c)
full_Vc = Vc_sol
sol = [full_a, full_b, full_c, full_Vc]
subs([eqn1, eqn2, eqn3, eqn4], [a, b, c, Vc], sol)
vpa(ans)
syms V 
k=1;
for T=300:1:400;
P=R*T*(1-full_c)/(V)-(full_a*full_b*(full_c^2)/(V^2))-(R*T/full_b)*log(1-((full_c*full_b)/V));
derivataprima= diff(P,V);
soluzione_tmp =vpasolve(derivataprima,V);
soluzione=soluzione_tmp(soluzione_tmp>0);
Vmin(k)=min(soluzione);
Pmin(k)=R*T*(1-full_c)/(Vmin(k))-(full_a*full_b*(full_c^2)/(Vmin(k)^2))-(R*T/full_b)*log(1-((full_c*full_b)/Vmin(k)));
Vmax(k)=max(soluzione);
Pmax(k)= R*T*(1-full_c)/(Vmax(k))-(full_a*full_b*(full_c^2)/(Vmax(k)^2))-(R*T/full_b)*log(1-((full_c*full_b)/Vmax(k)));
inver=fliplr(Pmax);
k=k+1;
unione_V=union(Vmin,Vmax);
unione_P= union(Pmin,inver,'stable');
plot(unione_V,unione_P,'k .')
axis ([0.01 1 -100 100])
fplot(P)
hold on
end