Friday, July 20, 2012

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. AthertonFerromagnetic 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. BlundellMagnetism 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! ;) 

15 comments:

Fredy Waller said...

thank you very much
ans my question it s
How can fitting Jiles Atherton with ode45 ?

Fredy Waller said...

plz

polt hysteresis lopp from Jiles Atherton Model

eddie said...

ode23, ode45 etc are solvers for systems of ordinary differential equations! In this post I was just interested in showing how to cast an implicit equation into an explicit form. No ODE here involved.

eddie said...

Give me (better via email) the details of the system you want to solve.

Anonymous said...

Hi! I'd like to plot hysteresis, as well. I'm interested in modeling a ferromagnetic material with a coercivity ~ 100 Oe. Do you have any tips on how to go about this? Thanks!

eddie said...

These references should be of some help :)

http://www.sciencedirect.com/science/article/pii/0304885386900661

http://ieeexplore.ieee.org/xpl/login.jsp?tp=&arnumber=5114077&url=http%3A%2F%2Fieeexplore.ieee.org%2Fxpls%2Fabs_all.jsp%3Farnumber%3D5114077

http://ieeexplore.ieee.org/xpl/login.jsp?tp=&arnumber=1062594&url=http%3A%2F%2Fieeexplore.ieee.org%2Fxpls%2Fabs_all.jsp%3Farnumber%3D1062594

Anonymous said...

Thanks for the reply! I'll try these--not having much success with MATLAB; hopefully Mathematica makes it simpler.....

Kush mudgal said...

Sir, i want to model the SMA actuator through jiles atherton modelling. plzzzz.. help me.. i am unable to do it. I am working on a project.

eddie said...

Seriously, I don't have much time to work on my hobbys this period, otherwise I would have been glad to help you.

Kush mudgal said...

:(

Anonymous said...

Hi, Thanks for this informative post. Do you have a typo with the parenthesis?
(@(y) y - A*(coth(B*x(i)+C*y) - 1./(B*x(i)+C*y))

i.e. the coth operator - shouldn't it apply to the larger expression? as in

(@(y) y - A*(coth(B*x(i)+C*y - 1./(B*x(i)+C*y)))

Thanks,
Deepak

eddie said...

Sorry for having been late!
Anyway, you were right, there was a typo in the post image showing the implicit function, not in the matlab code!
Thanks for your comment!

Yusufcan Uz said...

eddie I wanna ask something about implicitly functions. can you give your email adress?

eddie said...

for the email address...
take this:
learn-one-thing-a-day
remove the dashes
and add the gmail DOT com

PS
I am very busy with a new job, so I don't have much time for my blog anymore :(

ms86 said...

Hey Eddie

Hope you're still not too busy at your job (as you posted in your previous comment)

I'm trying to fit magnetization of superparamagnetic particles with something similar to yours, except my model involves a sum or an integration of a size distribution times the langevin function. I'm really not certain how to do it.. :S

Might you have any suggestions?