This microsite is associated to the paper
L. Jaulin (2024), Asymptotically minimal interval contractors based on the centered form; Application to the stability analysis of linear time-delayed differential equations, Acta Cybernetica.

This paper proposed a contractor which is minimal when dealing with narrow boxes. The method is based on the centered form classically used by interval algorithms combined with a Gauss Jordan band diagonalisation preconditioning. As an illustration, we propose to draw the set of all parameter vectors of a characteristic equation of a dynamical system which have a zero.

A video with an illustration of the convergence

Resolution with classic interval propagation

from codac import *
from vibes import vibes
from collections import deque
import numpy as np
import time

def C_classic(X):
    while 1:
        X0=X
        C.contract(X)
        if (X.is_empty()) or (X0.max_diam()/X.max_diam()<1.01): return X

def SIVIA(C,X0, eps,col):
  stack =  deque([X0])
  lf = LargestFirst(eps/1.0)
  k = 0
  while len(stack) > 0:
    X = stack.popleft()
    X=C(X)
    if not(X.is_empty()):
        if (X.max_diam() < eps):
          vibes.drawBox(X[0].lb(),X[0].ub(),X[1].lb(),X[1].ub(),col)
        else:
            (X1,X2) = lf.bisect(X)
            stack.append(X1)
            stack.append(X2)
  vibes.axisEqual()

f1 = Function("p1", "p2", "w", "-w^2+2*w*sin(w*p1)+cos(w*p2)")
f2 = Function("p1", "p2", "w", "2*w*cos(w*p1)-sin(w*p2)")
C1=CtcFwdBwd(f1,Interval(0,0))
C2=CtcFwdBwd(f2,Interval(0,0))
C=C1&C2
vibes.beginDrawing()
vibes.setFigureSize(1000,1000)

X0=[[0,2],[2,4],[0,10]]
SIVIA(C_classic,IntervalVector(X0),  2**(-8),'blue[cyan]')
X=IntervalVector(X0)
vibes.drawBox(X[0].lb(),X[0].ub(),X[1].lb(),X[1].ub(), "yellow[transparent]" )
vibes.endDrawing()

Resolution with the new centered contractor

from codac import *
from vibes import vibes
from collections import deque
import numpy as np
from numpy.linalg import inv,det
from scipy.linalg import lu
import time



def generate_J(F1x,F2x):
    F = sp.Matrix([F1x,F2x])
    Ax = F.jacobian(x)
    Ax=sp.simplify(Ax)
    A = Ax.subs({x[i]: z[i] for i in range(3)})
    return sp2str(A[0,0]),sp2str(A[0,1]),sp2str(A[0,2]),sp2str(A[1,0]),sp2str(A[1,1]),sp2str(A[1,2])


def sp2str(E):
    F= str(E.func)
    if isinstance(E,sp.Integer) or isinstance(E, sp.Symbol): X=str(E)
    elif F=="":
        X="("+sp2str(E.args[0])+")"
        for i in range(1,len(E.args)):              X=X+"+("+sp2str(E.args[i])+")"
    elif F=="":
        X="("+sp2str(E.args[0])+")"
        for i in range(1,len(E.args)):              X=X+"*("+sp2str(E.args[i])+")"
    elif F=="" and str(E.args[1])=="2":
        X="sqr("+sp2str(E.args[0])+")"
    elif F=="":       X="("+sp2str(E.args[0])+")^("+sp2str(E.args[1])+")"
    elif F=="":     X="-("+sp2str(E.args[0])+")"
    elif F in {"sin",'cos','exp'}:                  X=F+"("+sp2str(E.args[0])+")"
    else: print("ERROR sp2str(E) : F=",F)
    return X


x = [sp.symbols(f'x[{i}]') for i in range(3)]
z = [sp.Symbol(f'z[{i}]') for i in range(3)]
za = [sp.Symbol(f'za[{i}]') for i in range(3)]
zb = [sp.Symbol(f'zb[{i}]') for i in range(3)]
m = [sp.Symbol(f'm[{i}]') for i in range(3)]



def generate_Gp(E1,E2):
    q00,q01,q10,q11 = sp.symbols('q00 q01 q10 q11')
    phip1,phip2 = sp.symbols('phi1 phi2')
    F = sp.Matrix([E1,E2])
    Ax = F.jacobian(x)
    Ax=sp.simplify(Ax)
    A = Ax.subs({x[i]: z[i] for i in range(3)})
    a11,a12,a13,a21,a22,a23=A[0,0],A[0,1],A[0,2],A[1,0],A[1,1],A[1,2]
    ap11=q00*a11+q01*a21
    ap12=q00*a12+q01*a22
    ap13=q00*a13+q01*a23
    ap21=q10*a11+q11*a21
    ap22=q10*a12+q11*a22
    ap23=q10*a13+q11*a23
    g1 =phip1+ap11*(x[0]-m[0])+ap12*(x[1]-m[1])+ap13*(x[2]-m[2])
    g2 =phip2+ap21*(x[0]-m[0])+ap22*(x[1]-m[1])+ap23*(x[2]-m[2])
    g1=g1.subs({z[i]:za[i] for i in range(3)})
    g2=g2.subs({z[i]:zb[i] for i in range(3)})
    return sp2str(g1),sp2str(g2)





def SIVIA(C,X0, eps,col):
    vibes.beginDrawing()
    vibes.setFigureSize(1000,1000)
    stack =  deque([X0])
    k = 0
    while len(stack) > 0:
      X = stack.popleft()
      X=C(X)
      if not(X.is_empty()):
          if (X.max_diam()