Empty sym: 0-by-1

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 ?

Answer 1

[Walter Roberson]

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.

[Alessio Falcone]

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

[Alessio Falcone]

[Alessio Falcone]

Thank you very much Mr. Roberson