1
$\begingroup$

I am very new to fitting functions in mathematica and I have to fit an improper integral to data. I can import .csv file, but I cannot make a model.

Here are the definitions of functions:

e[t_] := E^(-t^2/w^2)*E^(-I*t*o);
eS[t_] := Conjugate[e[t]];
inten[t_] = e[t]*eS[t]

And I have a model which is a sum of improper integrals - here are some of them. Where is my mistake? Is this the correct way to write models?

model = {
    4*NIntegrate[inten[u]*inten[u - t], {u, -\[Infinity], \[Infinity]}] + 
     4*NIntegrate[(inten[u] + inten[u - t])*
        Re[e[u]*eS[u - t]], {u, -\[Infinity], \[Infinity]}]}
$\endgroup$
2
$\begingroup$

Introduction

I copied your equations as presented by yohbs

e[t_] := Exp[-t^2/w^2 - I*t*o];
eS[t_] = Simplify@ComplexExpand@Conjugate[e[t]];
inten[t_] = Simplify@ComplexExpand[e[t]*eS[t]];
integrand = 
  Simplify@ComplexExpand[
    4*inten[u]*inten[u - t] + 2*inten[u]^2 + 
     2*Re[e[u]*e[u]*eS[u - t]*eS[u - t]] + 
     4*(inten[u] + inten[u - t])*Re[e[u]*eS[u - t]]];

and performed the integration

model = Integrate[integrand, {u, -∞, ∞}]

which resulted in

E^(-(t^2/w^2)) Sqrt[π]
  Abs[w] (2 + E^(t^2/w^2) + 4 E^(t^2/(4 w^2)) Cos[o t] + Cos[2 o t])

an equivalent expression is

Sqrt[π]Abs[w] (1 + (4 Cos[o t])/E^((3 t^2)/(4 w^2)) +
                   (2 + Cos[2 ot])/E^(t^2/w^2))

I imported your data from your csv file. Below is a plot of the data.

ListPlot[data, PlotRange -> {Automatic, {0, 10}}, PlotStyle -> Red]

Mathematica graphics

A careful look shows that the data is not quite symmetrical about zero.

The edges asymptote to a value that is approximately 1.2 (i.e., when w get's large).

An expanded view of the data close to the origin shows the cyclic nature.

ListLinePlot[data, PlotRange -> {{-20, 20}, {0, 10}}, PlotStyle -> Red]

Mathematica graphics

By counting the peaks observe that the range -14 to -2 represents five periods so the approximate frequency of this wave is:

Solve[ω (14 - 2) == 10. π, ω]

(* {{ω -> 2.61799}} *)

The maximum value of the data is about 8.

The model evaluated at zero becomes

model /. t -> 0

(* 8 Sqrt[π] Abs[w] *)

Which implies that w is approximately

Solve[8 Sqrt[π] w == 8., w]

(* {{w -> 0.56419}} *)

However, this value don't fit at all. w also controls the envelope around the data.

Plot[model /. {w -> 0.56, o -> 2.6}, {t, -10, 10}, 
 PlotStyle -> Black, PlotRange -> {All, Full}]

Mathematica graphics

The envelope is so thin that one can't see the cyclic data.

To get the correct envelope shape we need w to be approximately 80.

Below I have plotted the model data (two scales) with data adjusted so that the magnitude matches.

Module[
 {
  w = 80,
  o = 2.6,
  max,
  data2 = data
  },
 max = 8 Sqrt[\[Pi]] Abs[w];
 data2[[All, 2]] = data2[[All, 2]]*max/8;
 Column[{
   Show[
    Plot[E^(-(t^2/w^2)) Sqrt[\[Pi]]
       Abs[w] (2 + E^(t^2/w^2) + 4 E^(t^2/(4 w^2)) Cos[o t] + 
        Cos[2 o t]) ,
      {t, -150, 150},
     PlotStyle -> Black,
      PlotRange -> {Automatic, Full},
     ImageSize -> 400
     ],
    ListPlot[data2, PlotStyle -> Red,
      PlotRange -> {{-20, 20}, Full}]
    ],
   Show[
    Plot[E^(-(t^2/w^2)) Sqrt[\[Pi]]
       Abs[w] (2 + E^(t^2/w^2) + 4 E^(t^2/(4 w^2)) Cos[o t] + 
        Cos[2 o t]) ,
      {t, -150, 150},
     PlotStyle -> Black,
      PlotRange -> {{-20, 20}, Full},
     ImageSize -> 400
     ],
    ListLinePlot[data2, PlotStyle -> Red,
      PlotRange -> {{-20, 20}, Full}]
    ]
   }]
 ]

Mathematica graphics

Note that when we get the frequency approximately correct, o ≈ 2.6, there is a phase shift between the data and the model.

The amplitudes are off by a factor of around 140 and there appears to be a slight skew to the envelope (i.e, not symmetrical about the origin)

New Model

In the new model we remove the term Sqrt[π] Abs[w] and introduce a new parameter, scale. We introduce a second new parameter, offset to account for the asymptote values.

Finally a time shift parameter, ts is used to account for the asymmetry about the origin and the phase shift.

The new model is:

model2 = scale ( (4 Cos[o (t - ts)])/E^((3 (t - ts)^2)/(4 w^2)) +
                 (2 + Cos[2 o (t - ts)])/E^((t - ts)^2/w^2) )   + offset

or

Mathematica graphics

NonlinearModelFit is used fit the data (I supplied starting values from experimentation).

fit = NonlinearModelFit[data, model2,
           {{o, 2.6}, {offset, 1.15}, {scale, 1}, {ts, 11.5}, {w, 80}}, t]

fit["BestFitParameters"]

(* {o -> 2.60488, offset -> 1.27066, scale -> 0.977271, 
 ts -> 11.3441, w -> 78.7623} *)

Below is a plot of the results.

Column[{
  Show[
   Plot[fit[t] ,
     {t, -345, 345},
    PlotStyle -> Black,
     PlotRange -> {Automatic, Full},
    ImageSize -> 400
    ],
   ListPlot[data, PlotStyle -> Red,
     PlotRange -> {{-20, 20}, Full}]
   ],
  Show[
   Plot[fit[t],
     {t, -150, 150},
    PlotStyle -> Black,
     PlotRange -> {{-20, 20}, Full},
    ImageSize -> 400
    ],
   ListLinePlot[data, PlotStyle -> Red,
     PlotRange -> {{-20, 20}, Full}]
   ]
  }]

Mathematica graphics

The new model fits the data much better than the original model.

$\endgroup$
  • $\begingroup$ @10001000101111 See if the new model makes sense to you. Also you can look at examples for NonlinearModelFit. $\endgroup$ – Jack LaVigne Jun 14 '17 at 2:14
  • $\begingroup$ This it a great answer - thank you! I thought that the difficulty in fitting the model was caused by the lack of pre-exponetial factor and I think that your scale parameter fixes that. I've noticed the shift (asymmetry) and I tried to accommodate for it while fitting in Origin Lab (which I'm more familiar with). But your solution works great $\endgroup$ – 10001000101111 Jun 14 '17 at 11:44
2
$\begingroup$

1) Mathematica always tries to get the most generic solutions. However, often these are too generic. You should probably tell mathematica that t,w,o are real numbers, if they are. In that case inten has a very simple form, which is a Gaussian. For this you can use ComplexExpand[inten[t]], or set $Assumptions accordingly.

2) The integral can be done analytically, so there's no need to use NIntegrate.

Altogether:

$Assumptions = Element[w | t | o, Reals];
e[t_] := Exp[-t^2/w^2 - I*t*o];
eS[t_] = Simplify @ ComplexExpand @ Conjugate[e[t]];
inten[t_] = Simplify @ ComplexExpand[e[t]*eS[t]];
integrand = 
 4 Simplify @ ComplexExpand[
    inten[u]*inten[u - t] + (inten[u] + inten[u - t])*
      Re[e[u]*eS[u - t]]
    ]
model = Integrate[integrand, {u, -\[Infinity], \[Infinity]}]

Edit

As I wrote in the comments, the model is not a good model for the data and the fitting is very difficult due to the rapid oscillations, but it still gives something:

data = Import["https://fs05n3.sendspace.com/dl/4ff87464123001f76123bfd144d26228/593fff8f3d88ae41/sxu5n8/data.csv"]
fit = NonlinearModelFit[data, model, {w, o}, t]
Plot[fit[x],{x,-300,300}]

enter image description here

$\endgroup$
  • $\begingroup$ Thanks for answering! I have two questions: 1) I thought that NIntegrate is faster, isn't that true? 2) I'm trying to fit this model to data: (fit = NonlinearModelFit[data, model, {w, o}, t]), but it gives me an error "earch specification 97 without variables should be a list with 1 to 4 \ elements". What is going on? $\endgroup$ – 10001000101111 Jun 12 '17 at 20:09
  • $\begingroup$ (a) NIntegrate might be faster or slower, depends on what you're integrating. Here I'd use Integrate because you only need to preform the integration once in the beginning of the calculation and have an analytic result which is valid for all w, t, etc that you don't have to recalculate at each iteration. (b) you will need to provide more code and the data you want to fit to in order for us to help you with that. $\endgroup$ – yohbs Jun 12 '17 at 21:43
  • $\begingroup$ sendspace.com/file/uecg1r Here is the code. I thought that perhaps I should add another parameter in front of the gauss for fitting, but first it would have to give something other than error. $\endgroup$ – 10001000101111 Jun 12 '17 at 22:00
  • $\begingroup$ There is a problem in the notebook with the data. One is not able to get all of it. There is a printed output with a subset of the data showing. You need to take other measures to get the data to users. Maybe save it to an <xxx>.mx file and share that link. $\endgroup$ – Jack LaVigne Jun 13 '17 at 0:24
  • $\begingroup$ I would advise against that. Try to simplify your problem as much as possible by creating a minimal working example that shows the problem. The actual original data is of little use since this is probably a simple syntax problem $\endgroup$ – yohbs Jun 13 '17 at 0:28

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.