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

Improve this question
$\endgroup$
7
  • $\begingroup$ What is the definition of V4Sym? $\endgroup$
    – Bob Hanlon
    Mar 14, 2018 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, 2018 at 13:52
  • $\begingroup$ Edit the question to include the definition. $\endgroup$
    – Bob Hanlon
    Mar 14, 2018 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, 2018 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, 2018 at 20:59

1 Answer 1

Reset to default
4
$\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}}
Improve this answer
$\endgroup$
0

Your Answer

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

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