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 it's done!

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

**xx=linspace(-30, 30, 200); % m****a**gnetic 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)

edufit: a Matlab data fitting interface
I hope that you liked this post and if you did, consider to support my blog! ;)