# How can I speed up the solving of this specific equation

I'm trying to work on solving the Jiles-Atherton model of magnetic hysteresis, so I'm using the following equations, as found in this paper (Jaafar, 2013)

The mathematica code intended to solve 1 time step is as follows

Subscript[M, s] = 1.7*10^6;
a = 1000;
c = 0.1;
k = 2000;
\[Alpha] = 0.001;
H[t_]:=(4*10^-5 + 0.1 Sin[2 \[Pi] t])
dH[t_] := \[Pi]/5 Cos[2 \[Pi] t]

NSolve[{
He[t] == H[t] + \[Alpha] M[t],
Man[t] == Subscript[M, s] (Coth[He[t]/a] - a/He[t]),
(Mirr[t] - Mirr[t - 1])/(H[t] - H[t - 1]) == (Man[t] - Mirr[t])/(
k Sign[dH[t]]),
(M[t] - M[t - 1])/(
H[t] - H[t - 1]) == (1 - c) (Man[t] - Mirr[t])/(
k Sign[dH[t]] - \[Alpha] (Man[t] - Mirr[t])) +
c (Man[t] - Man[t - 1])/(H[t] - H[t - 1])
}, {M[t], Man[t], Mirr[t], He[t]}] /. t -> 10^-3


The code above is iterated from 0 to 1.25 seconds using a do loop. However, NSolve takes very long to solve even 1 time step. How can I speed this up? Any help is appreciated, thank you.

Update: Thanks to Alex trounev for the suggestion. The code I have used is as follows (quite a Naive solution):

NDSolve[{
NDSolve[{
He[t] == H[t] + \[Alpha] M[t],
Man[t] == Subscript[M, s] (Coth[He[t]/a] - a/He[t]),
M'[t] == (Man[t] - Mirr[t])/(k Sign[dH[t]]) dH[t],
M'[t] == ((1 - c) (Man[t] - Mirr[t])/(
k Sign[dH[t]] - \[Alpha] (Man[t] - Mirr[t]))) dH[t] + c Man'[t],
He[0] == H[0],
M[0] == 0,
Man[0] == 0,
Mirr[0] == 0
}, {M, Man, Mirr, He}, {t, 0, 10}]


However, this now returns the error: NDSolve::icfail: Unable to find initial conditions that satisfy the residual function within specified tolerances. Try giving initial conditions for both values and derivatives of the functions. What does this mean and how can I fix this error?

• I don't think that this system can be simpified, I ran it through both Simplify and FullSimplify and there is no simplification. Commented Aug 15, 2023 at 6:22
• We can't solve this system of delayed equations using NSolve. It could be better to solve original differential equation (9) using NDSolve. Commented Aug 15, 2023 at 10:28
• Please, pay attention that H[t] - H[t - 1]=0. In this regards it could be better to use H[t]-H[t-dt]. Commented Aug 15, 2023 at 14:30
• Thanks for the help. NDSolve returns the error: NDSolve::idelay: Initial history needs to be specified for all variables for delay-differential equations. I'm not sure what this means Commented Aug 16, 2023 at 3:13
• The equation Man[t] == Subscript[M, s] (Coth[He[t]/a] - a/He[t]) does not agree with your initial conditions at t -> 0. If you substitute the initial conditions in this equation, it reduces to 0 == 0.025332. That's the problem. Commented Aug 16, 2023 at 8:27

There are 2 differential equations only therefore solution to this problem is given by

Subscript[M, s] = 1.7*10^6;
a = 1000;
c = 0.1;
k = 2000;
\[Alpha] = 0.001;
H[t_] := (4*10^-5 + 0.1 Sin[2 \[Pi] t]);
dH[t_] := \[Pi]/5 Cos[2 \[Pi] t]; He[t_] := H[t] + \[Alpha] M[t];
Man[t_] := Subscript[M, s] (Coth[He[t]/a] - a/He[t]);

sol = NDSolve[{Mirr'[t] == (Man[t] - Mirr[t])/(k Sign[dH[t]]) dH[t],
M'[t] == ((1 -
c) (Man[t] -
Mirr[t])/(k Sign[dH[t]] - \[Alpha] (Man[t] - Mirr[t]))) dH[
t] + c Man'[t], M[0] == 0, Mirr[0] == 0}, {M, Mirr}, {t, 0, 1},
Method -> "ExplicitEuler", StartingStepSize -> 0.001,
MaxSteps -> Infinity];


Visualization

fig1={ParametricPlot[{H[t], M[t] /. sol[[1]]}, {t, 0, 1},
AspectRatio -> 1/2, PlotRange -> All, Frame -> True,
FrameLabel -> {"H", "M"}],
ParametricPlot[{H[t], M[t] /. sol[[1]]}, {t, 0, 1},
AspectRatio -> 1/2, PlotRange -> {{-.01, .01}, {-.5, .5}},
Frame -> True, FrameLabel -> {"H", "M"}]}


Update 1. We can reproduce data above using FDM as follows

H[t_] := (4*10^-5 + 1/10 Sin[2 \[Pi] t]);
dH[t_] := \[Pi]/5 Cos[2 \[Pi] t];

sol0 = Solve[{(Mirr[t] -
Mirr[t - dt]) == (Man[t] - Mirr[t])/(k Sign[dH[t]]) (H[t] -
H[t - dt]), (M[t] -
M[t - dt]) == ((1 -
c) (Man[t] -
Mirr[t])/(k Sign[dH[t]] - \[Alpha] (Man[t] - Mirr[t]))) (H[
t] - H[t - dt]) + c (Man[t] - Man[t - dt])}, {M[t], Mirr[t]}];

M1 = M[t] /. sol0[[1]];

Mirr1 = Mirr[t] /. sol0[[1]];

Subscript[M, s] = 17*10^5; dt = 5 10^-4;
a = 1000.;
c = 1/10.;
k = 2000.;
\[Alpha] = .001; M[-dt] = 0.; Mirr[-dt] = 0.;
Man[-dt] = N[Subscript[M, s] (Coth[H[-dt]/a] - a/H[-dt])]; Do[
He[t] = H[t - dt] + \[Alpha] M[t - dt];
Man[t] = Subscript[M, s] (Coth[He[t]/a] - a/He[t]); Mp = M1;
Mirrp = Mirr1; He[t] = H[t] + \[Alpha] Mp; M[t] = M1;
Mirr[t] = Mirr1, {t, 0, 1, dt}]


Visualization together with data shown in Figure1:

Show[fig1[[2]],ListPlot[Table[{H[t], M[t]}, {t, 0, 1, dt}],
PlotRange -> {{-.01, .01}, {-.5, .5}}, FrameLabel -> {"H", "M"},
Frame -> True, PlotStyle -> Red]]