1
$\begingroup$

Through the following code, we generate Tp1:

tempPV = (
  3 π^(2/3) + 6 6^(2/3) P π V^(2/3) - 
   6^(2/3) V^(
    2/3) (-3 + Sqrt[
      9 + (4 6^(2/3) π^(4/3) q^2)/(
       V^(4/3) β^2)]) β^2 + 
   3 6^(2/3) V^(
    2/3) β^2 Log[
     1/6 (3 + Sqrt[
        9 + (4 6^(2/3) π^(4/3) q^2)/(V^(4/3) β^2)])])/(
  6 6^(1/3) π^(4/3) V^(1/3));

βin = 0.01;
βfi = 100;
βst = 0.005;
Table[
  xlog[β] = Log[10, β]
  , {β, βin, βfi, βst}];

parap1 = {q -> 0.1, V4 -> 5000};
parap2 = {T2 -> 15, q -> 0.1, V2 -> 10000};
parap4 = {T4 -> 5, q -> 0.1, V4 -> 5000};

Table[
  pressp2[β] = 
   P /. Solve[(tempPV - T2 == 0) /. V -> V2 /. parap2, P][[1]];(*p1=
  p2*)
  pressp4[β] = 
   P /. Solve[(tempPV - T4 == 0) /. V -> V4 /. parap4, P][[1]];(*p3=
  p4*)
  Tp1[β] = 
   T1 /. Solve[(tempPV - T1 == 0) /. V -> V4 /. parap1 /. 
       P -> pressp2[β], T1][[1]];
  , {β, βin, βfi, βst}];

mi = Min[Table[Tp1[β], {β, βin, βfi, βst}]]
ma = Max[Table[Tp1[β], {β, βin, βfi, βst}]]

ListPlot[
 Table[{xlog[β], Tp1[β]}, {β, βin, βfi, βst}], 
 ScalingFunctions -> {Rescale[#, {mi, ma}, {0.`, 1.`}] &, 
   Rescale[#, {0.`, 1.`}, {mi, ma}] &}, Joined -> True, Frame -> True,
  FrameStyle -> Black, 
 BaseStyle -> {FontSize -> 14, PrintPrecision -> 11}, 
 FrameLabel -> {"\!\(\*SubscriptBox[\(log\), \(10\)]\) (β)", 
   "\!\(\*SubscriptBox[\(T\), \(1\)]\)"}, RotateLabel -> False, 
 PlotStyle -> {Blue, Thickness[0.006]}, 
 PlotRange -> {{-2, 2}, {mi, ma}}, Axes -> None, AspectRatio -> 0.8, 
 ImageSize -> 400, FrameTicks -> {{ticks, None}, {Automatic, None}}] 

The result is the following plot:

plot with fuzzy tail

There is a strange fluctuation for $log_{10}^{\beta}=1-2$. As it should be a smoothly decreasing function, what is the origin of these fluctuations? How to fix this possibly numerical error?

$\endgroup$
  • 3
    $\begingroup$ You should really consider radically refactoring your code to take advantage of vectorized operations. First of all, remove all the [beta] indices and do not iterate over the values of beta. For instance, Table[xlog[β] = Log[10, β], {β, βin, βfi, βst}]; should really just be xlog = Log10@Range[βin, βfi, βst];. Second, pre-calculate the symbolic solution of the Solve equations in your table ONLY ONCE and then plug in values of beta to get your vectors of results. This will save A LOT of time. $\endgroup$ – MarcoB May 15 at 18:22
  • $\begingroup$ @MarcoB Yes you are right. This saves a lot of time. $\endgroup$ – Soodeh Z. May 15 at 20:38
3
$\begingroup$

Here is a more idiomatic approach:

para = {q -> 1/10, V4 -> 5000, T2 -> 15, V2 -> 10000, T4 -> 5};
Psol = First@Solve[tempPV - T2 == 0 /. V -> V2 /. para, P];
T1sol = T1 /. First@Solve[tempPV - T1 == 0 /. V -> V4 /. para, T1] /. Psol;

Block[
  {tmin = T1sol /. β -> βfi, tmax = T1sol /. β -> βin},
  LogLinearPlot[
    (T1sol - tmin)/(tmax - tmin), {β, βin, βfi},
    PlotRange -> All, WorkingPrecision -> 20,
    Axes -> False, Frame -> True
  ]
]

plot from LogLinear plot, appropriately scaled


Here is a completely refactored version of your code that runs quite a bit faster even though the calculations are carried out at much higher arbitrary precision:

(* Re-defined these as arbitrary-precision numbers rather than machine-precision reals *)
βin = 1/100; βfi = 100; βst = 5/1000;

(* Vectorize and avoid Table *)
xlog = Log10@Range[βin, βfi, βst];

(* grouped all parameters together so you only have one location to change *)
(* also re-defined the values as arbitrary-precision numbers               *)
para = {q -> 1/10, V4 -> 5000, T2 -> 15, V2 -> 10000, T4 -> 5};

(* Solve all equations ONLY ONCE, then plug in values when needed          *)
pressp2 = P /. First@Solve[(tempPV - T2 == 0) /. V -> V2 /. para, P];
Tp1 = T1 /. First@Solve[(tempPV - T1 == 0) /. V -> V4 /. para, T1] /. P -> pressp2;

(* Plug in the values of beta *)
(* Using 20-digit arbitrary-precision numerical values; fewer digits still gives artifacts *)
Tp1values = Tp1 /. β -> N[Range[βin, βfi, βst], 20];

(* Rescale values of Tp1 to run between 0 and 1                       *)
(* Pair the rescaled values with the corresponding abscissa from xlog *)
rescaledPaired = Transpose@{xlog, Rescale[Tp1values]};

(* Generate plot *)
ListLinePlot[
  rescaledPaired,
  PlotRange -> All, PlotRangePadding -> Scaled[0.05],
  Axes -> False, Frame -> True
]

expected plot of a monotonically decreasing function, without the fuzzy tail


Note that pressp4 was never needed to generate your plot, so I removed it from above. Nevertheless it would be obtained as follows:

pressp4 = P /. First@Solve[(tempPV - T4 == 0) /. V -> V4 /. para, P];
| improve this answer | |
$\endgroup$
  • $\begingroup$ @MacroB Wow! How a nice code. I learnt a lot from it. Thanks a lot for your time and consideration. You helped me a lot. $\endgroup$ – Soodeh Z. May 15 at 20:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.