4
$\begingroup$

I have an analytic form of a function 

G[w_] = -I*z0*w/(1 + I*z0*w*a)

I know this function can be obtained (up to some overall normalization) from solving a differential equation numerically for small w 

g[z_] := 1 + 3 (z/z0)^4 - 4 (z/z0)^3 DE = x''[z] + D[g[z]/z^2, z]/(g[z]/z^2) x'[z] + (w^2)/(g[z])^2 x[z] == 0; pfun = Block[{z0 = 0.1, eps = 10^-2, zb = 10^-2}, ParametricNDSolveValue[{DE, x'[z0 - eps] == I (z0^2 w)/(6 eps^2) x[z0 - eps], x[zb] == 1}, - (x'[zb]/zb^2), {z, zb, z0 - eps}, {w} ] ]

But G[w] and pfun[w] are two complex functions. To fit them (up to over all normalization) I need to fit the real and the imaginary parts separately. To fix the overall normalization once and for all I compute the ratio for very small w(the normalization doesn't depend on a for very small w as one can see from the expression of G[w]). 

test[a_] := Evaluate[Im[pfun[0.0005]]/Im[-((I w *0.1 )/(1 + I w *0.1* a)) //. w -> 0.0005]]

test[a] is approximately 1080 for large range of a.

Now if I try to match the real and imaginary parts independently the value of a comes out to be quite different in two fits!

list1 = Table[{w, Evaluate[Re[pfun[w]]]}, {w, -2, 2, 0.1}] aValue1 = FindFit[list1, Re[-1080.13 (I*w*0.1 )/(1 + I*w*0.1*a)], a, w]

This gives a=9.49

L1 = ListPlot[list1]; Show[Plot[Re[ -1080.13 (I w*0.1 )/(1 + I w*0.1*a)] //. aValue1, {w, -5, 5}], L1]]

Real parts of <code>G[w]</code>

Where as fitting the imaginary part in similar fashion I get a different value of a.

list2 = Table[{w, Evaluate[Im[pfun[w]]]}, {w, -2, 2, 0.2}]; aValue2 = FindFit[list2, Im[-1080.13 (I*w*0.1)/(1 + I*w*0.1*a)], a, w]

This gives a=1.35

L2 = ListPlot[list2]; Show[Plot[Im[ -1080.13 (I w*0.1 )/(1 + I w*0.1*a)] //. aValue2, {w, -20, 20}], L2]]

Imaginary parts of <code>G[w]</code>

My question is what should be the best fit value of a?

Improve this question
$\endgroup$
3
  • $\begingroup$ I edited your code to fix the code-formatting a bit (use four-space indentations for non-inline code blocks instead of the back-ticks), but there were some funny formatting symbols that showed up, so I'm not entirely sure I didn't delete something important. Please check to make sure it looks correct! $\endgroup$ Commented Dec 16, 2015 at 21:06
  • $\begingroup$ @march Sorry about the format and thanks for editing! :) I'll check.. $\endgroup$ Commented Dec 16, 2015 at 21:08
  • 1
    $\begingroup$ See this post and others it links to for simultaneous fitting (overlapping or same parameters, different data sets). $\endgroup$ Commented Feb 22, 2016 at 17:33

1 Answer 1

Reset to default
3
$\begingroup$

You don't need to fit the imaginary and real parts of your data separately; indeed, fitting them together should give you more accurate results.

FindFit will happily deal with complex-valued functions of real paramters. You just need to give it assumptions on the values of those parameters, and to select an appropriate norm function to use in the minimization process.

For instance, with your definitions:

g[z_] := 1 + 3 (z/z0)^4 - 4 (z/z0)^3 DE = x''[z] + D[g[z]/z^2, z]/(g[z]/z^2) x'[z]8 + (w^2)/(g[z])^2 x[z] == 0; pfun = Block[ {z0 = 0.1, eps = 10^-2, zb = 10^-2}, ParametricNDSolveValue[ {DE, x'[z0 - eps] == I (z0^2 w)/(6 eps^2) x[z0 - eps], x[zb] == 1}, -(x'[zb]/zb^2), {z, zb, z0 - eps}, w ] ];

we can calculate a list of complex values to fit, for small $w$ values:

list = Table[{w, pfun[w]}, {w, -1, 1, 0.05}];

The parametric function throws a few warnings on boundary conditions so you should check for the consistency of your model, but nonetheless it returns a list of values we can work with.

We then set up a model of the form you specified above (G[w]) in your original code:

model = -const I z0 w/(1 + I z0 w a) /. z0 -> 1/10;

Rather than using an estimated multiplicative constant in the model, we include one in the model to determine through fitting (const); I also included the value of $z_0=0.1$, so the free parameters are const and a.

We then set up the fitting function. In addition to a standard fitting procedure, we add the following:

  1. specify that a and const are real-valued;
  2. select the 1-norm (Norm[#, 1]& which corresponds to Total[Abs[#]]&) as the norm function;
  3. give reasonable starting values for the parameters.

This results in: 

fit = FindFit[ list, {model, {a, const} ∈ Reals}, {{a, 10}, {const, 1000}}, w, NormFunction -> (Norm[#, 1] &) ] (* Out: {a -> 8.92677, const -> 1529.84} *)

Here is a plot of the real and imaginary parts of list, together with the corresponding fitted functions:

Show[ ListPlot[Transpose@ReIm[list[[All, 2]]], DataRange -> MinMax@list[[All, 1]]], Plot[Evaluate@ReIm[model /. fit], Evaluate@Flatten@{w, MinMax@list[[All, 1]]}] ]

results of fit

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.