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:

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); % 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) A simple matlab fitting interface

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

## 15 comments:

thank you very much

ans my question it s

How can fitting Jiles Atherton with ode45 ?

plz

polt hysteresis lopp from Jiles Atherton Model

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.

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

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!

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

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

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.

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

:(

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

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!

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

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 :(

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?

Post a Comment