Fitting real and imaginary parts of a function

Question

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

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

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

Answer 1

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