How to avoid the overshoot in Jiles Atherton hysteresis loop

Trying to plot a hysteresis model using jiles atherton hysteresis model in matlab. All the parameters are defined in the beginning.

Here is the code I have written,

clear

%% Given parameters

a = 512;

c = 0.75;

k = 570;

M_s = 1.7*10^6;

i_s = 10^-6;

alpha = 0.032;

u_0 = 4*pi*10^-7;

f = 1;

h_max = 8000;

%% Generate time vector

t_start = 0; %start time

t_end = 60; %end time

dt = 0.001; %time steps

t = t_start:dt:t_end;

%% Initial values

M_irev(1) = 0; %Initial irreversible magnetization

h = h_max*sin(2*pi*f*t); % Sinusoidal magnetic field strength

%% Time stepping loop

for i=2:length(t)

% Compute effective field

h_e = h(i); % Applied field and effective field considered equal

% Compute anhysteric magnetization

M_an = M_s*(coth(h_e/a) - (a/(h_e))); % Langevin function

% Compute derivative of magnetic field with respect to time

dh_e = (h(i) - h(i-1)) / dt;

delta = sign(dh_e);

% Compute irreversible magnetization using Euler method

M_irev(i) = ((M_an - M_irev(i-1)) / ((k * delta) - (alpha * (M_an - M_irev(i-1))))) * dh_e;

% Update magnetization

M(i) = c*M_an + (1-c)*M_irev(i);

end

%% Magnetic flux density

B = u_0 * (h + M);

%% Plot results

figure;

plot(h, B);

xlabel('Magnetic Field Strength');

ylabel('Magnetic Flux Density');

title('HB Graph');

grid on;

The problem I am facing is when it plots the curve, it shows overshoot at three points. Suspecting the problem is in the 'Langevin function'. Lower values of h at the end part creates bigger values thus it shows overshoot at those points. I think I need to use approximation (Taylor series or similar method) at those points so that when overshoot happens the code will avoid those overshoots and take the approximation. Can anyone show me how can I do that approximation?

[Josefina]

Read the recent work: "Blind Efficient Method for Optimizing Jiles-Atherton Model Parameters", TMAG, (2025). The study provides a detailed decription the the implementation of the JA model in MATLAB, together with a novel fast and blind method for retrieving the 5 standard model parameters for a given material.

You may contact the corresponding author for any help or collaboration.

Answer 1

[Mathieu NOE]

hello

I am just trying to remove the spikes (with filloutiers) and further smooth the curve with a touch of smoothdata

look at these 2 lines I added at the plot section

Bs= filloutliers(B,'linear','movmedian',15, 'ThresholdFactor', 1.1);

Bs= smoothdata(Bs,'sgolay',50);

code :

clear
%% Given parameters
a = 512;
c = 0.75;
k = 570;
M_s = 1.7*10^6;
i_s = 10^-6;
alpha = 0.032;
u_0 = 4*pi*10^-7;
f = 1;
h_max = 8000;
%% Generate time vector
t_start = 0; %start time
t_end = 60; %end time
dt = 0.001; %time steps
t = t_start:dt:t_end;
%% Initial values
M_irev(1) = 0; %Initial irreversible magnetization
h = h_max*sin(2*pi*f*t); % Sinusoidal magnetic field strength
%% Time stepping loop
for i=2:length(t)
    % Compute effective field
    
    h_e = h(i); % Applied field and effective field considered equal
    % Compute anhysteric magnetization
    M_an = M_s*(coth(h_e/a) - (a/(h_e))); % Langevin function
    % Compute derivative of magnetic field with respect to time
    dh_e = (h(i) - h(i-1)) / dt;
    delta = sign(dh_e);
    % Compute irreversible magnetization using Euler method
    M_irev(i) = ((M_an - M_irev(i-1)) / ((k * delta) - (alpha * (M_an - M_irev(i-1))))) * dh_e;
    
    % Update magnetization
    M(i) = c*M_an + (1-c)*M_irev(i);
end
%% Magnetic flux density
B = u_0 * (h + M);
%% Plot results
figure;
Bs= filloutliers(B,'linear','movmedian',15, 'ThresholdFactor', 1.1);
Bs= smoothdata(Bs,'sgolay',50);
plot(h, B, h, Bs);
xlabel('Magnetic Field Strength');
ylabel('Magnetic Flux Density');
title('HB Graph');
grid on;

[Mathieu NOE]

we could also do a cycle by cycle averaging and get the one cycle averaged curve , but at the end the result may look quite the same (almost , stil there is a small discontinuity, we can put a bit of smoothing there)

for the time being, this is the result without smoothing :

%% Given parameters
a = 512;
c = 0.75;
k = 570;
M_s = 1.7*10^6;
i_s = 10^-6;
alpha = 0.032;
u_0 = 4*pi*10^-7;
f = 1;
h_max = 8000;
%% Generate time vector
t_start = 0; %start time
t_end = 60; %end time
dt = 0.001; %time steps
t = t_start:dt:t_end;
%% Initial values
M_irev(1) = 0; %Initial irreversible magnetization
h = h_max*sin(2*pi*f*t); % Sinusoidal magnetic field strength
%% Time stepping loop
for i=2:length(t)
    % Compute effective field
    
    h_e = h(i); % Applied field and effective field considered equal
    % Compute anhysteric magnetization
    M_an = M_s*(coth(h_e/a) - (a/(h_e))); % Langevin function
    % Compute derivative of magnetic field with respect to time
    dh_e = (h(i) - h(i-1)) / dt;
    delta = sign(dh_e);
    % Compute irreversible magnetization using Euler method
    M_irev(i) = ((M_an - M_irev(i-1)) / ((k * delta) - (alpha * (M_an - M_irev(i-1))))) * dh_e;
    
    % Update magnetization
    M(i) = c*M_an + (1-c)*M_irev(i);
end
%% Magnetic flux density
B = u_0 * (h + M);

%% cycle by cycle averaging 
dtzc = 1/f;
tzc = t(1):dtzc:t(end);
NN = numel(tzc)-1;
hh = 0;
BB = 0;

for k = 1:NN 
    tstart = tzc(k);
    tstop = tzc(k+1);
    
    tt = linspace(tstart,tstop,360); % 360 => 1 degree resolution
    hi = interp1(t,h,tt);
    Bi = interp1(t,B,tt);
    hh = hh + hi;
    BB = BB + Bi;
    
end

% divide by N to obtain the averaged data
    hh = hh/NN;
    BB = BB/NN;

%% Plot results
figure;
plot(h, B, hh, BB);
xlabel('Magnetic Field Strength');
ylabel('Magnetic Flux Density');
title('HB Graph');
grid on;

[Mathieu NOE]

just added now a bit of smoothing on top

see the effect here : smoothed curve is yellow

%% Given parameters
a = 512;
c = 0.75;
k = 570;
M_s = 1.7*10^6;
i_s = 10^-6;
alpha = 0.032;
u_0 = 4*pi*10^-7;
f = 1;
h_max = 8000;
%% Generate time vector
t_start = 0; %start time
t_end = 60; %end time
dt = 0.001; %time steps
t = t_start:dt:t_end;
%% Initial values
M_irev(1) = 0; %Initial irreversible magnetization
h = h_max*sin(2*pi*f*t); % Sinusoidal magnetic field strength
%% Time stepping loop
for i=2:length(t)
    % Compute effective field
    
    h_e = h(i); % Applied field and effective field considered equal
    % Compute anhysteric magnetization
    M_an = M_s*(coth(h_e/a) - (a/(h_e))); % Langevin function
    % Compute derivative of magnetic field with respect to time
    dh_e = (h(i) - h(i-1)) / dt;
    delta = sign(dh_e);
    % Compute irreversible magnetization using Euler method
    M_irev(i) = ((M_an - M_irev(i-1)) / ((k * delta) - (alpha * (M_an - M_irev(i-1))))) * dh_e;
    
    % Update magnetization
    M(i) = c*M_an + (1-c)*M_irev(i);
end
%% Magnetic flux density
B = u_0 * (h + M);

%% cycle by cycle averaging 
dtzc = 1/f;
tzc = t(1):dtzc:t(end);
NN = numel(tzc)-1;
hh = 0;
BB = 0;

for k = 1:NN 
    tstart = tzc(k);
    tstop = tzc(k+1);
    
    tt = linspace(tstart,tstop,360); % 360 => 1 degree resolution
    hi = interp1(t,h,tt);
    Bi = interp1(t,B,tt);
    hh = hh + hi;
    BB = BB + Bi;
    
end

% divide by N to obtain the averaged data
    hh = hh/NN;
    BB = BB/NN;

Bs= smoothdata(BB,'loess',20);
%% Plot results
figure;
plot(h, B, hh, BB, hh, Bs);
xlabel('Magnetic Field Strength');
ylabel('Magnetic Flux Density');
title('HB Graph');
grid on;

[Mathieu NOE]

hello again

problem solved ?

[Md Tanbir Siddik]

Hi,

Thank you for your effort in finding solution. actually it is not about smoothing the signal. The main goal is to plot the figure showed below,

the problem in the code is to add estimation when there is overshoot. The overshoot happens due to dividing by very less value of 'h_e' at the end of Langevin function. To tackle this problem, I need to use an approximation (e.g. Taylor syries approximation) when it starts overshooting. Maybe I should use "if" condition somewhere and implement approximation method when overshoot starts but I exactly can't figure out where and how to implement.

[Mathieu NOE]

I am not sure it's the Langevin function which causes the issue

I traced down a few quantities (so I had to mdify a bit your code

if I plot Langevin function output (M_an), there is not much issue to be seen , except one single peak)

it seems to me there is a bigger issue in the Euler equation , if you trace down the numerator and denominator, it shows a lot of spikes , so there is some discontinuities in your process that must be corrected

%% Given parameters
a = 512;
c = 0.75;
k = 570;
M_s = 1.7*10^6;
i_s = 10^-6;
alpha = 0.032;
u_0 = 4*pi*10^-7;
f = 1;
h_max = 8000;
%% Generate time vector
t_start = 0; %start time
t_end = 60; %end time
dt = 0.001; %time steps
t = t_start:dt:t_end;
%% Initial values
M_irev(1) = 0; %Initial irreversible magnetization
h = h_max*sin(2*pi*f*t); % Sinusoidal magnetic field strength
%% Time stepping loop
for i=2:length(t)
    % Compute effective field
    
    h_e = h(i); % Applied field and effective field considered equal
    % Compute anhysteric magnetization
    M_an(i) = M_s*(coth(h_e/a) - (a/h_e)); % Langevin function
    % Compute derivative of magnetic field with respect to time
    dh_e(i) = (h(i) - h(i-1)) / dt;
    delta(i) = sign(dh_e(i));
    % Compute irreversible magnetization using Euler method
    num(i) = (M_an(i) - M_irev(i-1))* dh_e(i);
    den(i) = ((k * delta(i)) - (alpha * (M_an(i) - M_irev(i-1))));
%     M_irev(i) = ((M_an(i) - M_irev(i-1)) / ((k * delta) - (alpha * (M_an(i) - M_irev(i-1))))) * dh_e(i);
    M_irev(i) = num(i) / den(i);
    
    % Update magnetization
    M(i) = c*M_an(i) + (1-c)*M_irev(i);
end
%% Magnetic flux density
B = u_0 * (h + M);
%% Plot results
figure;
plot(h, B);
% plot(h, B, hh, BB, hh, Bs);
xlabel('Magnetic Field Strength');
ylabel('Magnetic Flux Density');
title('HB Graph');
grid on;

%% more plots (for debug purposes)
figure;
plot(M_an);
title('M-an');

figure;
plot(delta);
title('delta');

figure;
plot(num);
title('M-irev num');

figure;
plot(den);
title('M-irev den');

[Md Tanbir Siddik]

Hi,

Thank you for your continuous support. The code with interpolation and smoothing solves the problem for now. I somehow overlooked that response when I was replying. It seems when I further decrease the step size from 0.0001, the curve shows some wave at those overshooting points. which later creates larger wave in the magnetostriction curve as well. But 0.0001 step size is fine for me till now. Here is the code that I am using for now,

clear

%% Given parameters

a = 512;

c = 0.75;

k = 570;

M_s = 1.7*10^6;

l_s = 10*10^-6;

alpha = 0.032;

u_0 = 4*pi*10^-7;

f = 1;

h_max = 8000;

%% Generate time vector

t_start = 0; %start time

t_end = 60; %end time

dt = 0.0001; %time steps

t = t_start:dt:t_end;

%% Initial values

M_irev(1) = 0; %Initial irreversible magnetization

h = h_max*sin(2*pi*f*t); % Sinusoidal magnetic field strength

%% Time stepping loop

for i=2:length(t)

% Compute effective field

h_e = h(i); % Applied field and effective field considered equal

% Compute anhysteric magnetization

M_an = M_s*(coth(h_e/a) - (a/(h_e))); % Langevin function

% Compute derivative of magnetic field with respect to time

dh_e = (h(i) - h(i-1)) / dt;

delta = sign(dh_e);

% Compute irreversible magnetization using Euler method

M_irev(i) = ((M_an - M_irev(i-1)) / ((k * delta) - (alpha * (M_an - M_irev(i-1))))) * dh_e;

% Update magnetization

M(i) = c*M_an + (1-c)*M_irev(i);

end

%% Magnetic flux density

B = u_0 * (h + M);

%% cycle by cycle averaging

dtzc = 1/f;

tzc = t(1):dtzc:t(end);

NN = numel(tzc)-1;

hh = 0;

BB = 0;

MM = 0;

for k = 1:NN

tstart = tzc(k);

tstop = tzc(k+1);

tt = linspace(tstart,tstop,360); %360 => 1 degree resolution

hi = interp1(t,h,tt);

Bi = interp1(t,B,tt);

Mi = interp1(t,M,tt);

hh = hh + hi;

BB = BB + Bi;

MM = MM + Mi;

end

% divide by N to obtain the averaged data

hh = hh/NN;

BB = BB/NN;

MM = MM/NN;

Bs= smoothdata(BB,'loess',20);

hs= smoothdata(hh,'loess',20);

Ms= smoothdata(MM,'loess',20);

%% Compute magnetostriction

l_ms = l_s*((Ms/M_s).^2);

%% Plot results

%Magnetic Field Strength vs Magnetic Flux Density

subplot(3,1,1);

plot(hh, Bs, "k");hold on;

%plot(h, B, hh, BB, hh, Bs); %Validate the estimation

xlabel('Magnetic Field Strength (H)');

ylabel('Magnetic Flux Density (B)');

title('HB Graph');

grid on;

%Magnetostriction As Function of Magnetic Flux Density

subplot(3,1,2);

plot(Bs,l_ms, "g");

xlabel('Magnetic Flux Density (B)');

ylabel('Magnetostriction (l_ms)');

title('Magnetostriction As Function of Magnetic Flux Density');

grid on;

%Magnetostriction As Function of Magnetic Field Strength

subplot(3,1,3);

plot(hs,l_ms, "b");

xlabel('Magnetic Field Strength (H)');

ylabel('Magnetostriction (l_ms)');

title('Magnetostriction As Function of Magnetic Field Strength');

grid on;

[Mathieu NOE]

ok, glad you have now a work around !