This microsite is associated to the paper

L. Jaulin (2023). Modelisation of a rolling disk with Sympy

Code to generate the evolution function of the rolling disk

from sympy import *


def Symb_build():
    t = symbols('t')
    m,g,r= symbols('m g r')
    c1,c2 = Function('c1')('t'),Function('c2')('t')
    dc1,dc2 = Function('dc1')('t'),Function('dc2')('t')
    ddc1,ddc2 = Function('ddc1')('t'),Function('ddc2')('t'),
    ψ,θ,φ = Function('ψ')('t'),Function('θ')('t'),Function('φ')('t')
    dψ,dθ,dφ = Function('dψ')('t'),Function('dθ')('t'),Function('dφ')('t')
    ddψ,ddθ,ddφ = Function('ddψ')('t'),Function('ddθ')('t'),Function('ddφ')('t')
    λ1,λ2 = Function('λ1')('t'),Function('λ2')('t')


    def subsdiff(E):
        E=E.subs({diff(dc1,t):ddc1,diff(dc2,t):ddc2,  diff(dφ,t):ddφ, diff(dθ,t):ddθ, diff(dψ,t):ddψ})
        E=E.subs({diff(c1,t):dc1,  diff(c2,t):dc2,    diff(φ,t):dφ,   diff(θ,t):dθ,   diff(ψ,t):dψ})
        return simplify(E)

    def Reuler(φ,θ,ψ):
        Rφ = Matrix([ [1,0,0],[0,cos(φ),-sin(φ)],[0,sin(φ),cos(φ)]])
        Rθ = Matrix([ [cos(θ),0,sin(θ)],[0,1,0],[-sin(θ),0,cos(θ)]])
        Rψ = Matrix([ [cos(ψ),-sin(ψ),0],[sin(ψ),cos(ψ),0],[0,0,1]])
        return Rψ*Rθ*Rφ

    def wr(R):
        W=Transpose(R)*diff(R,t)
        return Matrix([[-W[1,2]],[W[0,2]],[-W[0,1]]])

    def Lagrangian(q,dq):
        c1,c2,φ,θ,ψ=list(q)
        c=Matrix([[c1],[c2],[r*cos(θ)]])
        dc=diff(c,t)
        R=Reuler(φ,θ,ψ)
        Ep=m*g*c[2]
        I=Matrix([[1/2*m*r**2,0,0],[0,1/4*m*r**2,0],[0,0,1/4*m*r**2]])
        Ek=1/2*m*(dc.dot(dc))+(1/2)*wr(R).dot(I*wr(R))
        L=subsdiff(Ek-Ep)
        return Matrix([L])


    q=Matrix([c1,c2,φ,θ,ψ])
    dq=Matrix([dc1,dc2,dφ,dθ,dψ])
    ddq=Matrix([ddc1,ddc2,ddφ,ddθ,ddψ])

    L=Lagrangian(q,dq)

    Q=subsdiff(diff(L.jacobian(dq),t)-L.jacobian(q))


    A=Matrix([[1,0,-r*sin(ψ),-r*cos(ψ)*cos(θ),r*sin(ψ)*sin(θ)],
              [0,1, r*cos(ψ),-r*sin(ψ)*cos(θ),-r*cos(ψ)*sin(θ)]])
    τ=λ1*A[0,:]+λ2*A[1,:]
    a=A*dq
    da=diff(a,t)
    da=subsdiff(da)
    da=simplify(da)

    S=Matrix([da,*list(Q-τ)])
    M=S.jacobian([λ1,λ2,ddq])
    b=M*Matrix([λ1,λ2,ddq])-S
    λ1,λ2,ddc1,ddc2,ddφ,ddθ,ddψ=list((M.inv()*b))

    ddφ=simplify(ddφ)
    ddθ=simplify(ddθ)
    ddψ=simplify(ddψ)

    return lambdify((c1,c2,φ,θ,ψ,dφ,dθ,dψ,m,g,r),(dc1-a[0],dc2-a[1],dφ,dθ,dψ,ddφ,ddθ,ddψ))

symb=Symb_build()

Code for the simulation of the rolling disk

from roblib import *


def f(x):
    c1,c2,φ,θ,ψ,dφ,dθ,dψ=list(x.flatten())
    dx=symb(c1,c2,φ,θ,ψ,dφ,dθ,dψ,m,g,r)
    return array([[dx[0]],[dx[1]],[dx[2]],[dx[3]],[dx[4]],[dx[5]],[dx[6]],[dx[7]]])


def draw_wheel3D(ax,x,y,z,φ,θ,ψ,r=1,col='blue',size=1):
    M=tran3H(x,y,z)@eulerH(φ,θ,ψ)@wheel3H(r)
    draw3H(ax,M,col,False,1)


def draw(ax,x,r,col):
    c1,c2,φ,θ,ψ=x[0,0],x[1,0],x[2,0],x[3,0],x[4,0]
    c3=r*cos(θ)
    draw_wheel3D(ax,c1,c2,c3,φ,θ,ψ,r,col)
    draw_axis3D(ax,0,0,0,eye(3,3))

def SimuRollingDisk(x):
    ax = Axes3D(figure())
    dt = 0.02
    draw(ax,x,r,"blue")
    for t in arange(0,10,dt):
        x=x+dt*(0.25*f(x)+0.75*(f(x+(2/3)*dt*f(x))))
        pause(0.01)
        clean3D(ax,-2,14,-8,8,0,16)
        draw(ax,x,r,"black")
        c1,c2,φ,θ,ψ=x[0,0],x[1,0],x[2,0],x[3,0],x[4,0]




m,g,r = 5,9.81,1
c1,c2,φ,θ,ψ,dφ,dθ,dψ= 2,0, 0,0.1,0,2.5,0,0 
x=array([[c1],[c2],[φ],[θ],[ψ],[dφ],[dθ],[dψ]])
SimuRollingDisk(x)