An improved parametric model of hysteresis loop
Rostislav V. Lapshin1, 2
1

Solid Nanotechnology Laboratory, Institute of Physical Problems, Zelenograd, Moscow, 124460, Russian

Federation
2

Department of Photosensitive Nano and Microsystems, Moscow Institute of Electronic Technology, Zelenograd,

Moscow, 124498, Russian Federation
E-mail: rlapshin@gmail.com
A number of improvements have been added to the existing analytical model of hysteresis loop defined in parametric form. In particular, three phase shifts are included in the model, which permits to tilt the hysteresis loop
smoothly by the required angle at the split point as well as to smoothly change the curvature of the loop. As a result, the error of approximation of a hysteresis loop by the improved model does not exceed 1%, which is several
times less than the error of the existing model. The improved model is capable of approximating most of the known
types of rate-independent symmetrical hysteresis loops encountered in the practice of physical measurements. The
model allows building smooth, piecewise-linear, hybrid, minor, mirror-reflected, inverse, double and triple loops.
One of the possible applications of the model developed is linearization of a probe microscope piezoscanner. The
improved model can be found useful for the tasks of simulation of scientific instruments that contain hysteresis
elements.
Keywords: hysteresis, hysteresis loop, parametric form, piecewise-linear hysteresis loop, play, non-ideal relay,
double hysteresis loop, triple hysteresis loop, area of hysteresis loop, piezoscanner linearization, scanning probe
microscope, SPM, AFM, STM, simulation

I. INTRODUCTION
The phenomenon of hysteresis is widespread in nature, it is often met in many fields of science and engineering including instruments used in scientific research.1, 2, 3 There are a number of quite complicated analytical models describing this phenomenon.4 One of the simpler ones is the analytical model suggested in Ref. 1. With that
model, a family of hysteresis loops is described by the following parametric equations
x (α ) = a cosm α +b x sin n α ,
y (α ) =b y sin α ,

(1)

where α is a real parameter (α=0…2π); a is x coordinate of the split point
(see Fig. 1); bx, by are the saturation point coordinates; m is an integer odd
number (m=1, 3, 5, …) defining the curvature of the hysteresis loop; n is
an integer number defining the type of the hysteresis loop and its curvature. With n=1, the Leaf loop type is formed; with n=2 – the Crescent
(Boomerang), and with n=3 – the Classical.
The key distinctive feature of the model (1) is its simplicity. The
model is intuitive, it allows quickly creating hysteresis loops of a required
type and easily determining their parameters a, bx, by, m, n. The model (1)
covers most of the known types of rate-independent symmetrical hystere1

sis loops. The 1.5-6% approximation accuracy of the model (1) is quite

Fig. 1. Hysteresis loops of Leaf (n=1),
Crescent (Boomerang, n=2), and Classical (n=3) types. All the three loops
have the same area S.

An improved parametric model of hysteresis loop
enough for most practical tasks. However, there are cases when a higher accuracy is required. The improved
model5 suggested in the article approximates hysteresis loops with an error 1% or less.
To simplify constructions, analysis, and identification of hysteresis loops, Mathcad® worksheets (MathSoft,
USA) are attached to the article (see Ref. 6) as supplementary materials. In order not to overcomplicate the article,
the derivation of some formulae is omitted. The detailed derivation can be found in the supplementary materials.

II. DESCRIPTION OF THE IMPROVED MODEL
The main modification to the previous hysteresis model (1) is the introduction of phase shifts ∆α1, ∆α2, ∆α3
x (α ) = aˆ cosm (α + ∆α1 ) + bˆx sin n (α + ∆α 2 ),
y (α ) = by sin(α + ∆α 3 ),

(2)

where â , b̂x are corrected parameters of a, bx, respectively.

A. Smooth hysteresis loops
First, let us consider the influence of each of the three phase shifts.

1. Loop tilting by the phase shift ∆α1
The phase shift ∆α1 allows tilting a hysteresis loop by changing slope
angle β =π/2-θ of a tangent to the loop at the split point a (see Fig. 2). In
the model suggested earlier, loop tilting was available by rotation of coordinate system by angle θ and predistortion of the loop parameters a, bx, by
Fig. 2. Tilting hysteresis loop of the
Classical type with the phase shift
∆α1. Tilting by angle ±θ at the split
point a is provided with a shift by
m ∆α1. Loop tilting results in an increase in the loop area S. The area of
the loop with a positive slope is equal
to the area of the loop with a negative
slope.

by their rotation in the opposite direction. Since the introduction of the
phase shift ∆α1 leads to a change of the coordinates of the points of splitting a and saturation bx, a correction of coordinates of these points is required. In the improved model (2), the corrected parameters â , b̂x are
found from the following simple system of equations composed for the
split point α=0 and the saturation point α=π/2 (∆α2=∆α3=0)

a cos m 0 + b x sin n 0 = aˆ cos m (0 + ∆α1 ) + bˆx sin n 0,
a cos m

π
2

+ b x sin n

π

π
π

= aˆ cos m  + ∆α1  + bˆx sin n ,
2
2
2


(3)

whence
aˆ =

a
,
cos ∆α1
m

(4)

bˆx =b x + a tan m ∆α1
can be easily determined. Leaf type loops (n=1) with m=1 do not depend on the phase shift ∆α1. The phase shift

∆α1 required for loop tilting by the preset angle θ at the split point a (α=0) is calculated by the formula6
∆α1 = − arctan

by tan θ
am

.

(5)

Unsplit loops (a=0) cannot be tilted with the phase shift ∆α1.

2. Changing loop curvature by the phase shift ∆α2
The phase shift ∆α2 allows for changing the curvature of a hysteresis loop (see Fig. 3). Unlike the parameter m,
2

R. V. Lapshin
the phase shift ∆α2 provides a
continuous change of the curvature of the loop. Like the phase
shift ∆α1, in case of the phase
shift ∆α2, corrected coordinates
â , b̂x of the split and the saturation points, respectively, should
be determined. To do that, a sys-

Fig. 3. Continuous change in the curvature of the hysteresis loop by the
phase shift ∆α2. The loop area S increases with an increase in the phase
shift ∆α2.

Fig. 4. Continuous change in the curvature of the hysteresis loop by the
phase shift ∆α3. The loop area S decreases with an increase in the phase
shift ∆α3.

tem of equations composed for
the split point α=0 and the saturation point α=π/2 should be solved
(∆α1=∆α3=0)

a cos m 0 + b x sin n 0 = aˆ cos m 0 + bˆx sin n (0 + ∆α 2 ),
a cos m

π
2

+ b x sin n

π
2

= aˆ cos m

π

π

+ bˆx sin n  + ∆α 2 ,
2
2


(6)

whence

aˆ = a −b x tan n ∆α 2 ,
bˆx =

bx
.
cosn ∆α 2

(7)

Leaf type loops (n=1) with m=1 do not depend on the phase shift ∆α2.

3. Changing loop curvature by the phase shift ∆α3
Like the phase shift ∆α2, the phase shift ∆α3 allows for continuous change of the curvature of a hysteresis loop
(see Fig. 4). Composing equations for the split point α=0 and the saturation point α=π/2, like it has been done
above, one can obtain (∆α1=∆α2=0)
a cos m 0 + b x sin n 0 = aˆ cos m (0 − ∆α 3 ) + bˆx sin n (0 − ∆α 3 ),
a cos m

π
2

+ b x sin n

π

π

π

= aˆ cos m  − ∆α 3  + bˆx sin n  − ∆α 3 ,
2
2

2


(8)

whence the corrected parameters can be found as follows
aˆ =

a cos n ∆α 3 + b x sin n ∆α 3
,
sin m + n ∆α 3 + cos m + n ∆α 3

b cos m ∆α 3 − a sin m ∆α 3
bˆx = x m + n
.
sin ∆α 3 + cos m + n ∆α 3

(9)

The shape of Leaf type loops (n=1) with m=1 does not depend on the phase shift ∆α3.
The effect of the phase shift ∆α3=∆α is opposite to the one of ∆α2=∆α (compare Fig. 4 with Fig. 3). Moreover,
the loop built with ∆α2=∆α is somewhat different from the loop built with ∆α3=-∆α. The difference grows as |∆α|
increases. The phase shift ∆α3 causes an offset of parameter α by value ∆α3. As a result, parameter α becomes
equal to -∆α3 at the split point a and to π/2-∆α3 at the saturation point b.

3

An improved parametric model of hysteresis loop
4. Using several phase shifts simultaneously
In the general case when all the three phase shifts ∆α1, ∆α2, and ∆α3 are used, the system of equations composed for the split point α=0 and for the saturation point α=π/2 will look as follows
a cos m 0 + b x sin n 0 = aˆ cos m (0 + ∆α1 − ∆α 3 ) + bˆx sin n (0 + ∆α 2 − ∆α 3 ),
a cos m

π
2

+ b x sin n

π

π

π

= aˆ cos m  + ∆α1 − ∆α 3  + bˆx sin n  + ∆α 2 − ∆α 3 .
2
2

2


(10)

Solving the system (10), the searched for corrected parameters are found

aˆ =

a cos n (∆α 2 − ∆α 3 ) −b x sin n (∆α 2 − ∆α 3 )
,
sin m (∆α1 − ∆α 3 ) sin n (∆α 2 − ∆α 3 ) + cos m (∆α1 − ∆α 3 ) cosn (∆α 2 − ∆α 3 )

bˆx =

a sin m (∆α1 − ∆α 3 ) + b x cos m (∆α1 − ∆α 3 )
.
sin m (∆α1 − ∆α 3 ) sin n (∆α 2 − ∆α 3 ) + cos m (∆α1 − ∆α 3 ) cos n (∆α 2 − ∆α 3 )

(11)

By substituting the obtained parameters â , b̂x into (2), it is easy to show that with m=1, loops of the Leaf type
(n=1) do not depend on the phase shifts ∆α1, ∆α2 and the shape of these loops does not depend on the phase shift
∆α3.
Although the effect of the phase shift ∆α3 is opposite to the effect of the phase shift ∆α2, they do not completely neutralize each other when used jointly ∆α2=∆α3=∆α (∆α1=0). When the two phase shifts are set equal to
each other ∆α2=∆α3=∆α, their joint action is opposite to the action of the phase shift ∆α1=∆α. Therefore, in order
to tilt a hysteresis loop in the split point, the phase shifts ∆α2=∆α3=-∆α can be used instead of the phase shift
∆α1=∆α. Since the phase shifts ∆α2 and ∆α3=-∆α2 produce similar effect, the practical application of the model (2)
may involve only two of them: ∆α1 and ∆α2 (see Section III) or ∆α1 and ∆α3.
It is noteworthy that with ∆α1=∆α2=∆α3=∆α≠0, where ∆α is an arbitrary real number, the corrected parameters â , b̂x of the improved model (2) degenerate into the parameters a, bx of the original model (1), respectively.
Thus, in the case under consideration, the hysteresis loops built by the model (2) coincide in shape with the loops
built by the model (1), but the parameter α of the former is shifted by ∆α as compared to the latter. The conclusion
about the degeneration of the model (2) into the model (1) also comes out of the following consideration. Since the
effect of the joint usage of the equal phase shifts ∆α2=∆α3=∆α is opposite to the effect of the phase shift ∆α1=∆α,
the joint usage of three identical phase shifts ∆α1=∆α2=∆α3=∆α should lead to the formation of loop (1).

5. Hysteresis loops of Bat and Astro types
In Ref. 1, the hysteresis loop of Bat (Butterfly) type was obtained by taking the absolute value of generating
function y(α), i. e., ~
y (α ) = b y sin α . A similar type of loop can also be obtained by raising sin(α) in y(α) to an even
power k=2, 4, …, i. e.,
~
x (α ) = x (α ),
~
y (α ) =b sin k α .

(12)

y

The Bat type hysteresis loop with k=2 is shown in Fig. 5(a). Odd powers k=3, 5, … will yield loops of Astro7 type.
The hysteresis loop of Astro type with k=3 is shown in Fig. 5(b).

6. Arithmetic operations on hysteresis loops
The hysteresis loops based on models (1) and (2) can be added, subtracted, multiplied by a number, raised to a

4

R. V. Lapshin
power, and used as arguments in some functions
x (α ) = ∑ Ai xiki (α ),
i

(13)

y (α ) = ∑ Bi yili (α ),
i

where xi(α), yi(α) are expressions (1) or (2) of the ith loop; Ai, Bi are real factors; ki, li are positive integer and usually odd powers. The arithmetic operations (13) defined on the loops allow changing the shape and curvature of
hysteresis loops as well as making double loops. Addition of loops of different types is admissible.

(a)

(b)

Fig. 5 Hysteresis loop (a) Bat (Butterfly), (b) Astro.
An example of the addition of two hysteresis loops of different types
– the Leaf loop (a1=a/2, b1x=0.27, b1y=by/2, m1=3, n1=1) and the Classical
loop (a2=a-a1, b2x=bx-b1x, b2y=by-b1y, m2=n2=3) is shown in Fig. 6. The
parameters a, bx, by, m, θ of the resulting loop are given in the top left

Fig. 6. Tilted Classical loop (green) as
a result of the addition of two hysteresis loops – the Leaf loop (red) and the
Classical loop (blue). The Classical
loop tilted by the same angle θ=15° at
the split point a by means of the phase
shift ∆α1 is shown for comparison in
the dotted line.

corner of the figure. It should be noted that by adding with a Leaf type
loop, Classical loops can be tilted at the split point a by angle θ (β =π/2-θ is a slope angle of tangent at point a). By
taking into account the preset angle θ, the loop parameter b1x (or b1y or b2y or by) can be found according to the formula6

b1 x = (b1 y + b2 y )tan θ = by tan θ .

(14)

B. Piecewise-linear hysteresis loops
Within the model under consideration, the simplest way of building
piecewise-linear loops is as follows. The period T=2π of parameter α is
divided into the required number k of intervals, where k is an even integer
(k≥4). Continuous values α in formula (2) are then replaced with values
changing with step T/k so obtaining the x, y coordinates of the points of a
piecewise-linear hysteresis loop. Those coordinates are then connected
with line segments. An example of piecewise-linear loops of three types
Leaf, Crescent (Boomerang), and Classical built for k=12 (six-linear loop,
k/2=6) is given in Fig. 7.

Fig. 7. Piecewise-linear
hysteresis
loops of Leaf (n=1), Crescent (Boomerang, n=2), and Classical (n=3)
types. The number of division intervals k=12 (six-linear loop).
5

An improved parametric model of hysteresis loop
1. Using trapezoidal pulses instead of sine and cosine functions
Another method of obtaining piecewise-linear hysteresis loops consists in replacing the sine and cosine functions in the generating functions x(α) and y(α) of the model (2) by trapezoidal pulses trps and trpc of unit amplitude,
respectively,
x (α ) = aˆ trp cm (α + ∆α 1 ) + bˆx trp sn (α + ∆α 2 ),
y (α ) = b y trp s (α + ∆α 3 ),

(15)

where trpc(α)=trps(α+T/4); T is the period of pulses. The expressions defining the trapezoidal pulses trp of unit
amplitude are given in Ref. 1. The simplest loops are obtained with the zero phase shifts ∆α1=∆α2=∆α3=0, when
the upper d and the lower D bases of the trapezoidal pulses relate as D=3d (T=d+D=4d). In this case, the shape of
loops (15) is independent of the parameter m.
The formulae for the corrected parameters â and b̂x of the piecewise-linear loops (15) are derived using the
same logic as in the case of smooth loops. The working formulae for those parameters can be obtained from formulae (4), (7), (9), (11) valid for smooth loops (2) by simply replacing the sine and the cosine functions with trapezoidal pulses trps and trpc, respectively.
Fig. 8 shows piecewise-linear loops of three types: Leaf (also known
as two-linear loop or Play without “Whiskers”), Crescent (Boomerang),
and Classical, which are built according to formulae (15). The loops Crescent and Classical are hybrid loops since they combine a rectilinear sections and curvilinear ones. If bx=0, the loops (15) degenerate into rectangular loop (Non-ideal Relay without Whiskers).
Tilting the hybrid Classical loop in the split point a by means of the
phase shift ∆α1 is not always applicable since a jump of derivative on the
curvilinear section of the loop appears with such tilting. Loop tilting by
Fig. 8. Piecewise-linear
hysteresis
loops of Leaf (two-linear loop, Play
without Whiskers, n=1), hybrid Crescent (hybrid Boomerang, n=2), and
hybrid Classical (n=3) types built on
trapezoidal pulses.

means of rotation of coordinate system by angle θ allows to keep smoothness of the curvilinear section of the loop. However, an ordinary rotation
results in tilting of the horizontal straight sections. Tilting the loop by the
preset angle θ at the split point while keeping the straight sections horizontal allows the following transformation6

a[x (α ) y (T 4 ) − x (T 4 ) y (α )] + bx x (0 ) y (α )
,
x (0 ) y (T 4 )
y (α ) = y (α ),

x (α ) =

(16)

where x, y are functions (15) in which the loop parameters a and bx are replaced with the corresponding predistorted
values a and bx

a x = a cos θ ,
bx = bx cos θ − b y sin θ .

(17)

Transformations (17) rotate the split and the saturation points around the origin of coordinates by angle -θ. Determination of a predistorted value of the parameter by is not necessary because when tilting the loop through transformations (16), no changing of coordinate y(α) occurs (to avoid tilting of the horizontal sections of the loop).
6

R. V. Lapshin

Fig. 9. Hybrid Classical hysteresis (a)
(b)
loops with the specified slope at the Fig. 10. Piecewise-linear hysteresis loops of Leaf type (two-linear loop, Play
split point. The loops are built on with Gain without Whiskers) built on triangular pulses.
trapezoidal pulses.
Fig. 9 shows the Hybrid Classical hysteresis loops8 tilted through transformations (16). The smooth hysteresis loops (2) can also be tilted by transformations (16).6

2. Using triangular pulses instead of sine and cosine functions
Besides the trapezoidal pulses trp, triangular tri and rectangular rect pulses, which are particular cases of the
trapezoidal pulses,1 can be used in formulae (15). For example, loop of Leaf type can be obtained on the base of the
triangular pulses tri according to the following equations (m=n=1, ∆α1=∆α2=∆α3=0)

x (α ) = a tri c α + bx tri s α ,
y (α ) = by tris α ,

(18)

where tris, tric are triangular pulses of unit amplitude [tric(α)=tris(α+T/4)]. Piecewise-linear loop Leaf built according to formulae (18) is presented in Fig. 10(a). This loop is a variety of the loop Play with Gain without Whiskers.
The angles β and γ shown in the figure are determined by the formulae

β = arctan
γ = arctan

by
bx − a
by
bx + a

,
(19)
.

Using formulae (19), equations (18) can be expressed, if necessary, through the angles β and γ. With bx=a, loop
(18) is degenerating into the loop variety Non-ideal Relay with Gain without Whiskers.
Another method of description of Leaf type loop using triangular pulses is based on the following system of
equations6

x (α ) = bx tris α ,
 a

y (α ) = by 
tric α + tris α .
 a − bx


(20)

Fig. 10(b) shows piecewise-linear loop Leaf built by formulae (20) with bx>2a. This loop is also a variety of the
loop Play with Gain without Whiskers. The angles β and γ are determined according to the formulae6

β = arctan
γ = arctan

by
bx − a

,

by (bx − 2a )
bx (bx − a )

(b − 2a ) tan β .
= arctan x

(21)

bx

7

An improved parametric model of hysteresis loop

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

Fig. 11. Piecewise-linear hysteresis loops: (a) Play with Gain, (b) Play with Gain without Whiskers, (c) Play with
Attenuation, (d) Play with Attenuation without Whiskers, (e) Play (backlash), (f) Play without Whiskers, (g) Nonideal Relay with Gain, (h) Non-ideal Relay with Gain without Whiskers, (i) Non-ideal Relay with Attenuation, (j)
Non-ideal Relay with Attenuation without Whiskers, (k) Non-ideal Relay (Schmitt trigger), (l) Non-ideal Relay
without Whiskers (rectangular loop) built on triangular pulses. Loops (b)-(l) are special cases of the loop (a) Play
with Gain.
Using formulae (21), equations (20) can be expressed, if necessary, through the angles β and γ. With bx=2a, loop
(20) is degenerating into the loop Play without Whiskers. With bx<2a, loop (20) is degenerating into the loop Play
with Attenuation without Whiskers.
In the present study, general equations describing a piecewise-linear hysteresis loop Play with Gain [see
Fig. 11(a)] are obtained6
x (α ) =b x tri s α ,
  a tan β tric α +b x tan γ −b y
 1
y (α ) = 2(b y −b x tan γ ) H r 
+ tri s α  −  + b x tan γ tri s α ,
 2
  (a −b x ) tan β + b x tan γ

(22)

where Hr is a real (non-ideal) unit step function defined as follows

0, t < 0,

t
H r (t ) =  , 0 ≤ t ≤ t f ,
t f
1, t > t ,
f


8

(23)

R. V. Lapshin
where tf is a “front duration” of the real step function. The parameter tf is determined by loop parameters a, bx, by, β,

γ according to the following formula6
tf =

2(bx tan γ − by )

(a − bx ) tan β + bx tan γ

,

(24)

and, vice versa, any of the listed parameters can be determined by this formula by fixing front duration tf of the step
function. With tf→0 (β→90°), the function Hr degenerates into the ideal unit step function H (Heaviside function).
In that case, loops of Non-ideal Relay are formed [see Figs. 11(g)-(l)].
The loops shown in Figs. 11(b)-(l) are the particular cases of loop (22). Equations (22) are universal since they
allow to obtain the whole set of piecewise-linear loops of Play and Relay. Step-by-step derivation of formulae (22)
is given in the supplementary materials.6
Moving the splitting parameter a in formulae (20) into the argument of the triangular pulse function, other
equations describing the piecewise-linear hysteresis loop Play with Gain can be obtained6

x (α ) =b x tri s α ,
 

T
a tan β  
T
a tan β  
 +b tan γ tri s α .
y (α ) = (b y −b x tan γ ) tri s  α +  1 −
− tric  α + 1 −
 b tan β −b   x

8  b x tan β −b y  
8
 
x
y

 


(25)

The angle β is determined according to the formula6

β = arctan

2by − bx tan γ
bx

.

(26)

Equations (25) are less universal than equations (22). Nevertheless, loops similar to ones shown in Figs. 11(a)-(f)
can be built with the help of equations (25).

C. Double hysteresis loops
1. Linking loops at the saturation point b
To make the description of double hysteresis loops more accurate, the method suggested in Ref. 1 has been
improved. In particular, according to the improved method, the equations for a double loop non-self-crossing in the
origin of coordinates (0-shaped loop) are as follows (α=0…2π)

π


x (α ) = x  2α − sgn (π − α ) − ∆α 3  + bx sgn(π − α ),
2


π


y (α ) = y  2α − sgn (π − α ) − ∆α 3  + by sgn (π − α ),
2



(27)

where sgn(π-α)=(π-α)/|π-α| is the signum function. Double loop (27) is formed by linking two loops in the saturation point b where the generating function y(α) reaches its maximum value by. Fig. 12(a) shows an example of double non-self-crossing loop built by formulae (27).
According to (27), movement along the double loop starts at the point (0, 0), goes on counter-clockwise first
by the top loop then by the bottom one and finishes at the point (0, 0). The double loop (27) differs from the one
used before in that the starting point of the second loop coincides with the end point of the first one, and the end
point of the second loop with the starting point of the first one. Thus, the double loop (27) can be drawn continuously, “without lifting the pencil from the paper”, as the parameter α changes from 0 to 2π.
The loop with a self-crossing in the origin of coordinates (8-shaped loop) can be built according to the formu-

9

An improved parametric model of hysteresis loop

(a)

(b)

(c)

Fig. 12. Double (a), (c) smooth, (b) piecewise-linear hysteresis loops formed by linking two loops (a), (c) Classical,
(b) Play with Gain in the (a), (b) saturation point b, (c) point xmax. In the linking point, the loops can be both nonself-crossing and self-crossing.
lae (α=0…2π)


π
x (α ) = x  2α −  sgn(π − α ) − ∆α 3  + bx sgn(π − α ),
2



π


y (α ) = y  2α −  sgn(π − α ) − ∆α 3  + b y sgn(π − α ).
2



(28)

The loops (27) and (28) have no appearance differences.
Double non-self-crossing and self-crossing piecewise-linear loops are built by formulae (27) and (28), respectively, replacing π with T/2. An example of double non-self-crossing loop Play with Gain built by formulae (27)
with ∆α3=0 is shown in Fig. 12(b). Self-crossing loop has the same appearance. Double loops can be built similarly
from all other loops [see Figs. 11(b)-(l)] derived from the loop Play with Gain (22).

2. Linking loops at the point xmax
Besides linking two loops in the saturation point b, the loops can be linked at the point, where generating function x(α) reaches its maximum xmax. In order to determine the value of parameter αmax for which x(αmax)=xmax, the
equation dx (α ) dα = 0 |α =α max should be numerically solved. In the expanded form it can be written as follows
− aˆm sin(α max + ∆α1 )cos m −1 (α max + ∆α1 ) + bˆx n cos(α max + ∆α 2 )sin n −1 (α max + ∆α 2 ) = 0.

(29)

Double non-self-crossing loops linking in the xmax point are built according to the formulae

π
π


x (α ) = x 2α − sgn(π − α ) − + α max  + x (α max ) sgn(π − α ),
2
2


π
π


y (α ) = y  2α − sgn(π − α ) − + α max  + y (α max ) sgn(π − α ).
2
2



(30)

Double self-crossing loops linking in the xmax point are built according to the formulae


π
π
x (α ) = x  2α −  sgn(π − α ) − + α max  + x (α max )sgn(π − α ),
2
2



π
π


y (α ) = y  2α −  sgn(π − α ) − + α max  + y (α max )sgn(π − α ).
2
2



(31)

A double loop built by formulae (30) is shown in Fig. 12(c). A double loop built by formulae (31) has the same
appearance.

10

R. V. Lapshin

(a)

(b)

Fig. 14. Triple self-crossing hysteresis
Fig. 13. Double non-self-crossing hysteresis loops of Propeller type formed loop of the Classical type formed as a
by pinching a loop at the point of coordinate origin by means of (a) exponen- result of a foldover.
tiation, (b) zero splitting a and a phase shift ∆α2.

3. Pinching a loop in the origin

of coordinates
Double non-self-crossing loops of Propeller type can be formed by pinching a loop in the origin of coordinates
[see Fig. 13(a)]. The pinching is achieved by raising the generating functions x(α), y(α) to odd powers k and l according to expressions (13). The parameter αd can be found from equation y (α d ) = b y 2 written for a half-height
of the saturation point by

α d = arcsin

1
l

2 Bbyl −1

.

(32)

Solving the equation
x (α d ) − x (π − α d ) = 2a d ,

(33)

a split value a can be determined such that the split value ad of the double loop at the half-height of the saturation
point by is equal to the specified value. The analytical solution of equation (33) for k=3 is as follows
3

a=

(

A4 4ad + 4 4 A2bx6 sin 6 n α d + ad2

(

) − 4 A b sin
2

2 2
x

2n

2 A cosm α d 3 A2 4ad + 4 4 A2bx6 sin 6 n α d + ad2

)

αd

.

(34)

There is another method of building double non-self-crossing loops by means of loop pinching. This method
implies setting the split parameter a to zero and the phase shift ∆α2 (or ∆α3) to a non-zero value. With such loops,
the shift ∆α2 (or ∆α3) serves as a split parameter a. Fig. 13(b) shows an example of a double non-self-crossing
loop9, 10 obtained by pinching. The parameter αd can be found from equation y(αd)=by/2 written for a half-height of
the saturation point by. Solving numerically the equation
x (α d ) − x (π − α d ) = 2a d ,

(35)

it is possible to determine such a value of the phase shift ∆α2 (or ∆α3) that the split value ad of the double loop at
the half-height of the saturation point by is equal to the specified value.

D. Triple hysteresis loops
According to the improved model (2), the triple self-crossing loops10, 11 are obtained by setting up a phase shift
-∆α2 (or +∆α3) that “squeezes” the loop so tight that a foldover appears. Fig. 14 shows an example of a triple loop
of the Classical type formed as a result of a foldover.
11

An improved parametric model of hysteresis loop
In addition, a triple loop can be “assembled” of three loops – one central loop and two outside. Triple loops
linked at the saturation points b are built according to the formulae (α=0…2π)

 

π
3π 

x (α ) = (rect α + rect (α − π ))x1  3α −  + (1 − rect α − rect (α − π )) x2  ± 3α −
 + (b1 x + b2 x )sgn(π − α ),
2
2 

 

3π
 
π

y (α ) = (rect α + rect (α − π )) y1  3α −  + (1 − rect α − rect (α − π )) y2  ± 3α −
2
2

 

(36)



 + (b1 y + b2 y )sgn (π − α ),



where x1(α), y1(α) are equations (2) of the central loop; x2(α), y2(α) are equations (2) of the outside loops; b1x, b1y
are coordinates of the saturation point of the central loop; b2x, b2y are coordinates of the saturation points of the
outside loops; rect α=H(α)-H(α-π/3) is a π/3-wide rectangular pulse.
When assembling loops (36), the condition γ1=γ2 is valid, where γ1, γ2 are slope angles of tangents to the
unsplit (a1=0) central loop and the unsplit (a2=0) outside loops, respectively, at the saturation point b1. The slope
angles γ1, γ2 of the tangents are defined by the formulae6

γ 1 = arctan
γ 2 = arctan

b1 y
n1b1 x

,

b2 y
n2b2 x

(37)
.

Non-self-crossing loop will appear when the argument of functions x2(α), y2(α) is used with the plus sign in
equations (36) and self-crossing
loop – when it is used with the
minus sign. Non-self-crossing and
self-crossing loops have the same
appearance. The types of the central and the outside loops can be
diverse. Fig. 15(a) shows a triple
non-self-crossing

smooth

loop

(a=0.2, bx=2.0, by=2n2b2xtanγ1+b1y)
(a)

(b)

Fig. 15. Triple (a) smooth, (b) piecewise-linear hysteresis loops formed by
linking in the saturation point b of three loops (a) Classical, (b) Play with
Gain. In the linking points, the loops can be both non-self-crossing and selfcrossing.

built by equations (36). The loop
consists of central Classical loop
(1)

(a1=a,

b1x=1.0,

b1y=1.0,

m1=n1=3) and a pair of outside

Classical loops (1) (a2=0.1, b2x=(bx-b1x)/2, b2y=n2b2xtanγ1, m2=n2=3). The formulae for triple loops linked in the xmax
points are presented in the supplementary materials.6
Triple non-self-crossing and self-crossing piecewise-linear loops are built by formulae (36) by replacing π by
T/2. Fig. 15(b) shows a triple non-self-crossing piecewise-linear loop (a=0.4, bx=2.0, by=2.0) built by equations
(36). The loop consists of central Play with Gain loop (22) (a1=a, b1x=1.0, b1y=1.0, β 1=77°, γ1=17°) and a pair of
outside Play with Gain loops (22) (a2=0.2, b2x=(bx-b1x)/2, b2y=(by-b1y)/2, β 2=77°, γ2=γ1). Triple loops can be built
similarly from all other piecewise-linear loops [see Figs. 11(b)-(l)], which are the particular cases of the loop Play
with Gain.
Triple loops (36) are convenient in the cases when single smooth or hybrid hysteresis loops having long
12

R. V. Lapshin
“whiskers”12, 13 are required. Recall
that long whiskers in the model (1)
can be obtained by increasing
power m.6 However, at the same
time, a change of loop curvature
occurs that is not always necessary.
Fig. 16(a) shows the simulation of
a single smooth loop Classical with
long

whiskers

(a=0.2,

bx=2.0,

by=(bx-b1x)tanγ1+b1y) by using the
triple

non-self-crossing

smooth

hysteresis loop. The triple loop

(a)

(b)

Fig. 16. Simulation of a single loop with long whiskers (a) Classical, (b)
Hybrid Classical by means of the triple non-self-crossing hysteresis loop.
Whiskers are the outside pair of (a) smooth, (b) piecewise-linear unsplit
loops of Leaf type.

consists of the central Classical
loop (1) (a1=a, b1x=1.0, b1y=1.0, m1=n1=3) and whiskers which are formed from a pair of the outside unsplit loops
of Leaf type (1) (a2=0, b2x=(bx-b1x)/2, b2y=(by-b1y)/2, m2=3, n2=1) oriented at the angle γ1.
Fig. 16(b) shows an example of building hybrid Classical loop with whiskers6 (a=0.4, bx=1.0, by=1.0) of three
loops. The loop consists of the central hybrid Classical loop (16) (a1=a, b1x=0.3, b1y=by, n1=5, θ1=5°) and whiskers
which are formed from a pair of the outside unsplit piecewise-linear loops Play with Gain without Whiskers (18)
(a2=0, b2x=(bx-b1x-a)/2, b2y=0). Piecewise-linear loops with whiskers of any length are built directly by equations
(22).

E. Area of a hysteresis loop
The area of a hysteresis loop characterizes energy losses in a piezoelectric/ferromagnetic material while applying an alternating electric/magnetic field. To find the hysteresis loop area, the following general formula is used
S=

1 
dy (α )
dx (α ) 
− y (α )
 x (α )
dα .
∫
2 
dα
dα 

(38)

By substituting the expressions x(α), y(α) of the improved model (2) along with their derivatives into (38) and integrating, one can obtain
m −1
m +1
n −1
n +1
 aˆ cos(∆α − ∆α ) 
m − 1 2  bˆx sin(∆α 2 − ∆α 3 ) 
n − 1 2  
1
3 
2 −
2 −
+
S=
2
mC
C
2
nC
C
m −1
m +1 
n −1
n +1  πb y ,


2
2
2m
2n




 

(39)

where C lk = l! [k!(l − k )!] are binomial coefficients. Formula (39) is valid for any positive integer m and n.
The loop tilt produced by a phase shift ∆α1 (m>1, ∆α2=∆α3=0) results in an increase in the loop area S (see
Fig. 2). The area S of a hysteresis loop increases with an increase in the phase shift ∆α2 (see Fig. 3). With an increase in the phase shift ∆α3, the area S decreases (see Fig. 4).
According to (38), as a hysteresis loop is being scaled as x (α ) = Ax (α ) , y (α ) = By (α ) , its area S changes
proportionally to the product of the scale factors A and B, i. e., S = ABS . A series of three hysteresis loops, where
the areas are doubling with scaling, is shown in Fig. 17(a).
It follows from formula (39) that with m=1 and ∆α2=∆α3, loop area

13

An improved parametric model of hysteresis loop
S = πaby .

(40)

Thus, under the specified conditions the loops will have the same area regardless of their type, saturation bx and
phase shifts ∆α provided that the parameters a and by are the same. Fig. 17(b) shows an example of loops of the
same area having different phase shifts ∆α1. The formula of areas of the Leaf loops with m=1 and ∆α1=∆α3 looks
the same as (40); the areas of these loops are also independent of the values of saturation bx and phase shifts ∆α.
From formula (39) it also follows that with m>1 and ∆α2=∆α3, loop area
m −1

m +1

m −1 2
C m +1
2
S= m
πab y .
2 cos m −1 (∆α 1 − ∆α 3 )
2mC m −21 −

(41)
Thus, under the specified conditions the loops will have the same
area S regardless of their type and
saturation bx provided that the parameters a, by and m are the same
[see Fig. 17(c)].
It follows from formula (41)
that in case of ∆α2=∆α3=0, the area
of a loop having a negative slope

(a)

(b)

(∆α1>0) at the split point equals to
the area of a loop having a positive
slope (∆α1<0) at the split point,
i. e., S + ∆α = S − ∆α . The inequal1

ity

S

∆α1 − ∆α 3 > 0

1

> S ∆α − ∆α
1

3 =0

also

follows from formula (41). Thus,
for any loop type in case of
∆α2=∆α3=0 (m>1), the energy
losses in the tilted loops (|∆α1|>0,

(c)
(d)
Fig. 17. The hysteresis loop area (a) is proportional to the product of the
scale factors A and B along axes x and y, respectively, (b) is invariant to the
phase shifts ∆α1 (restrictive conditions are: m=1, ∆α2=∆α3), (c), (d) is invariant to the saturations bx (restrictive condition is: ∆α2=∆α3).

see Fig. 2) are greater than in the
“upright” ones (∆α1=0).
By substituting the zero phase
shifts ∆α1=∆α2=∆α3=0 into (39),

the following formula is obtained
m −1
m +1

m − 1 2  1
S =  2mCm −21 −
Cm +1 m πaby

2
2



(42)

by which the area of the loop built according to the model (1) can be calculated. Formula (42) provides a solution
for any positive integer m.
There is an inaccuracy in Ref. 1 in formula (27) and in the accompanying text to this formula. Formula (27) is
valid for any n (even or odd). Since formula (27) does not have a solution with m=1, formula (42) presented in this
14

R. V. Lapshin
Table. Average relative approximation error δ
teresis loop.

(%) of the original and the improved parametric models of hys-

Loop Leaf,

Loop Classical,

Loop Classical,

Fig. 18(a)

Fig. 18(b)

Fig. 18(c)

Original

0.8

2.9

1.7

Improved

0.5

1.0

1.0

Model

paper should be used instead of it. By the same reason instead of formula (28) from Ref. 1, the formula presented in
the supplementary materials6 should be used.
Since n does not enter into formula (42), all three loop types Leaf, Crescent, and Classical of the model (1)
have the same loop area S provided that the other parameters (a, by, and m) of these loops are the same (see Fig. 1).
Since parameter bx does not enter into formula (42) as well, all the loops (1) having the same parameters a, by, and
m but different bx have the same area S [see Fig. 17(d)].
It is noteworthy that with ∆α1=∆α2=∆α3=∆α≠0, where ∆α is an arbitrary real number, formula (39) also takes
the form (42). In this case, loops built according to models (1) and (2) will have equal areas. This is because under
the above conditions, these loops are of exactly the same shape (see Section II.A.4).
The area of piecewise-linear hysteresis loops is determined by the general formula (38).6 The easiest way to
determine the area of the piecewise-linear hysteresis loops is using the formula of area of a parallelogram. In that
case, the working formula for the area of the loop Play with Gain (22) [see Fig. 11(a)] is as follows6

[

]

 (a − bx ) tan γ + by tan β a sin β sin γ
S = 4a 
−
tan β − tan γ
sin(β − γ )



.


(43)

The area of any piecewise-linear loop [see Figs. 11(b)-(l)] derived from the loop Play with Gain (22) can be determined by formula (43).

III. APPLICATION OF THE IMPROVED MODEL
The developed model enables building smooth, piecewise-linear, hybrid, minor, mirror-reflected, inverse, double and triple loops. Using the phase shifts ∆α1, ∆α2, ∆α3 as additional parameters allows reducing the approximation error of hysteresis loops by several times.
Error analysis conducted according to the method described in Ref. 1 has shown that average relative approximation error of the improved model does not exceed 1% (see Table).6 For the comparison to be correct, the hysteresis loops from Ref. 1 were used during the error determination. Fig. 18 shows approximating loops built according to the existing (1) and the improved (2) models overlaid with the experimental hysteresis loops taken from
Ref. 1. It is easy to see that the approximation accuracy of the improved model is notably higher than the approximation accuracy of the existing model.
Like the original model, the improved one can be applied for correction of the distortions caused by hysteresis
of piezomanipulators [see Fig. 18(a)] of the scanning probe microscope (SPM). To do this, phase-shifting elements
providing the phase shifts ∆α1, ∆α2 (or ∆α1, ∆α3) should be added to the existing hardware suggested in Ref. 1.
Calculation of derivatives, searching for the harmonically linearized transfer function of a hysteresis element and
the inverse function in the improved model can be done similarly to Ref. 1.

15

An improved parametric model of hysteresis loop

(a)

(b)

(c)

Fig. 18. Approximation of real hysteresis loops of types (a) Leaf, (b), (c) Classical. The red uneven loop is experimental; the blue dotted loop is the existing model; the green smooth one is the suggested improved model. Approximation error of the improved model does not exceed 1%.

IV. SUMMARY
Improvements have been made to the existing model of a hysteresis loop built on parametric equations. As a
result, the approximation accuracy has increased several times. The possibility is demonstrated that various types of
hysteresis loops can be manipulated by means of summation/subtraction operations. A general formula has been
derived for building piecewise-linear hysteresis loops of Play and Non-ideal Relay types as well as their variations
widely used in simplified models of hysteresis phenomena. Several new formulae describing various types of double and triple loops have been obtained. A more general expression has been suggested to determine the hysteresis
loop area, which allows for evaluation of energy losses in piezoelectric/ferromagnetic materials. The hysteresis
loop model developed is especially suitable for solving the tasks of simulation of cyclically operating instruments
that include hysteresis elements.14

ACKNOWLEDGMENTS
Ministry of Education and Science of the Russian Federation supported this work. I thank Oleg E. Lyapin,
Larisa B. Sharova, and David W. Waddell for critical reading of the manuscript and checking of the supplementary
materials; Assoc. Prof. Eugene A. Fetisov for support and stimulation.

References
1

R. V. Lapshin, Analytical model for the approximation of hysteresis loop and its application to the scanning tun-

neling microscope, Rev. Sci. Instrum. 66, 4718 (1995) (available at www.niifp.ru/staff/lapshin/en/#articles).
2

B. Mokaberi and A. A. G. Requicha, Compensation of scanner creep and hysteresis for AFM nanomanipulation,

IEEE Trans. Autom. Sci. Eng. 5, 197 (2008).
3

B. Graffel, F. Müller, A.-D. Müller, and M. Hietschold, Feedforward correction of nonlinearities in piezoelectric

scanner constructions and its experimental verification, Rev. Sci. Instrum. 78, 053706 (2007).
4

V. Hassani, T. Tjahjowidodo, and T. N. Do, A survey on hysteresis modeling, identification and control, Mech.

Syst. Sig. Process. 49, 209 (2014).
5

R. V. Lapshin, Improved approximating model of hysteresis loop for the linearization of a probe microscope pie-

zoscanner, Proc. of the XIX Russian Symposium on Scanning Electron Microscopy and Analytical Methods of

16

R. V. Lapshin
Investigation of Solids, pp. 154-155, Chernogolovka, Russian Federation, June 1-4, 2015 (in Russian, available at
www.niifp.ru/staff/lapshin/en/#reports).
6

Supplementary materials: R. V. Lapshin, Hysteresis loop, Mathcad worksheet, version-date 03.15.2017;

R. V. Lapshin, Hysteresis loop, Readable Mathcad worksheet, 640 pp., 2017 (last versions are available at
www.niifp.ru/staff/lapshin/en/#downloads).
7

Y. Shimizu, K. Matsuda, M. Mizutani, K. Nishimura, T. Kawabata, S. Ikeno, Y. Hishinuma, and S. Aoyama,

Superconducting properties of MgB2 particle impregnated with Mg-based alloys, Mater. Trans. 52, 272 (2011).
8

S. H. Baek and C. B. Eom, Reliable polarization switching of BiFeO3, Phil. Trans. R. Soc. A 370, 4872 (2012).

9

Y. V. Pershin and M. Di Ventra, Memory effects in complex materials and nanoscale systems, Adv. Phys. 60, 145

(2011).
10

S. Sawada, T. Yamaguchi, H. Suzuki, and F. Shimizu, Experimental studies on phase transitions in ferroelectric

{N(CH3)4}2ZnCl4, J. Phys. Soc. Jpn. 54, 3129 (1985).
11

K. Dumesnil, C. Dufour, Ph. Mangin, and A. Rogalev, Location of magnetic springs in exchange-coupled super-

lattices, ESRF Highlights 2002, pp. 43-45 (2002).
12

H. Khanduri, M. C. Dimri, H. Kooskora, I. Heinmaa, G. Viola, H. Ning, M. J. Reece, J. Krustok, and R. Stern,

Structural, dielectric, magnetic, and nuclear magnetic resonance studies of multiferroic Y-type hexaferrites,
J. Appl. Phys. 112, 073903 (2012).
13

F. J. Muñoz Lara, A. B. Gaspar, D. Aravena, E. Ruiz, M. Carmen Muñoz, M. Ohba, R. Ohtani, S. Kitagawa, and

J. A. Real, Enhanced bistability by guest inclusion in Fe(II) spin crossover porous coordination polymers, Chem.
Commun. 48, 4686 (2012).
14

S. A. Agafonov and V. A. Matveev, Dynamics of a balanced rotor under the action of an elastic force with a

hysteresis characteristic, Mech. Solids 47, 160 (2012).

17

