import numpy as np

# parametri fisici
m = 1.0     # kg
g = 9.81    # m/s^2
L = 1.0     # m

def virtual_work(theta, delta_theta):
    """
    Lavoro virtuale del peso per una trave incernierata.
    """
    delta_yG = (L / 2) * np.cos(theta) * delta_theta
    return -m * g * delta_yG

def equilibrium_condition(theta):
    """
    Coefficiente di delta_theta nel lavoro virtuale.
    Deve essere zero per l'equilibrio.
    """
    return m * g * (L / 2) * np.cos(theta)

# test su vari angoli
angles = np.linspace(-np.pi, np.pi, 9)

print("theta [deg] | coefficiente lavoro virtuale")
print("-------------------------------------------")
for th in angles:
    print(f"{th*180/np.pi:8.2f} | {equilibrium_condition(th):8.2f}")

theta [deg] | coefficiente lavoro virtuale
-------------------------------------------
 -180.00 |    -4.91
 -135.00 |    -3.47
  -90.00 |     0.00
  -45.00 |     3.47
    0.00 |     4.91
   45.00 |     3.47
   90.00 |     0.00
  135.00 |    -3.47
  180.00 |    -4.91

import numpy as np

# Parametri fisici
m = 1.0      # massa della trave [kg]
g = 9.81     # gravità [m/s^2]
L = 1.0      # lunghezza della trave [m]
k = 50.0     # costante molla [N/m]

def virtual_work(theta, delta_theta):
    """
    Lavoro virtuale totale della trave con molla al punto destro
    """
    # lavoro del peso al baricentro
    delta_L_peso = - m * g * (L / 2) * np.cos(theta) * delta_theta
    
    # lavoro della molla al punto destro (verticale)
    delta_L_molla = - k * L * np.sin(theta) * L * np.cos(theta) * delta_theta
    
    # lavoro totale
    return delta_L_peso + delta_L_molla

def equilibrium_condition(theta):
    """
    Coefficiente di delta_theta nel lavoro virtuale totale.
    Equilibrio: coefficiente = 0
    """
    return m * g * (L / 2) * np.cos(theta) + k * L**2 * np.sin(theta) * np.cos(theta)

# Risoluzione numerica dell'equilibrio
from scipy.optimize import fsolve

# guess iniziale piccola
theta_guess = 0.0
theta_eq = fsolve(equilibrium_condition, theta_guess)[0]

print(f"Angolo di equilibrio (rad): {theta_eq:.4f}")
print(f"Angolo di equilibrio (deg): {np.degrees(theta_eq):.2f}")

# Verifica: lavoro virtuale a theta_eq
delta_theta_test = 1e-6
lv = virtual_work(theta_eq, delta_theta_test)
print(f"Lavoro virtuale al punto di equilibrio: {lv:.6e}")

Angolo di equilibrio (rad): -0.0983
Angolo di equilibrio (deg): -5.63
Lavoro virtuale al punto di equilibrio: -8.470329e-22

import numpy as np 

def compute_virtual_displacement(partial_derivatives, delta_q):
    """
    Calculate the virtual displacement.
        :param partial_derivatives: Derivatives of the position vector
         with respect to generalized coordinates.
        :param delta_q: Changes in generalized coordinates.
    :return: Virtual displacement.
    """
    return np.dot(partial_derivatives, delta_q)
    
def principle_of_virtual_work(forces, virtual_displacements):
    """
    Compute the virtual work for given forces and virtual displacements.
        :param forces: List or array of forces acting on the system.
        :param virtual_displacements: List or array of virtual displacements.
    :return: Total virtual work.
    """
    return np.dot(forces, virtual_displacements)
    
def apply_lagrange_multipliers(forces, constraint_directions, lambdas):
    """
    Adjust forces with constraint directions via Lagrange multipliers.
        :param forces: Forces acting on the system.
        :param constraint_directions: Directions of constraint forces.
        :param lambdas: Lagrange multipliers.
    :return: Adjusted forces.
    """
    return forces - np.dot(lambdas, constraint_directions)

def determine_equilibrium(config, generalized_forces, constraints):
    """
    Algorithmic implementation to find equilibrium using virtual work.
        :param config: Initial configuration of the system.
        :param generalized_forces: Generalized forces of the system.
        :param constraints: Constraint equations of the system.
    :return: Equilibrium configuration.
    """
    while True:
        delta_q = np.random.rand(len(config)) # Example small changes
        virtual_displacements = compute_virtual_displacement(np.eye(len(config)), delta_q)
        virtual_work = principle_of_virtual_work(generalized_forces, virtual_displacements)

        if np.isclose(virtual_work, 0, atol=1e-5):
            break
        # Simple update example: adjust to achieve zero virtual work
        config -= 0.01 * virtual_work
    return config

# Example setup for testing the implementation
initial_config = np.array([1.0, 1.0, 1.0])
generalized_forces = np.array([5.0, -3.0, 2.0])
constraints = np.array([0.0]) # Placeholder for no constraints
equilibrium_config = determine_equilibrium(initial_config, generalized_forces, constraints)
print("Equilibrium Configuration:", equilibrium_config)

Equilibrium Configuration: [-4822.81810632 -4822.81810632 -4822.81810632]

import numpy as np

# ------------------------------------------------------------
# PARAMETRI FISICI DELLA RUOTA
# ------------------------------------------------------------
m = 1.0        # massa
R = 0.5        # raggio
I = 0.5 * m * R**2   # momento d'inerzia per disco pieno

# ------------------------------------------------------------
# VINCOLO NON OLONOMO
#   dot_x - R * dot_theta = 0
# ------------------------------------------------------------
def nonholonomic_constraint(q_dot):
    """
    Restituisce il valore del vincolo C(q_dot) = dot_x - R dot_theta
    """
    dot_x, dot_theta = q_dot
    return dot_x - R * dot_theta


# ------------------------------------------------------------
# MATRICE DELLE DERIVATE DEL VINCOLO
#   C_qdot = [1,  -R]
#   serve per costruire il termine lambda * C'
# ------------------------------------------------------------
def constraint_jacobian():
    return np.array([1.0, -R])   # derivata del vincolo rispetto alle velocità


# ------------------------------------------------------------
# EQUAZIONI DI LAGRANGE CON MOLTIPLICATORE
#   m ddot_x =  lambda
#   I ddot_theta = - R lambda
#
#   scritte in forma matriciale:
#       M * q_ddot = lambda * C'
# ------------------------------------------------------------
def mass_matrix():
    return np.diag([m, I])


# ------------------------------------------------------------
# RISOLUZIONE DELL'EQUILIBRIO ISTANTANEO
# (per accelerazioni che soddisfano il vincolo)
#
#   Risolviamo il sistema:
#   [ M   -C'^T ] [ q_ddot ]   = [ 0 ]
#   [ C'    0   ] [ lambda ]     [ 0 ]
#
# È il sistema classico per vincoli non olonomi con Lagrange multipliers
# ------------------------------------------------------------
def solve_accelerations():
    M = mass_matrix()
    C = constraint_jacobian()

    # Costruzione della matrice aumentata:
    # [ M   -C^T ]
    # [ C     0  ]
    A = np.block([
        [M,           -C.reshape(-1,1)],
        [C.reshape(1,-1), np.zeros((1,1))]
    ])

    # Termine di destra (zero perché nessuna forza esterna)
    b = np.zeros(3)

    # Soluzione del sistema lineare
    solution = np.linalg.solve(A, b)

    q_ddot = solution[:2]
    lambda_value = solution[2]

    return q_ddot, lambda_value


# ------------------------------------------------------------
# INTEGRAZIONE TEMPORALE (EULERO)
# ------------------------------------------------------------
def simulate(dt=0.01, T=2.0):
    q = np.array([0.0, 0.0])        # [x, theta]
    q_dot = np.array([1.0, 1.0/R])  # velocità iniziale CONSISTENTE col vincolo

    xs = []
    thetas = []

    for _ in np.arange(0, T, dt):

        # Calcoliamo accelerazioni e lambda
        q_ddot, lam = solve_accelerations()

        # Aggiorniamo velocità e posizione
        q_dot = q_dot + q_ddot * dt
        q = q + q_dot * dt

        xs.append(q[0])
        thetas.append(q[1])

    return np.array(xs), np.array(thetas)


# ------------------------------------------------------------
# ESECUZIONE SIMULAZIONE
# ------------------------------------------------------------
x_values, theta_values = simulate()

print("x final:", x_values[-1])
print("theta final:", theta_values[-1])

x final: 2.0000000000000013
theta final: 4.000000000000003

Principio dei Lavori virtuali¶

Trave rigida incernierata a sinistra, soggetta a peso proprio¶

risulta in equilibrio (instabile) solo quando la trave è in verticale¶

Stessa trave, ma con una molla a destra¶

Esercizio del libro privo di senso¶

Ruota che rotola senza attrito - vincolo NON olonomo¶