0 个投票

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 ?

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

 采纳的回答

Walter Roberson
Walter Roberson 2021-10-19
Ran in:

1 个投票

Ran in:
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
Alessio Falcone 2021-10-19
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 2021-10-19
Thank you very much Mr. Roberson

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Calculus 的更多信息

产品

版本

R2021b

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by