# 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,
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,
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? – Ulrich Neumann Jun 1 '19 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 '19 at 12:55
• Ok. Try to solve the two odes in one step! – Ulrich Neumann Jun 1 '19 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 '19 at 13:53

ParametricNDSolve[{D[y[x], {x, 2}] ==Sign[y[x]]*(1 - Exp[-Sqrt[y[x]^2]]), y == 0,