3
$\begingroup$

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

(*Additional functions*)
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 interactionTable[temp] for temp=0.1

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

enter image description here

Is there any way to get rid of those peaks?

Many thanks! I appreciate any suggestions.

$\endgroup$
1
  • $\begingroup$ You may need higher WorkingPrecision in your numerical calculations. $\endgroup$
    – MarcoB
    Feb 20 at 2:23

1 Answer 1

4
$\begingroup$

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

enter image description here

ListLinePlot[gammaInteractionTable[1/10]]

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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