あなたは私の最後の希望です。私の大学には、私の質問に答えられる人はいません。私は6つのパラメータa0、a1、a2、b0、b1、b2に依存する非常に複雑な関数を持っています。これは、かなり新しい状態方程式によって計算された圧力、体積液体、体積蒸気のデルタを最小化します。NMinimize は非常に遅く、タイミングが非常に高いため、この式については考慮できませんでした。コードには、私のコードに関するいくつかの説明といくつかの問題があります。ひざまずいて、あなたが私を助けてくれることを祈っています。申し訳ありませんが、これらの方程式の構築に 4 か月を費やした後、テストできませんでした。そしてフラストレーションは日に日に増していきます！
Clear["Global`*"];
data = {{100., 34.376, 0.036554, 23.782}, {105., 56.377, 0.037143,
15.116}, {110., 88.13, 0.037768, 10.038}, {115., 132.21, 0.038431,
6.9171}, {120., 191.43, 0.039138, 4.9183}, {125., 268.76,
0.039896, 3.5915}, {130., 367.32, 0.040714, 2.6825}, {135.,
490.35, 0.0416, 2.0424}, {140., 641.18, 0.042569, 1.5803}, {145.,
823.22, 0.043636, 1.2393}, {150., 1040., 0.044825,
0.98256}, {155., 1295., 0.046165, 0.78568}, {160., 1592.1,
0.047702, 0.63206}, {165., 1935.1, 0.0495, 0.51014}, {170.,
2328.3, 0.051667, 0.41163}, {175., 2776.5, 0.054394,
0.33038}, {180., 3285.2, 0.058078, 0.26139}, {185., 3861.7,
0.063825, 0.19945}, {190., 4518.6, 0.079902, 0.12816}};
tvector = data[[All, 1]];(*K*)
pvector =
data[[All, 2]];(*KPa*)
vlvector = data[[All, 3]];(*L/mol*)
vvvector =
data[[All, 4]];
(*L/mol.*)
r = 8.314472;
tc = 190.56;
avvicinamento = Length[tvector] - 3;
trexp = Take[tvector, avvicinamento]/tc;
vlexp = Take[vlvector, avvicinamento];
vvexp = Take[vvvector, avvicinamento];
zeri = Table[i*0., {i, avvicinamento}];
pexp = Take[pvector, avvicinamento];
(*Function for calculation of Fugacity of CSD Equation*)
(*Function for calculation of Fugacity of CSD Equation*)
fug[v_, p_, t_, a_, b_] :=
Module[{y, z, vbv, vb, f1, f2, f3, f4, f}, y = b/(4 v);
z = (p v)/(r t);
vbv = Log[(v + b)/v];
vb = v + b;
f1 = (4*y - 3*y^2)/(1 - y)^2;
f2 = (4*y - 2*y^2)/(1 - y)^3;
f3 = (2*vbv)/(r t*b)*a;
f4 = (vbv/b - 1/vb)/(r t)*a;
f = f1 + f2 - f3 + f4 - Log[z];
Exp[f]]
(*g Minimize the equality of fugacity*)
g[p_?NumericQ, t_?NumericQ, a0_?NumericQ, a1_?NumericQ, a2_?NumericQ,
b0_?NumericQ, b1_?NumericQ, b2_?NumericQ] := Module[{},
a = a0*Exp[a1*t + a2*t^2];
b = b0 + b1*t + b2*t^2;
csd = a/(r*t*(b + v)) - (-(b^3/(64*v^3)) + b^2/(16*v^2) +
b/(4*v) + 1)/(1 - b/(4*v))^3 + (p*v)/(r*t);
vol = NSolve[csd == 0 && v > 0, v, Reals];
sol = v /. vol;
(*If[Length[sol]==1,Interrupt[];Print["Sol==1"]];*)
vliquid = Min[sol];
vvapor = Max[sol];
fl = fug[vliquid, p, t, a, b];
fv = fug[vvapor, p, t, a, b];
(*Print[{t,p,vol,Abs[fl-fv]}];*)
Abs[fl - fv]];
(*This function minimize the pcalc-pexp and vcalc-vexp *)
hope[a0_?NumericQ, a1_?NumericQ, a2_?NumericQ, b0_?NumericQ,
b1_?NumericQ, b2_?NumericQ] :=
Module[{},
pp[a0, a1, a2, b0, b1, b2] :=
Table[FindRoot[{g[p, tvector[[i]], a0, a1, a2, b0, b1, b2]},
{p,pvector[[i]]}],{i,avvicinamento}];
pressioni1 = pp[a0, a1, a2, b0, b1, b2];
pcalc = p /. pressioni1;
differenza = ((pcalc - pexp)/pexp)^2;
If[MemberQ[differenza, 0.],
differenza = zeri + RandomReal[{100000, 500000}];(*
First problem:
As I've FindRoot that finds the solutions equal to the starting \
point, I don't want these kind of solutions and with this method - \
+RandomReal[{100000,500000}] -
a keep away this solutions.Is it right? *)
deltap = Total[differenza],
differenzanonzero = Select[differenza, # > 0 &];
csd1[a_, b_, p_, t_] :=
a/(r*t*(b + v)) - (-(b^3/(64*v^3)) + b^2/(16*v^2) + b/(4*v) +
1)/(1 - b/(4*v))^3 + (p*v)/(r*t);(*Funzione CSD*)
volumi =
Table[NSolve[csd1[a, b, pcalc[[i]], tvector[[i]]], v, Reals], {i,
avvicinamento}];
soluzioni = v /. volumi;
vvcalc = Table[Max[soluzioni[[i]]], {i, avvicinamento}];
vlcalc = Table[Min[soluzioni[[i]]], {i, avvicinamento}];
deltavl = Total[((vlexp - vlcalc)/vlcalc)^2];
deltavv = Total[((vvexp - vvcalc)/vvcalc)^2];
deltap = Total[differenza];
Print[a0, " ", b0, " ", delta];
delta = 0.1*deltavl + 0.1*deltavv + deltap]];
NMinimize[{hope[a0, a1, a2, b0, b1, b2],
500 < a0 < 700 && -0.01 < a1 < -1.0*10^-5 && -10^-5 < a2 < -10^-7 &&
0.0010 < b0 < 0.1 && -0.0010 < b1 < -1.0*10^-5 &&
10^-9 < b2 < 10^-7}, {a0, a1, a2, b0, b1, b2}]
前もって感謝します！Mariano Pierantozzi 化学工学の博士課程の学生