# How to avoid the overshoot in Jiles Atherton hysteresis loop

11 views (last 30 days)
Md Tanbir Siddik on 25 Apr 2024
Commented: Mathieu NOE on 6 May 2024
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?

Mathieu NOE on 25 Apr 2024
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;
Md Tanbir Siddik on 6 May 2024
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 on 6 May 2024
ok, glad you have now a work around !

### Categories

Find more on Programming in Help Center and File Exchange

R2023b

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!