Simulations of complementarity systems using time-stepping

One of the useful outcomes of the theory of complementarity systems is a new family of methods for numerical simulation of discontinuous systems. Here we will demonstrate the essence by introducing the method of time-stepping. And we do it by means of an example.

Example 1 (Simulation using time-stepping) Consider the following discontinuous dynamical system in \mathbb R^2: \begin{aligned} \dot x_1 &= -\operatorname{sign} x_1 + 2 \operatorname{sign} x_2\\ \dot x_2 &= -2\operatorname{sign} x_1 -\operatorname{sign} x_2. \end{aligned}

The state portrait is in Fig. 1.

Show the code 
f₁(x) = -sign(x[1]) + 2*sign(x[2])
f₂(x) = -2*sign(x[1]) - sign(x[2])
f(x) = [f₁(x), f₂(x)]

using CairoMakie
fig = Figure(size = (800, 800),fontsize=20)
ax = Axis(fig[1, 1], xlabel = "x₁", ylabel = "x₂")
streamplot!(ax,(x₁,x₂)->Point2f(f([x₁,x₂])), -2.0..2.0, -2.0..2.0, colormap = :magma)
fig
Precompiling CairoMakie...
    481.2 ms  ✓ Graphics
    685.0 ms  ✓ FreeType
   1317.4 ms  ✓ Cairo
   1494.7 ms  ✓ FreeTypeAbstraction
   3137.2 ms  ✓ MathTeXEngine
  27960.0 ms  ✓ TiffImages
  66238.1 ms  ✓ Makie
  24447.9 ms  ✓ CairoMakie
  8 dependencies successfully precompiled in 120 seconds. 267 already precompiled.
Figure 1: State portrait of the discontinuous system

One particular (vector) state trajectory obtained by some default ODE solver is in Fig. 2.

Show the code 
using DifferentialEquations

function f!(dx,x,p,t)
    dx[1] = -sign(x[1]) + 2*sign(x[2])
    dx[2] = -2*sign(x[1]) - sign(x[2])
end

x0 = [-1.0, 1.0]
tfinal = 2.0
tspan = (0.0,tfinal)
prob = ODEProblem(f!,x0,tspan)
sol = solve(prob)

using Plots
Plots.plot(sol,xlabel="t",ylabel="x",label=false,lw=3)
Precompiling DifferentialEquations...
    505.6 ms  ✓ FiniteDiff
    530.9 ms  ✓ FiniteDiff → FiniteDiffSparseArraysExt
    506.3 ms  ✓ DifferentiationInterface → DifferentiationInterfaceFiniteDiffExt
    620.8 ms  ✓ FiniteDiff → FiniteDiffStaticArraysExt
    798.0 ms  ✓ FiniteDiff → FiniteDiffBandedMatricesExt
   1870.2 ms  ✓ ForwardDiff
    619.1 ms  ✓ RecursiveArrayTools → RecursiveArrayToolsForwardDiffExt
    633.7 ms  ✓ FastPower → FastPowerForwardDiffExt
    695.0 ms  ✓ PreallocationTools
    728.4 ms  ✓ DifferentiationInterface → DifferentiationInterfaceForwardDiffExt
    857.8 ms  ✓ ForwardDiff → ForwardDiffStaticArraysExt
    757.8 ms  ✓ NLSolversBase
    880.7 ms  ✓ BracketingNonlinearSolve → BracketingNonlinearSolveForwardDiffExt
    830.9 ms  ✓ LoopVectorization → ForwardDiffExt
    967.1 ms  ✓ NonlinearSolveBase → NonlinearSolveBaseForwardDiffExt
   1364.7 ms  ✓ LineSearches
   2342.9 ms  ✓ SparseDiffTools
   1015.5 ms  ✓ LineSearch → LineSearchLineSearchesExt
   1037.4 ms  ✓ NLsolve
   3505.2 ms  ✓ DiffEqBase
   1162.1 ms  ✓ SparseDiffTools → SparseDiffToolsPolyesterExt
   1102.6 ms  ✓ DiffEqBase → DiffEqBaseSparseArraysExt
   1139.2 ms  ✓ DiffEqBase → DiffEqBaseChainRulesCoreExt
   1993.4 ms  ✓ Optim
   1240.1 ms  ✓ DiffEqBase → DiffEqBaseDistributionsExt
    980.1 ms  ✓ NonlinearSolveBase → NonlinearSolveBaseDiffEqBaseExt
   2328.6 ms  ✓ DiffEqCallbacks
   2788.2 ms  ✓ OrdinaryDiffEqCore
   7509.8 ms  ✓ SimpleNonlinearSolve
   2320.5 ms  ✓ JumpProcesses
   2941.1 ms  ✓ DiffEqNoiseProcess
   1943.9 ms  ✓ SteadyStateDiffEq
   4953.1 ms  ✓ NonlinearSolveSpectralMethods
   1317.6 ms  ✓ OrdinaryDiffEqCore → OrdinaryDiffEqCoreEnzymeCoreExt
   1239.2 ms  ✓ SimpleNonlinearSolve → SimpleNonlinearSolveChainRulesCoreExt
   1187.5 ms  ✓ SimpleNonlinearSolve → SimpleNonlinearSolveDiffEqBaseExt
   1284.5 ms  ✓ NonlinearSolveSpectralMethods → NonlinearSolveSpectralMethodsForwardDiffExt
   9928.9 ms  ✓ Sundials
   1431.9 ms  ✓ OrdinaryDiffEqFunctionMap
   1612.9 ms  ✓ OrdinaryDiffEqPRK
   1742.2 ms  ✓ OrdinaryDiffEqExplicitRK
   9520.9 ms  ✓ NonlinearSolveQuasiNewton
   2250.3 ms  ✓ OrdinaryDiffEqStabilizedRK
   4801.9 ms  ✓ OrdinaryDiffEqTsit5
   2106.0 ms  ✓ OrdinaryDiffEqRKN
   2540.2 ms  ✓ OrdinaryDiffEqSSPRK
   1600.8 ms  ✓ OrdinaryDiffEqQPRK
   1826.6 ms  ✓ OrdinaryDiffEqSymplecticRK
   2109.4 ms  ✓ OrdinaryDiffEqHighOrderRK
   3072.3 ms  ✓ OrdinaryDiffEqLowOrderRK
   1988.2 ms  ✓ OrdinaryDiffEqFeagin
  15594.6 ms  ✓ NonlinearSolveFirstOrder
   2871.5 ms  ✓ OrdinaryDiffEqLowStorageRK
   1646.2 ms  ✓ OrdinaryDiffEqNordsieck
   2441.7 ms  ✓ NonlinearSolveQuasiNewton → NonlinearSolveQuasiNewtonForwardDiffExt
   3110.0 ms  ✓ OrdinaryDiffEqDifferentiation
   1991.4 ms  ✓ OrdinaryDiffEqAdamsBashforthMoulton
   4113.8 ms  ✓ BoundaryValueDiffEqCore
   4695.1 ms  ✓ OrdinaryDiffEqExtrapolation
   4251.7 ms  ✓ BoundaryValueDiffEqAscher
   9335.4 ms  ✓ NonlinearSolve
   5125.7 ms  ✓ BoundaryValueDiffEqMIRKN
  16299.2 ms  ✓ OrdinaryDiffEqRosenbrock
   3603.7 ms  ✓ NonlinearSolve → NonlinearSolveNLsolveExt
   3459.8 ms  ✓ NonlinearSolve → NonlinearSolveSundialsExt
   3924.5 ms  ✓ OrdinaryDiffEqNonlinearSolve
  29659.8 ms  ✓ OrdinaryDiffEqVerner
   3640.4 ms  ✓ OrdinaryDiffEqStabilizedIRK
   3954.6 ms  ✓ OrdinaryDiffEqPDIRK
   4964.9 ms  ✓ OrdinaryDiffEqSDIRK
   3596.8 ms  ✓ OrdinaryDiffEqIMEXMultistep
   1869.2 ms  ✓ OrdinaryDiffEqLinear
   3945.7 ms  ✓ OrdinaryDiffEqExponentialRK
   7492.2 ms  ✓ StochasticDiffEq
  13569.4 ms  ✓ OrdinaryDiffEqFIRK
   9873.5 ms  ✓ OrdinaryDiffEqBDF
  23604.3 ms  ✓ OrdinaryDiffEqDefault
  63691.8 ms  ✓ BoundaryValueDiffEqMIRK
   7942.9 ms  ✓ OrdinaryDiffEq
   4460.0 ms  ✓ BoundaryValueDiffEqShooting
   4699.9 ms  ✓ DelayDiffEq
  91590.7 ms  ✓ BoundaryValueDiffEqFIRK
   6814.6 ms  ✓ BoundaryValueDiffEq
   7975.9 ms  ✓ DifferentialEquations
  84 dependencies successfully precompiled in 142 seconds. 237 already precompiled.
Precompiling SciMLBaseMakieExt...
   4208.2 ms  ✓ SciMLBase → SciMLBaseMakieExt
  1 dependency successfully precompiled in 5 seconds. 305 already precompiled.
Precompiling DiffEqBaseUnitfulExt...
    955.6 ms  ✓ DiffEqBase → DiffEqBaseUnitfulExt
  1 dependency successfully precompiled in 2 seconds. 123 already precompiled.
Precompiling FileIOExt...
   1683.3 ms  ✓ Plots → FileIOExt
  1 dependency successfully precompiled in 2 seconds. 191 already precompiled.
Precompiling GeometryBasicsExt...
   1777.6 ms  ✓ Plots → GeometryBasicsExt
  1 dependency successfully precompiled in 2 seconds. 207 already precompiled.
Figure 2: Trajectory of the discontinuous system

We can also plot the trajectory in the state space, as in Fig. 3.

Show the code 
Plots.plot(sol[1,:],sol[2,:],xlabel="x₁",ylabel="x₂",label=false,aspect_ratio=:equal,lw=3,xlims=(-1.2,0.5))