2d cable based on segments

I'm trying to simulate in Python a 2d cable based on massless constant segments. It's something similar to a chain of pendulums. I would like to control one endnode with the mouse, while the rest of the chain is free. I got following code, but in the simulation, the whole cable falls freely and cannot hold first node with the mouse. Can anyone guess what I did wrong, please?

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

# Simulation parameters
num_nodes = 10
length = 1.0
dt = 0.01
total_time = 10.0
gravity = 9.81
segment_length = length / (num_nodes - 1)
constraint_iterations = 10

# Initialize positions and velocities
node_positions = np.zeros((num_nodes, 2))
node_velocities = np.zeros((num_nodes, 2))

# Set initial positions (straight line)
for i in range(num_nodes):
    node_positions[i, 0] = i * segment_length

# Function to update node positions
def update_positions(positions, velocities, dt, control_position):
    new_positions = positions.copy()
    new_velocities = velocities.copy()

    # Apply gravity to velocities
    new_velocities[1:, 1] -= gravity * dt  # Skip the controlled endpoint

    # Update positions with velocities
    new_positions[1:] += new_velocities[1:] * dt  # Skip the controlled endpoint

    # Set the controlled endpoint to the control_position
    new_positions[0] = control_position

    # Constraint iterations to maintain constant segment lengths
    for _ in range(constraint_iterations):
        for i in range(1, num_nodes):
            segment_vector = new_positions[i] - new_positions[i - 1]
            segment_length_current = np.linalg.norm(segment_vector)
            if segment_length_current != 0:
                correction_vector = (segment_length_current - segment_length) * (segment_vector / segment_length_current)
                if i < num_nodes - 1:  # Adjust both nodes if not the last segment
                    new_positions[i] -= correction_vector / 2
                    new_positions[i - 1] += correction_vector / 2
                else:  # Adjust only the previous node if the last segment
                    new_positions[i - 1] += correction_vector

    return new_positions, new_velocities

# Function to handle mouse motion
def on_mouse_move(event):
    if event.inaxes:
        global control_position
        control_position = np.array([event.xdata, event.ydata])

# Initialize control position
control_position = node_positions[0]

# Set up plot
fig, ax = plt.subplots()
ax.set_xlim(-0.1, length + 0.1)
ax.set_ylim(-length, length)
line, = ax.plot(node_positions[:, 0], node_positions[:, 1], 'o-', color='blue')

# Update function for animation
def animate(frame):
    global node_positions, node_velocities
    node_positions, node_velocities = update_positions(node_positions, node_velocities, dt, control_position)
    line.set_data(node_positions[:, 0], node_positions[:, 1])
    return line,

# Connect the mouse motion event
fig.canvas.mpl_connect('motion_notify_event', on_mouse_move)

# Create animation
ani = FuncAnimation(fig, animate, frames=int(total_time / dt), interval=dt*1000, blit=True)

plt.xlabel('X Position')
plt.ylabel('Y Position')
plt.title('2D Cable Simulation with Cursor Control')
plt.show()

It is not really a coding problem. There is no problem with your code (at least, no coding problem that could explain why it seems not attached to your mouse. Of course, one could debate about global, or the absence of vectorization of the constraint). I mean, no problem with mouse event listener. The 1st need position is indeed attached to the mouse position.

The problem is more in the maths/modelisation of constraint.

There are 3 problems in the way constraints are implemented.

1. Last node free falls

Firstly, as is, the last node is never altered by constraints. Which means that it falls freely, whatever occurs with the rest of the chain. Its position is segment_length*(num_nodes-1) - 0.5*0.91*t**2 (free fall from its initial position. Computed using the integration scheme of the lines new_velocities[1:, 1] -= gravity * dt and new_positions[1:] += new_velocities[1:] * dt

Plus, it doesn't really make sense: why, for the last segment, only the before last node should be impacted?

If one segment should be different, it should be the first, maybe, but no reason to have a special behaviour for the last segment. Especially if there isn't for the first (the first in not only the symmetric of the last, so even if there were a physical reason to have a special treatment for the last — and I don't think there is —, that reason should apply for the first also. And even more, the first segment at least have the reason that it is the one attached to the mouse) I feel that that special treatment comes from confusion between segments and nodes. That with some way of doing the same thing, you may be tempted to iterate nodes, apply some correction due to attachment with previous and with next node, but for last segment, that hasn't a next node. But here, each segment is just a segment, with no relation with previous or next segment (except the implicit one that they share a node). And you already cover the fact that last node hasn't a next node, by iterating segments from 1 to num_nodes (so 1 segment less than nodes).

That last node alone, falling freely toward y=-∞, was dragging all the others with it, because of the constraints. Since the constraint apply dissymetrically with, that means that the before last node is dragged +1/2 force upward to try to get closer to the before-before-last node. But dragged -1 downward to try to get closer to the last one, that free falls toward y=-∞. Plus, that +1/2 coefficient are not absolute: they apply to a correctionVector that way longer for the last segment (since it is trying to catch up with a free falling node), than with other segments (that have correctionVector symmetric). So, before last node goes to -∞ too. And therefore so does before-before-last node. And therefore so does before-before-before-last node. Etc.

So

       for i in range(1, num_nodes):
            segment_vector = new_positions[i] - new_positions[i - 1]
            segment_length_current = np.linalg.norm(segment_vector)
            if segment_length_current != 0:
                correction_vector = (segment_length_current - segment_length) * (segment_vector / segment_length_current)
                new_positions[i] -= correction_vector / 2
                new_positions[i - 1] += correction_vector / 2

2. First node goes not strictly at mouse position

If you get rid of that special treatment, you have something more moderate. But still, the chain falls downward. At least it remains a chain, with the constraints. But it falls downward. Because the adjustment for the constraint is a compromise. If you consider 1st segment, you move downward the 1st node position to do half of the extra-distance with the 2nd node, and the 2nd node upward to do the other half. But that is not what it should do. The 1st node should not be engaged in any compromise. its position is non-negociable.

So, it is the 1st segment that should have a special treatment: only the 2nd node move upward, the 1st node doesn't go down.

        for i in range(1, num_nodes):
            segment_vector = new_positions[i] - new_positions[i - 1]
            segment_length_current = np.linalg.norm(segment_vector)
            if segment_length_current != 0:
                correction_vector = (segment_length_current - segment_length) * (segment_vector / segment_length_current)
                if i==1:
                    new_positions[i] -= correction_vector
                else:
                    new_positions[i] -= correction_vector / 2
                    new_positions[i - 1] += correction_vector / 2

3. Ever increasing speed.

With this, you get something where no segment is freefalling. The first segment remain attached to the mouse. And the other remain linked with it. Yet, they fall. They don't free fall, but they fall. And the segments get longer and longer.

Why, because what your constraint (with my corrections) do probably what you intended to do

1. Compute speed and position as if each node were free falling
2. Correct node position to enforce length constraint
3. Because each node (but the first) has some contradictory correction to do, one step of this is not enough. So, you redo step 2 some times, until it converges toward positions where all constraints are respected.

10 is apparently not enough. This is anyway a Zeno situation. You have an error (length longer than they should) proportional to speed*dt. And then each step reduce that error, until it is small enough. And 10 seems not enough.

Reason why I say "seems", is because if you do 100 steps instead, it seems to work perfectly.... for a while. And one may wonder why a scheme that works well when mouse is at position, say, (0.5, 0.5) at the beginning (meaning that the 10 corrections are sufficient so that the length error are negligible from the error created at each time steps), are not sufficient after a while (after 2 minutes or so, clearly, even if you start still from a position (0.5, 0.5), correction are not sufficient.

Reason for that is that you haven't corrected speed.

So, at the beginning of simulation, you compute a quite slow speed because of the free fall. That slow speed put the nodes a little be too low, but 10 corrections of their position are enough to put each of them where they belong. But the speed is untouched. So it continues to accelerate as if it were freefalling.

After 2 minutes, speed is g*120 (so MACH 3, roughly). And you don't see it, because the corrections prevent the nodes to move. But at each timestep, what occurs, is, each node is moved downward a few dozens meters, because of their MACH 3 speed. And then moved up a few dozens meters, because of the 10 corrections.

But when speed is really uncanny, the 10 corrections are not enough. And if you do 1000 corrections steps, it will take 10 minutes before it is isn't enough, but well, if at each timestep each nodes goes downward 10 miles, then upward 9.99 miles, it starts to really show.

So, long story short, you need to take into account that speed is also broken by your corrections (as in real life after all. If your cable was a literal chain, not only each link of the chain is forced to keep its length, but also, it speeds is stopped by the same phenomenon (reality is more that each link is elastic, and that correction vector applies only to speed, not to position. But that requires a very small dt unless the cable is really elastic)

Since speed is derivative of the position, that means correcting speed by a factor 1/dt (it is quite strange, since that means derivating something that is the result of an integration. but, well, it is what comes with that hard-enforced constraint)

    # Constraint iterations to maintain constant segment lengths
    for _ in range(constraint_iterations):
        for i in range(1, num_nodes):
            segment_vector = new_positions[i] - new_positions[i - 1]
            segment_length_current = np.linalg.norm(segment_vector)
            if segment_length_current != 0:
                correction_vector = (segment_length_current - segment_length) * (segment_vector / segment_length_current)
                if i==1:
                    new_positions[i] -= correction_vector
                    new_velocities[i] -= correction_vector/dt
                else:
                    new_positions[i] -= correction_vector / 2
                    new_velocities[i] -= correction_vector/dt/2
                    new_positions[i - 1] += correction_vector / 2
                    new_velocities[i-1] += correction_vector/dt/2

There is some numerical error. But it is stable. And even if you have an insufficient number of correction steps, the consequence is only an apparent temporay elasticity. But nothing cumulative.

4. damping

I could add a fourth point. That there is no friction anywhere here. So you get a perpetual motion (for example, if you make a pendulum motion)

But firstly, there is a an implicit involuntary one: triangular inequality implies a small loss of energy during correction of non straight-aligned segments.

And secondly, I guess that would have been the next step. (quite simple to do: speed is corrected by a factor <1)