I'd like to fit a trap potential
V4Sym[x_, a_, L_, c1_]
to a combination of two (or possibly more, but for now two) Gaussians:
doubleGauss[x_, in1_, in2_, x1_, x2_, σ1_, σ2_] :=
in1 * Exp[-(x - x1)^2/σ1^2] + in2*Exp[-(x - x2)^2/σ2^2]
After a quick search, I came acrosss this solution, which looked promising. However, if I try this, I get the following:
L := 10^-5
Lε := (L - 10^-8*L)/2
a := 10^8
c1 := 1
doubleGauss[x_, in1_, in2_, x1_, x2_, σ1_, σ2_] :=
in1*Exp[-(x - x1)^2/σ1^2] + in2*Exp[-(x - x2)^2/σ2^2]
V4Sym[x_, a_, L_, c1_]:=c1*(Coth[a*(L+x)/2]*Coth[a*(L - x)/2])^2 - 2*a^2*(Sech[a*(L+x)/2]^2 + Sech[a*(L - x)/2]^2)
distance := Integrate[(V4Sym[x, a, L, c1] - doubleGauss[x, in1, in2, x1, x2, σ1, σ2])^2,
{x, -Lε, Lε}, Assumptions -> (in1 | in2 | x1 | x2) ∈ Reals]
{min, sol} = NMinimize[distance, {in1, in2, x1, x2, σ1, σ2}]
The output I then get is
PolynomialGCD::lrgexp: Exponent is out of bounds for function PolynomialGCD.
a few quick remarks: L, a and c1 are a set of parameters for which I've chosen values in the range of what I'll need in the end. Lε is for the integration boundaries since the potential diverges at |x|=L/2.
Does anyone have some hints on what I'm doing wrong?
EDIT: forgot the definition of the potential.
EDIT2: sorry, made a mistake in defining the potential here (I had defined it implicitly in my .nb and did not want to include that step here since it did not seem relevant) the actual potential is
V4Sym[x_, a_, L_, c1_] :=
c1*(Coth[a*(L/2 + x)]*Coth[a*(L/2 - x)])^2 -
2*a^2*(Sech[a*(L/2 + x)]^2 + Sech[a*(L/2 - x)]^2)
So, as a whole, to plot it, we have
L := 10^-5
L\[CurlyEpsilon] := (L - 10^-8*L)/2
a := 10^8
c1 := 1
V4Sym[x_, a_, L_, c1_] :=
c1*(Coth[a*(L/2 + x)]*Coth[a*(L/2 - x)])^2 -
2*a^2*(Sech[a*(L/2 + x)]^2 + Sech[a*(L/2 - x)]^2)
Plot[V4Sym[x, a, L, c1], {x, -L\[CurlyEpsilon], L\[CurlyEpsilon]}]
I apologize, my potential seems to behave weirdly for this small parameter L (although I can't seem to make out why exactly). Let's for now stick with L=c1=1. This is what the shape of the potential looks like in this case:
V4Sym[x_, a_, L_, c1_] :=
c1*(Coth[a*(L/2 + x)]*Coth[a*(L/2 - x)])^2 -
2*a^2*(Sech[a*(L/2 + x)]^2 + Sech[a*(L/2 - x)]^2)
Plot[{Evaluate[V4Sym[x, 1, 1, 1]], V4Sym[x, 5, 1, 1],
V4Sym[x, 10, 1, 1], V4Sym[x, 30, 1, 1]}, {x, -0.5, 0.5},
PlotLegends -> {"a=1", "a=5", "a=10", "a=30"},
PlotRange -> {-1500, 150}, PlotLabel -> "L=1, c1=1"]
For values of a higher than 30, the shape stays the same, it only gets more extreme.
Now, if I try to approximate this rather well-behaved (apart from the boundaries) potential with my combination of two Gaussians, I still run into some problems:
L = 1;
L\[CurlyEpsilon] = (L - 10^-7*L)/2;
a = 30;
c1 = 1;
doubledistance[in1_?NumberQ, in2_?NumberQ, x1_?NumberQ,
x2_?NumberQ, \[Sigma]1_?NumberQ, \[Sigma]2_?NumberQ] :=
NIntegrate[(V4Sym[x] -
doubleGauss[x, in1, in2, x1,
x2, \[Sigma]1, \[Sigma]2])^2, {x, -L\[CurlyEpsilon],
L\[CurlyEpsilon]}]
{min, doublesol} = NMinimize[
doubledistance[in1, in2, x1, x2, \[Sigma]1, \[Sigma]2], {in1, in2,
x1, x2, \[Sigma]1, \[Sigma]2}]
yields the output
During evaluation of In[266]:= NIntegrate::inumr: The integrand (-0.281773 E^(-12.6195 x^2)-0.990791 E^(-8.36999 x^2)+V4Sym[x])^2 has evaluated to non-numerical values for all sampling points in the region with boundaries {{-(9999999/20000000),9999999/20000000}}. >>
During evaluation of In[266]:= NIntegrate::inumr: The integrand (-0.281773 E^(-12.6195 x^2)-0.990791 E^(-8.36999 x^2)+V4Sym[x])^2 has evaluated to non-numerical values for all sampling points in the region with boundaries {{-(9999999/20000000),9999999/20000000}}. >>
During evaluation of In[266]:= NIntegrate::inumr: The integrand (-0.281773 E^(-12.6195 x^2)-0.990791 E^(-8.36999 x^2)+V4Sym[x])^2 has evaluated to non-numerical values for all sampling points in the region with boundaries {{-(9999999/20000000),9999999/20000000}}. >>
During evaluation of In[266]:= General::stop: Further output of NIntegrate::inumr will be suppressed during this calculation. >>
During evaluation of In[266]:= NMinimize::nnum: The function value NIntegrate[(V4Sym[x]-doubleGauss[x,0.990791,0.281773,-0.345651,0.281501])^2,{x,-L[CurlyEpsilon],L[CurlyEpsilon]}] is not a number at {in1,in2,x1,x2,[Sigma]1,[Sigma]2} = {0.990791,0.281773,0.300512,0.930797,-0.345651,0.281501}. >>
During evaluation of In[266]:= NMinimize::nnum: The function value NIntegrate[(V4Sym[x]-doubleGauss[x,0.990791,0.281773,-0.345651,0.281501])^2,{x,-L[CurlyEpsilon],L[CurlyEpsilon]}] is not a number at {in1,in2,x1,x2,[Sigma]1,[Sigma]2} = {0.990791,0.281773,0.300512,0.930797,-0.345651,0.281501}. >>
During evaluation of In[266]:= NMinimize::nnum: The function value NIntegrate[(V4Sym[x]-doubleGauss[x,0.990791,0.281773,-0.345651,0.281501])^2,{x,-L[CurlyEpsilon],L[CurlyEpsilon]}] is not a number at {in1,in2,x1,x2,[Sigma]1,[Sigma]2} = {0.990791,0.281773,0.300512,0.930797,-0.345651,0.281501}. >>
During evaluation of In[266]:= General::stop: Further output of NMinimize::nnum will be suppressed during this calculation. >>
During evaluation of In[266]:= Set::shape: Lists {min,doublesol} and NMinimize[doubledistance[in1,in2,x1,x2,[Sigma]1,[Sigma]2],{in1,in2,x1,x2,[Sigma]1,[Sigma]2}] are not the same shape. >>
Out[266]= NMinimize[ doubledistance[in1, in2, x1, x2, [Sigma]1, [Sigma]2], {in1, in2, x1, x2, [Sigma]1, [Sigma]2}]
V4Sym
? $\endgroup$V4Sym
throughout and, apart from the boundaries at|x|=L/2
, where it has singularities, it is well-behaved. $\endgroup$