Best strategies to solve an ODE with piecewise functions using OrdinaryDiffEq

Specific Domains Numerics
, ,
1

I want to solve an ODE problem involving some piecewise functions and I would like to ask for advice on what could be the best strategy to solve them.

Physically, the ODEs I want to solve represent a contact element sticking and sliding according to Coulomb’s friction law and with variable normal load. The code for the forces is the following

using OrdinaryDiffEq, Plots

normal(xn) = xn > 0.0 ? xn : 0.0

function tangential!(w, xt, N)
    w_init = w[1]

    if (N > 0.0) && (abs(xt - w_init) < N)
        T = xt - w_init
        return T

    elseif (N > 0.0) && (abs(xt - w_init) ≥ N)
        sg = sign(xt - w_init)
        @. w = xt - sg * N
        T = sg * N
        return T

    else
        @. w = xt
        T = 0.0
        return T

    end
end

where w is a mutating 1-element vector that takes into account the “memory” of the friction state. Then, the ODE is defined and solved as

function f!(du, u, p, t)
    w = p
    N = normal(u[2])
    T = tangential!(w, u[1], N)

    du[1] = u[3]
    du[2] = u[4]
    du[3] = 0.75sin(t) - T
    du[4] = -N + 0.5
end

begin
    u0 = [0.0, 2.0, 0.0, 0.0]
    p = [0.0]
    tspan = (0, 40π)
    prob = ODEProblem(f!, u0, tspan, p)

    sol = solve(prob, Midpoint())

    plot(sol)
end

where the Midpoint method has been used. The idea is to find the best strategy to solve this simple problem to, in the end, use multiple of these contacts in a more complex setup. I have thought of the following strategies using the OrdinaryDiffEq package:

  1. Low order explicit methods
  2. High order explicit methods using callbacks
  3. Implicit methods

The use of low order explicit methods is easy to implement since the code above just runs perfectly.

High order explicit methods cannot be used directly, since our function does not have continuous derivatives. Therefore, the only way I can see the use of a high order scheme is to add callbacks to the problem. Using callbacks I could stop the integration when we move from one section to the piecewise function to another and, since each section is infinitely differentiable, the use of high order schemes could give good results.

Finally, I found that I could only make implicit methods work by implementing them myself specifically suited for this problem. If I try to execute an implicit scheme with my system I get errors when it tries to automatically differentiate, which I suspect is expected due to the mutations that are going on in my function. Therefore, I do not know if using an implicit method is viable in my problem with the interface of DifferentialEquations and I would like if someone could give some advice on this issue.

2

You probably want to look into the methods for dynamical ODE problems as this is a partitioned (second order) ODE.

Set autodiff=false, like Rodas5P(autodiff=false).

https://diffeq.sciml.ai/stable/basics/faq/#I-get-Dual-number-errors-when-I-solve-my-ODE-with-Rosenbrock-or-SDIRK-methods

Though with your equation I don’t see why one would use an implicit method. Nothing sticks out to me as stiff there.

See https://diffeq.sciml.ai/stable/basics/faq/#My-Problem-Has-Discontinuities-and-is-Unstable-/-Slow or

For your case, you callbacks and a decreased interp_points probably makes sense. Did you try it?

3

I tried doing this, and it works fine, though it gives very different results. I do not know if, when solving the system implicitly, the mutations I have in my system can induce some error by getting some unexpected mutation.

You are right that this equation is not stiff by itself. This system will actually be part of a much larger system trying to model the oscillations of a structure with friction at its contact interfaces. By our experience with previous algorithms in Python, solving the whole system implicitly gave as the fastest integration times, while using IMEX methods had a more restrictive dt so that is why I am interested in the possibility of fully solving the system implicitly.

Thank you for pointing this out. I tried to specify the problem as a dynamical ODE and got a x3 faster solving times using the same numerical scheme. Interestingly, defining it as a SecondOrderODEProblem gives just slightly slower times (0.1 ms seconds more).

I haven’t tried it, since I think it could get a little messy once you have several elements that behave in this piecewise manner. Inspired by your link above, I have tried to find a smooth version of the normal and tangential forces, but it does not seem to be easy since one piecewise function depends on the other. Now I have

normal_smooth(xn) = 1 / 65 * log(1 + exp(65xn))

for the normal force and

function tangential_smooth!(w, xt, N)
    w_init = w[1]
    if (N > 0.0) && (abs(xt - w_init) < N)
        T = xt - w_init
        return T
    elseif (N > 0.0) && (abs(xt - w_init) ≥ N)
        sg = sign(xt - w_init)
        @. w = xt - sg * N
        T = sg * N
        return T
    end
end

for the tangential force, where the last condition of the if-else block was removed but I do not find an easy way to approximate the rest of tangential_smooth! to a better behaved function. Using this I get very close results to my original setup but it probably would be optimal to have this other function as a differentiable function.

4

Sounds like it’s worth trying RKC methods like ROCK2/ROCK4?

Note that Python’s slowness tends to cause a bias towards “more vectorized” methods, giving a bias towards implicit methods simply because they use more linear algebra which has comparatively less overhead than BLAS 1 operations. I wouldn’t extrapolate too far.

1 Like
5

I actually did not know about the existence of this methods, so thank you for pointing it out. I will try them out to see what I can get.