How to fit IMPLICIT functions with Matlab! A worked-out example

In Matlab you can perform implicit curve fit using either "nlinfit" or "lsqcurvefit" (the last one only if you have the optimization toolbox) but you need to re-write the implicit function into an explicit one solving it numerically (i.e. using "fzero").
THERE IS NO WAY to fit directly using an implicit expression!

To learn how to successfully fit using an implicit function you can follow the instructions on this Matlab page page or to read and use my worked-out example.

I will show you how my implementation of the Langevin function with a field Weiss (that is an implicit equation). This example is interesting by itself becasue it addresses the topic of how to model an hysteresis loop.
The equation is taken from the paper by D.C. Jiles and D.L. Atherton,  Ferromagnetic Hysteresis, IEEE Transaction on Magnetics, Vol. Mag-19, No. 5, September 1983.

For a general reference about the Langevin equation (and generally for magnetism) I recommend the nice (and inexpensive!) book by S. Blundell,  Magnetism in Condensed Matter.

So let's start!

So we want to create a fitting function in Matlab using this IMPLICIT equation:

As you see the above equation is implicit due to the term alpha*M.

The first step is to transform the equation into an EXPLICIT one. As I wrote in the beginning, it is necessary just to use the fzero function.
I named this function  LangevinWithWeissField.m (if you are lazy, you can download the file .m directly from this link)

%---------------------------------------------------------------------------------------
function y=LangevinWithWeissField(p, x)% Implicit function: Langevin function with a Weiss field.
% y = A*(coth(B*x(i)+C*y) - 1./(B*x(i)+C*y))
% p is the parameter vector
% Physical meaning of the parameters
% x is the magnetic field H
% y is the magnetization M corresponding to n*<m>z where n is the number of
% magnetic moments,  and <m>z is the projected magnetic moment along the
% field direction (called z) of the total moment mu
%A = amplitude
%B = mu0 m/(kB*T) where kB is the Boltzman constant and T is the temperature
%C is a factor giving the Weiss field (that is proportional to the total
%magnetization y, i.e. M=n*<m>z, note that this term includes also the previous factor
% mu0 m /(kB*T)

% assign the parameters ...

A=p(1);
B=p(2);
C=p(3);

y=zeros(size(x)); % define a vector to allocate the magnetization values
NN=length(x); % total length of the field vector x, i.e. B

opt = optimset('display','off');
% I out off all the messages coming from fzero. If something goes wrong, change this option to 'off'
% to see at which x values fzero failed.

for i=1:NN   
   y(i)=fsolve(@(y) y - A*(coth(B*x(i)+C*y) - 1./(B*x(i)+C*y)) , 0.0001, opt);
   % Here 0.0001 is our starting point to find the solution around 0.
end

end % close the function
%---------------------------------------------------------------------------------------


Once you have saved the code in your Matlab working directory, all the work is done!

You can plot the Langevin function with the parameter you like, for example (from the Matlab command window)

xx=linspace(-30, 30, 200); % magnetic field
p=[10, 2, 0.1]; %parameters
yy=LangevinWithWeissField(p, xx); % magnetization values
plot(xx, yy, '*-r'); % plot

Now let's see how to "fit" the y-data we just created, i.e. the vector yy.
I use the  LSQCURVEFIT function as done in this mathworks' link.
The code to use in the command line is

lsqcurvefit(@(params, xdata) LangevinWithWeissField(params, xdata),[1 1 1], xx, yy)

In the previous line, the vector  [1 1 1] represents the initial guess of the fitting parameters.
After a while (on my computer it took 54 seconds) the lsqcurvefit gives you the optimal fit parameters, that corresponds to the original values stored in the p vector, namely [10 2 0.1].


To conclude, here is two things to have to keep in mind.
1) It's VERY important to use a good set of values as initial guess for your fit.  With bad (=random) combination of values, your fit could not converge!

2) Using implicit functions can very easily make the fitting SLOW, terribly SLOW!
So, try to avoid their use whenever it is possible.

Any question or comment, as always, is welcome!


UPDATE
You could also be interested in these other posts

(A)  How to hold fitting parameters in nlinfit (MATLAB) 
(B) A simple matlab fitting interface
 


I hope that you liked this post and if you did, consider to support my blog! ;)