Springs

This implementation is using ipywidgets to add some interactivity, so please use the binder link if you want to try it out, otherwise keep reading and enjoy!

#collapse
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy.integrate import solve_ivp
from ipywidgets import interactive

Defining a system with 3 masses connected with 2 springs:

N=3
xs0 = np.linspace(0,1,N)
ys0 = np.zeros(N)
k=1.0  # spring coefficient
m=4.0  # masses

def plot_system(ax,xs=[],ys=[]):
    ln, = ax.plot(xs, ys, '-ro', ms=20, markerfacecolor='black', markeredgecolor='black')
    return ln

fig,ax=plt.subplots()
plot_system(ax, xs=xs0, ys=ys0);

Mesure the original length of the springs:

lx0 = np.diff(xs0)
ly0 = np.diff(ys0)

Define a method to calculate the spring forces in x, and y direction:

def calculate_spring_forces(xs,ys):
    dxs = np.diff(xs)
    dys = np.diff(ys)
    fxs = k*(dxs-lx0)
    fys = k*(dys-ly0)
    return fxs,fys

Define a method to transfer the spring forces to the connected masses. A spring create a force that acts on both masses that it is connected to in the oposite direction:

def spring_force_to_point(fs):
    
    fs_ = np.concatenate([[0],fs,[0]])
    f_1 = fs_[0:-1]
    f_2 = fs_[1:]
    return f_1, f_2

def calculate_point_forces(fxs,fys):
    
    fx_left, fx_right = spring_force_to_point(fxs)
    fy_bottom, fy_top = spring_force_to_point(fys)
    
    fx = -fx_left + fx_right
    fy = -fy_bottom + fy_top
    
    return fx,fy
    

Here is a plot showing the spring force acting on the masses:

#collapse
def f(x0=-0.3,y0=0.5,x1=0.5,y1=0.0,x2=1.3,y2=-0.0):
    
    fig,ax=plt.subplots()
    xs = [x0,x1,x2]
    ys = [y0,y1,y2]
    
    fxs,fys = calculate_spring_forces(xs=xs, ys=ys)
    
    fx,fy = calculate_point_forces(fxs=fxs, fys=fys)
    
    for x,y,fx_,fy_ in zip(xs,ys,fx,fy):
        ax.arrow(x,y,fx_,fy_, head_width=0.05, head_length=0.1,)
    
    plot_system(ax,xs=xs,ys=ys)
    
    
    ax.set_xlim(-0.5,1.5)
    ax.set_ylim(-1,1)
    plt.show()
    ax.set_title('Point forces')

    
interactive_plot = interactive(f, x0=(-0.5,1.0), x1=(-0.5,1.0), x2=(-0.5,1.0), y0=(-1.0,1.0), y1=(-1.0,1.0), y2=(-1.0,1.0))
output = interactive_plot.children[-1]
output.layout.height = '350px'
interactive_plot

This system will now be simulated, where the following state space model will be used:

def update(t, states):
    
    xs = states[0:N]
    ys = states[N:2*N]
    
    dxs = states[2*N:3*N]
    dys = states[3*N:]
    
    
    fxs,fys = calculate_spring_forces(xs=xs, ys=ys)
    fx,fy = calculate_point_forces(fxs=fxs, fys=fys)
    
    ddxs = fx/m
    ddys = fy/m
    
    dstates = np.concatenate([dxs,dys,ddxs,ddys])
    
    return dstates

A simulation where the masses have been moved from their original positions will now be conducted:

px = np.random.rand(N)
py = np.random.rand(N)

xlim=[-0.5,1.5]
ylim=[-1,1]

xs0 = xlim[0] + px*np.diff(xlim)
ys0 = ylim[0] + py*np.diff(ylim)

dxs0 = np.zeros(N)
dys0 = np.zeros(N)

states_0 = np.concatenate([xs0,ys0,dxs0,dys0])


interval=1/25
seconds=40
Ns=int(40/interval)
t = np.linspace(0,40,Ns)
result = solve_ivp(fun = update,t_span=[t[0],t[-1]],t_eval=t, y0 = states_0)

The animation function in matplotlib can be used to visualize the solution:

#hide_output
#collapse
from matplotlib.animation import FuncAnimation
from matplotlib import animation, rc
from IPython.display import HTML

fig, ax = plt.subplots()
ln, = ax.plot([], [], '-ro', ms=20, markerfacecolor='black', markeredgecolor='black')

def init():
    ax.set_xlim(-0.5,1.5)
    ax.set_ylim(-1,1)
    return ln,

def update_plot(frame):
    
    xs = result.y[0:N,frame]
    ys = result.y[N:2*N,frame]
        
    ln.set_data(xs, ys)
    return ln,

anim = FuncAnimation(fig, update_plot, frames=np.arange(Ns),
                    init_func=init, blit=True, interval=interval)
#anim

And there is the video:

#collapse
plt.rcParams["animation.html"] = "jshtml"
HTML(anim.to_jshtml())

</input>

Once Loop Reflect