# Peak Appearances in Plot Outputs from Optimization Routines

Dear fellow researchers and Mathematica enthusiasts,

I've been working on a computational physics problem and encountered an interesting issue regarding the appearance of peaks in my plots, which I generated using NMinimize. The core of my work involves evaluating certain integrals and their contribution to a potential function, which is then minimized to find equilibrium states under specific parameters.

For clarity, I've abstracted my code as follows:

(*Define the energy density function*)

Clear[EnergyDensity];
EnergyDensity[eta_, mu_, T_] := (eta^3)/(6  Pi) + (eta/(2  Pi))  T^2  (PolyLog[2, -Exp[-(eta - mu)/T]] + PolyLog[2, -Exp[-(eta + mu)/T]]) + (T^3/(2  Pi))  (PolyLog[ 3, -Exp[-(eta - mu)/T]] + PolyLog[3, -Exp[-(eta + mu)/T]]);

(*Optimize the second integral*)

Clear[SecondIntegral];

SecondIntegral[eta_, mu_, T_] := -(T/(4  Pi))  ((eta/T) + Log[1 + Exp[-(eta - mu)/T]] + Log[1 + Exp[-(eta + mu)/T]]);

(*Simplify the third integral using the properties of PolyLog and Log \
functions*)
Clear[ThirdIntegral];

ThirdIntegral[eta_, mu_, T_] := (T^2/(4  Pi))  ((eta/T)  Log[(1 + Exp[(eta + mu)/T])/(1 + Exp[(eta - mu)/T])] + PolyLog[2, -Exp[(eta + mu)/T]] - PolyLog[2, -Exp[(eta - mu)/T]]);

(*effective potential function*)

Clear[EffectivePotential];

EffectivePotential[mu_, v0_, v_, eta_, delta_, lambda_, n_, T_] := -((delta  v0^2)/(2  lambda)) + 2  EnergyDensity[eta, mu, T] + 4  delta  eta^2  SecondIntegral[eta, mu, T] + 4  v ThirdIntegral[eta, mu, T] + (8  (delta  lambda/n)  eta^2  (SecondIntegral[eta, mu, T])^2 - 4  (delta  lambda/n)  (ThirdIntegral[eta, mu, T])^2);

Clear[PotentialBar, InteractionEnergy, GammaInteractionEnergy, TotalPotential, Solution, EtaBar, K0Bar];

PotentialBar[mu_, k0_, eta_, lambda_, T_] :=  4  lambda  ThirdIntegral[eta, mu + k0, T];

InteractionEnergy[mu_, k0_, eta_, lambda_, T_] := -4  eta  SecondIntegral[eta, mu + k0, T];

GammaInteractionEnergy[mu_, k0_, eta_, lambda_, T_] := -4  ThirdIntegral[eta, mu + k0, T];

TotalPotential[mu_, k0_, v0_, eta_, delta_, lambda_, n_, T_] := EffectivePotential[Abs[mu + k0], v0, v0 - k0, Abs[eta], delta, lambda, n, T];

(*Solution using the NMinimize*)

Solution[mu_, lambda_, n_, T_] := NMinimize[TotalPotential[mu, k0, PotentialBar[mu, k0,eta, lambda, T], eta, 1, lambda, n, T], {eta, k0}][[2]];

(*Values for eta and k0*)
EtaBar[mu_, lambda_, n_, T_] := eta /. Solution[mu, lambda, n, T];
K0Bar[mu_, lambda_, n_, T_] := k0 /. Solution[mu, lambda, n, T];

(*Generate tables for interaction energies*)
interactionTable[temp_] := Table[{x, InteractionEnergy[x, K0Bar[x, Pi, 3, temp], Abs[EtaBar[x, Pi, 3, temp]], Pi, temp]}, {x, 0, 9, 0.05}];

gammaInteractionTable[temp_] := Table[{x, GammaInteractionEnergy[x, K0Bar[x, Pi, 3, temp], Abs[EtaBar[x, Pi, 3, temp]], Pi, temp]}, {x, 0, 9, 0.05}];

ListLinePlot[interactionTable[0.1]]
ListLinePlot[gammaInteractionTable[0.1]]



The plots I got are in the images below:

-Plot of interactionTable[temp] for temp=0.1:

-Plot of gammaInteractionTable[temp] for temp=0.1:

Is there any way to get rid of those peaks?

Many thanks! I appreciate any suggestions.

• You may need higher WorkingPrecision in your numerical calculations. Feb 20 at 2:23

Specify a WorkingPrecision for NMinimize to avoid using machine precision.

(*Solution using the NMinimize*)

Solution[mu_, lambda_, n_, T_] :=
NMinimize[
TotalPotential[mu, k0, PotentialBar[mu, k0, eta, lambda, T], eta, 1,
lambda, n, T], {eta, k0},
WorkingPrecision -> 10][[2]];

(*Values for eta and k0*)
EtaBar[mu_, lambda_, n_, T_] := eta /. Solution[mu, lambda, n, T];
K0Bar[mu_, lambda_, n_, T_] := k0 /. Solution[mu, lambda, n, T];


Use exact values for constant values

(*Generate tables for interaction energies*)
interactionTable[temp_] :=
Table[{x,
InteractionEnergy[x, K0Bar[x, Pi, 3, temp], Abs[EtaBar[x, Pi, 3, temp]],
Pi, temp]}, {x, 0, 9, 1/20}];

gammaInteractionTable[temp_] :=
Table[{x,
GammaInteractionEnergy[x, K0Bar[x, Pi, 3, temp],
Abs[EtaBar[x, Pi, 3, temp]], Pi, temp]}, {x, 0, 9, 1/20}];

ListLinePlot[interactionTable[1/10]]


ListLinePlot[gammaInteractionTable[1/10]]