# Modelling Hysteresis with a Differential Equation

I want to implement the bulk ferromagnetic hysteresis model (mostly the Jiles-Atherton Model), see http://drum.lib.umd.edu/bitstream/1903/6043/1/PhD_99-1.pdf page 44 equation (30).

The needed equations are $$\delta = \text{sgn}(\dot{H})\\[1ex] H_\text{eff}=H+\alpha M \\[1ex] M_\text{an}(H_\text{eff})= M_\text{sat}\cdot\mathcal{L}(H_\text{eff}/a)\\[1ex] \delta_\text{M}=\begin{cases}0 & \dot{H}<0\text{ and }M_\text{an}-M>0\\0 & \dot{H}>0\text{ and }M_\text{an}-M<0\\1&\text{else}\end{cases}$$ and the main equation is $$\frac{\text{d}M}{\text{d}H}=\frac{\frac{k\delta c}{\mu_0}\frac{\text{d}M_\text{an}}{\text{d}H_\text{eff}}+\delta_\text{M}(M_\text{an}-M)}{\frac{k\delta}{\mu_0}-\delta_\text{M}(M_\text{an}-M)\alpha-\frac{k\delta \alpha c}{\mu_0}\frac{\text{d}M_\text{an}}{\text{d}H_\text{eff}}}$$

Short explanation: $\delta$ takes into account wheter we are at the upper hysteresis branch or the lower one, $H_\text{eff}$ is needed to shift the Langevinfunction $\mathcal{L}(x)$ to achieve a hysteresis, $M_\text{an}$ is the anhysteretic Magnetization and $\delta_\text{M}$ is necessary to prevent unphysical behaviour.

## How I tried to implement it

There were two approaches I tried: The first one is based on ParametricNDSolve, and, as I understand it, the values need to go from $0$ to $H_\text{max}$ (this is no problem and works) then down to $-H_\text{max}$ and up again to $H_\text{max}$ (this is the part I'm unsure how to achieve).

My first guess was to define the function piecewise, with an implemented shift (Beware monster ugly function below!):

µ0 = 4 π*10^-7;
BulkFerromagneticHysteresisModelComplete = ParametricNDSolveValue[{Mges'[Happlied] ==
Piecewise[{{((k*c/µ0)*
Msat (a/(Happlied + α*Mges[Happlied])^2 -
Csch[(Happlied + α*Mges[Happlied])/a]^2/a) +
Piecewise[{{0,
Msat*Piecewise[{{#/3,
Abs[#] < 10^-5}, {Coth[#] - 1/(#),
True}}] &@((Happlied + α*Mges[Happlied])/
a) - Mges[Happlied] > 0}, {1,
True}}]*(Msat*
Piecewise[{{#/3, Abs[#] < 10^-5}, {Coth[#] - 1/(#),
True}}] &@((Happlied + α*Mges[Happlied])/
a) - Mges[Happlied]))/((k/µ0) -
Piecewise[{{0,
Msat*Piecewise[{{#/3,
Abs[#] < 10^-5}, {Coth[#] - 1/(#),
True}}] &@((Happlied + α*Mges[Happlied])/
a) - Mges[Happlied] > 0}, {1,
True}}]*(Msat*
Piecewise[{{#/3, Abs[#] < 10^-5}, {Coth[#] - 1/(#),
True}}] &@((Happlied + α*Mges[Happlied])/
a) - Mges[Happlied])*α - (k*
c*α/µ0)*
Msat (a/(Happlied + α*Mges[Happlied])^2 -
Csch[(Happlied + α*Mges[Happlied])/a]^2/a)),
Happlied <
20000}, {((-k*c/µ0)*
Msat (a/(-Happlied +
40000 + α*Mges[-Happlied + 40000])^2 -
Csch[(-Happlied +
40000 + α*Mges[-Happlied + 40000])/a]^2/a) +
Piecewise[{{0,
Msat*Piecewise[{{#/3,
Abs[#] < 10^-5}, {Coth[#] - 1/(#),
True}}] &@((-Happlied +
40000 + α*Mges[-Happlied + 40000])/a) -
Mges[-Happlied + 40000] > 0}, {1,
True}}]*(Msat*
Piecewise[{{#/3, Abs[#] < 10^-5}, {Coth[#] - 1/(#),
True}}] &@((-Happlied +
40000 + α*Mges[-Happlied + 40000])/a) -
Mges[-Happlied + 40000]))/((-k/µ0) -
Piecewise[{{0,
Msat*Piecewise[{{#/3,
Abs[#] < 10^-5}, {Coth[#] - 1/(#),
True}}] &@((-Happlied +
40000 + α*Mges[-Happlied + 40000])/a) -
Mges[-Happlied + 40000] > 0}, {1,
True}}]*(Msat*
Piecewise[{{#/3, Abs[#] < 10^-5}, {Coth[#] - 1/(#),
True}}] &@((-Happlied +
40000 + α*Mges[-Happlied + 40000])/a) -
Mges[-Happlied + 40000])*α - (-k*
c*α/µ0)*

Msat (a/(-Happlied +
40000 + α*Mges[-Happlied + 40000])^2 -
Csch[(-Happlied +
40000 + α*Mges[-Happlied + 40000])/a]^2/a)),
20000 < Happlied < 60000},{((k*c/µ0)*
Msat (a/(Happlied +
40000 + α*Mges[Happlied - 80000])^2 -
Csch[(Happlied +
40000 + α*Mges[Happlied - 80000])/a]^2/a) +
Piecewise[{{0,
Msat*Piecewise[{{#/3,
Abs[#] < 10^-5}, {Coth[#] - 1/(#),
True}}] &@((Happlied -80000 + α*Mges[Happlied -80000])/a) -
Mges[Happlied -80000] > 0}, {1,
True}}]*(Msat*
Piecewise[{{#/3, Abs[#] < 10^-5}, {Coth[#] - 1/(#),
True}}] &@((Happlied -80000 + α*Mges[Happlied -80000])/a) -
Mges[Happlied -80000]))/((k/µ0) -
Piecewise[{{0,
Msat*Piecewise[{{#/3,
Abs[#] < 10^-5}, {Coth[#] - 1/(#),
True}}] &@((Happlied -
80000 + α*Mges[Happlied - 80000])/a) -
Mges[Happlied - 80000] > 0}, {1,
True}}]*(Msat*
Piecewise[{{#/3, Abs[#] < 10^-5}, {Coth[#] - 1/(#),
True}}] &@((Happlied -
80000 + α*Mges[Happlied -80000])/a) -
Mges[Happlied -80000])*α - (k*
c*α/µ0)*

Msat (a/(Happlied -
80000 + α*Mges[Happlied - 80000])^2 -
Csch[(Happlied -
80000 + α*Mges[Happlied -80000])/a]^2/a)),
60000 < Happlied < 100000}}],
Mges == 0.0001}, {Mges}, {Happlied, 0, 100000}, {Msat, a,
k, α, c}];


Unfortunately this doesn't work (I suspect because of Mges[Happlied-80000] respectively Mges[-Happlied+40000]? If yes, how can I workauround this?).

EDIT: To be more precise, compiling the above shown works, proceeding with

modelcomplete[Msat_, a_, k_, α_, c_][Happlied_] := Through[BulkFerromagneticHysteresisModelComplete[Msat, a,
k, α, c][Happlied], List] /; And @@ NumericQ /@ {Msat, a, k, α, c, Happlied};


leads after several seconds (it looks as if it is already done for some time with compiling) to the error:

Mathematica Kernel for Windows has stopped working
Problem Event Name: APPCRASH
Application Name: MathKernel.exe
...
Fault Module Name: mathdll.dll
...


Another possible approach is described in Computer-aided Analysis of Electric Machines: A Mathematica Approach from Vlado Ostovic, but the book is from the early 90s, so the code presented there doesn't work - the main idea is to use a periodic function to achieve the hyteresis effect. It seems I got the virgin curve working there as well:

Ms = 2*10^6;
a = 1100;
kon = 0.0004;
alpha = 1.6*10^3;
mu0 = N[4*Pi*10^(-7)];
c = 0.2;
Hmax = 8000;
omega = N[2*Pi];

eq = {Mges'[t] == (Ms (Coth[(Hmax*Cos[2 Pi t] + alpha Mges[t])/a] -
a/(Hmax*Cos[2 Pi t] + alpha Mges[t])) -
Mges[t])/((1 + c) (kon/mu0 -
alpha (Ms (Coth[(t + alpha Mges[t])/a] -
a/(Hmax*Cos[2 Pi t] + alpha Mges[t])) - Mges[t]))) +
c/(1 + c) (Ms (a/(Hmax*Cos[2 Pi t] + alpha Mges[t])^2 -
Csch[(Hmax*Cos[2 Pi t] + alpha Mges[t])/a]^2/a)), Mges == 0}

tests = NDSolve[eq, {Mges}, {t, 0, 1}]

Plot[Evaluate[Mges[t] /. tests], {t, 0, 1}, PlotRange -> All]


But instead of the Mges[t] it would be necessary to use Mges[Hmax*Cos[2 Pi t]], sadly Mathematica crashes if I try to use this expression.

EDIT: To be more precise here: If I use Mges[Hmax*Cos[2 Pi t]] throughout eq I get the following error:

NDSolve::derivs: No derivatives of dependent variables were found in the equations. NDSolve is designed to solve differential or differential algebraic equations. Use NSolve or FindRoot to numerically solve algebraic equations. >>


I suspect this is due to Mathematica not recognizing Mges[Hmax*Cos[2 Pi t]] as $M(x)$ but as $M(f(x))$.

Does anyone have any idea how to overcome these obstacles? Thanks a lot in advance, I really appreciate any help!

• John, could you please clarify what you mean by “doesn’t work” for the first method? And for the second method do you mean the kernel crashes? – Verbeia Apr 26 '14 at 16:56
• @Verbeia I tried to clarify the admittedly unclear formulation, see the two edits. Thanks for pointing out! – John Apr 26 '14 at 17:19