# Removing fluctuations in plot, possibly due to numerical precision

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][];(*p1=
p2*)
pressp4[β] =
P /. Solve[(tempPV - T4 == 0) /. V -> V4 /. parap4, P][];(*p3=
p4*)
Tp1[β] =
T1 /. Solve[(tempPV - T1 == 0) /. V -> V4 /. parap1 /.
P -> pressp2[β], T1][];
, {β, β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: 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?

• 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. – MarcoB May 15 at 18:22
• @MarcoB Yes you are right. This saves a lot of time. – Soodeh Z. May 15 at 20:38

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
]
] 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
] 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];

• @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. – Soodeh Z. May 15 at 20:42