lsqcurvefit does not match data

Hi, I'm trying to fit measured data to a custom nonlinear function. However, the fit doesn't match the data, although there's some resemblance. I'm aware that the starting guess is quite important, but I don't really have "good" guess. I have attatched my code below.

% Armature displacement
xdata = linspace(-1.9e-3,1.9e-3,19);
% Magnetic saturation as a function of displacement
ydata = [5.87464279998885 4.73039540152481 4.67862244835332 3.93109398658849 ...
    3.07326104723833 2.82766772242832 2.14048524598602 1.69605944153972 1.68390931280163 1.64257093300466...
    1.68390931280163 1.69605944153972 2.14048524598602 2.82766772242832 3.07326104723833 3.93109398658849 ...
    4.67862244835332 4.73039540152481 5.87464279998885];

% Air gap length;
D = 1.9e-3;

%(edit) current
i = linspace(-0.0065,0.0065,19);

% Guess
x0 = [19.4614 -3.9480 0.6650];

sat = @(s,x) s(1).*(i + s(end).*(x)./D).^2 + s(2).*(i + s(end).*(x)./D).^4;
[s,resnorm,redidual] = lsqcurvefit(sat,x0, xdata, ydata);

plot(x/D,ydata)
hold on;
plot(x/D,sat(s,xdata))
title('Magnetic saturation vs. displacement')
xlabel('Normalized displacement, x/D')
ylabel('Saturation')

As you can see, the fitted curve deviates a bit from the measured curve. I have tried many different starting guesses, but the fitted curve always seems to deviate from the measured.

I think I might be missing something. Any help is much appreciated.

Cheers!

[Matt J]

What is 'i'?

[Sansi Rajalingam]

oops!

i = linspace(-0.0065,0.0065,19);

[David Goodmanson]

HI Sansi

I don't quite see what is going on here. ydata is a symmetric curve about zero, to 14 decimal places. So one half or the other of that curve cannot be experimental data. Two of your fitting functions, (x/D)^2 and (x/D)^4, are also symmetric about zero. But since the fitting function i is antisymmetric about zero, it can't contribute to the fit.

It is not recommended to use i as a variable, since in Matlab that use conflicts with i = sqrt(-1). If you use i as a variable, then later on in a program if you want to use i as the usual complex quantitiy, its value won't be right. There is not a conflict for people who always use 1i in place of i, but the less said about those kinds of people the better.

[Matt J]

There is not a conflict for people who always use 1i in place of i, but the less said about those kinds of people the better.

What??? How can you denigrate the use of 1i when the documentation clearly tells us that it is better???

[David Goodmanson]

Hi Matt,

It was easy, actually. But I don’t think that I denigrated 1i; I only denigrated everybody anywhere who uses it. In a more serious vein, although Mathworks recommends 1i, I prefer i because I think it simply looks a lot better and more compact in code, as well as being in accord with the notation used in mathematics and physics for several hundred years. I don't know why people might prefer 1i over i, but hey, if they do, go for it.

One consideration that could outweigh the extra keystroke and cumbersome look of 1i is the speed advantage mentioned in the tip. Something like

N = 1e8
x = 1:N
y = i*x
% or
y = exp(i*x)

shows no noticeable speed difference, whether you use i or 1i. On the other hand, with

tic
for k = 1:1e8
 z = 2*i;
end
toc

then using 1i is about 50 times faster. The result seems to hinge on how many times you have to query i or 1i to come up with a value. But in my experience it seems to be much more common to have an expression like i*(1:N) or i*rand(1,N) where i or 1i is only queried once.

Answer 1

[Matt J]

Your model is clearly wrong, assuming that the given data is valid. In your current model, saturation is always 0 when displacement is 0, but according to your ydata, saturation=2 at displacement=0. Adding a fourth offset parameter to the model greatly improves the fit,

x0=[nan,nan, nan, 0.6650]; %initial guess of nonlinear parameter only

 qdata=i + x0(end).*(xdata)./D;
 pol=polyfit(qdata,ydata,4);
 
x0(1:3)=pol([5,3,1]); %derive initial guesses from polyfit of the linear coefficients

sat = @(s,x) s(1) + s(2).*(i + s(end).*(x)./D).^2 + s(3).*(i + s(end).*(x)./D).^4;
[s,resnorm,redidual] = lsqcurvefit(sat,x0, xdata, ydata);

[Sansi Rajalingam]

Now that you point it out, I see that the model will never fit the data. I should probably have also mentioned that the measurements were only conducted for positive displacements, and the assumption was that the saturation would be symmetric for negative displacements. I now see that this assumption is not valid.

Many thanks.