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.
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.
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 ))