1
$\begingroup$

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

Shape of V4Sym for a=1 (blue), a=5 (orange), a=10 (green) and a=30 (red)

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

$\endgroup$
  • $\begingroup$ What is the definition of V4Sym? $\endgroup$ – Bob Hanlon Mar 14 '18 at 12:45
  • $\begingroup$ Sorry, forgot about that: 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) $\endgroup$ – GeorgF Mar 14 '18 at 13:52
  • $\begingroup$ Edit the question to include the definition. $\endgroup$ – Bob Hanlon Mar 14 '18 at 13:55
  • $\begingroup$ I don't think that's the problem, I've made several plots of V4Sym throughout and, apart from the boundaries at |x|=L/2, where it has singularities, it is well-behaved. $\endgroup$ – GeorgF Mar 14 '18 at 19:31
  • 1
    $\begingroup$ I would adivse you to try using sums rather than integrals here. There’s no point is calculating them for such a simple minimization. $\endgroup$ – Vsevolod A. Mar 14 '18 at 20:59
3
$\begingroup$

Make the constants using Set (=) rather than SetDelayed. Make distance an explicit function of its parameters, and only defined when they have numeric values. Use NIntegrate to avoid symbolic integrals that might get into trouble with exponents.

doubleGauss[x_, in1_, in2_, x1_, 
  x2_, \[Sigma]1_, \[Sigma]2_] := 
 in1*Exp[-(x - x1)^2/\[Sigma]1^2] + in2*Exp[-(x - x2)^2/\[Sigma]2^2]

L = 10^-5;
L\[CurlyEpsilon] = (L - 10^-8*L)/2;
a = 10^8;
c1 = 1;

V4Sym[x_] := 
 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[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]}]

In[9]:= {min, sol} = 
 NMinimize[
  distance[in1, in2, x1, x2, \[Sigma]1, \[Sigma]2], {in1, in2, x1, 
   x2, \[Sigma]1, \[Sigma]2}]

Out[9]= {1.05066267514*10^-13, {in1 -> 0.980323726315, 
  in2 -> 1.00030736963, x1 -> -1.67814472669, 
  x2 -> 0.0197762596417, \[Sigma]1 -> -0.397188653907, \[Sigma]2 -> 
   0.976862432515}}
$\endgroup$

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.