# Using ParametricNDSolve to solve a hysteretic ODE (Wolfram Mathematica 10.0)

I am trying to solve an ordinary differential equation using ParametricNDSolve and I want to fit the solution to some data to get the value of the parameter. The problem I have is that the ODE is hysteretic meaning that the equation itself is different depending on the value of the boundary condition. I am having trouble implementing this behaviour. Here is the code so far:

First some constants:

b0 = 0.5; bac = 0.001; w = 1; mu0 = 4 Pi*10^-7; jc =
3*10^8; width = 20;


and

bc[d_] := Sqrt[mu0*jc*b0*d]; t0 = 2 Pi/w;


The ODE looks like this:

sol1 = ParametricNDSolve[{D[y[x], {x, 2}] ==
Sign[y[x]]*(1 - Exp[-Sqrt[y[x]^2]]), y[0] == 0,
y'[width/2] == bac/bc[d]*Sin[t0/4]}, y, {x, 0, width/2}, {d},
Method -> {"DiscontinuityProcessing" -> False}]


and it works fine, I get the solution with d as a parameter. The problem is with the next bit of the equation, when the hysteresis comes into play (when the boundary condition is y'[width/2] == bac/bc[d]*Sin[t], t>t0/4).

The equation for t>t0/4 should look something like this:

 ParametricNDSolve[{D[y2[x], {x,2}] == ((2 - Exp[-y[d][x] /. sol1])*
Exp[(y2[x] - y[d][x] /. sol1)/(2 - Exp[-y[d][x] /. sol1])] - 1),y2[0]==0,
y2'[width/2] == bac/bc[d]*Sin[t0/4 + (5*t0)/(4*100)]}, y2, {x, 0, width/2},d]


where y[d][x] /. sol1 should be the solution of the first ODE. But when I execute this I get an error

ParametricNDSolve::dsfun:ParametricFunction ... cannot be used as a function.



I do not know whether ParametricNDSolve does not accept a parametric function in the ODE? Ideally I would like to have sol2 as a parametric function with d as parameter that I can use in FindFit on some data to extract the value of d. Any info would be much appreciated, thanks!

• Is it correct that the integration range of the second step ends at width/2? Jun 1, 2019 at 11:47
• Hi Ulrich, yes the integration stops at width/2. The difference from the first ODE is in the boundary condition for y2' (i. e. y2'[width/2]).
– John
Jun 1, 2019 at 12:55
• Ok. Try to solve the two odes in one step! Jun 1, 2019 at 13:14
• Hi Ulrich, the problem with solving both equations in one step is that the term -y[d][x] /. sol1 in the second equation is the solution of the first equation. Can this be implemented in ParametricNDSolve? Thanks!
– John
Jun 1, 2019 at 13:53

Try to solve it in one step:

ParametricNDSolve[{D[y[x], {x, 2}] ==Sign[y[x]]*(1 - Exp[-Sqrt[y[x]^2]]), y[0] == 0,
y'[width/2] == bac/bc[d]*Sin[t0/4]
, D[y2[x], {x,2}] == ((2 - Exp[-y[x] ])*Exp[(y2[x] - y[x] )/(2 - Exp[-y[x] ])] - 1), y2[0] == 0,
y2'[width/2] == bac/bc[d]*Sin[t0/4 + (5*t0)/(4*100)]
}
, {y, y2}, {x, 0, width/2}, {d},Method -> {"DiscontinuityProcessing" -> False}]

• Hi Ulrich, thank you very much, this solves it! It is a bit embarrassing how simple your solution is! Cheers, John
– John
Jun 1, 2019 at 14:48
• You are Welcome . I like simple solutions... Jun 1, 2019 at 17:35