7. Elastic string model#
7.1. Forces on a string#
Imagine a thin cord stretched between a fixed support and your hand. Move the held end and let go: it moves, changes shape, and eventually settles. Which forces make that happen?
Gravity pulls every small piece downward.
Tension is the pulling force transmitted along the string. A stretched piece tries to shorten, pulling on the pieces beside it. It cannot push along its length.
Bending rigidity is the stiffness resisting a change of curvature.
Internal friction dissipates energy.
Supports exert reaction forces that hold the attached ends in place.
Length constraints keep the represented pieces of material from collapsing into points.
If a piece has more material than the distance it spans, it becomes slack rather than acting like a compressed metal spring. Its material length remains present in bends and folds. We leave out contact, knots, and breaking.
7.2. Tension#
It’s helpful to think of elastic string as a spring. In 1676 Robert Hooke found through experimentation, that the force needed to extend a spring by some distance scales linearly with respect to that distance and called it Hooke’s law. Let’s define tension \(T\) as the force that needs to be applied to a string to extend it to length \(\ell\). Conversly, this the force with which string pulls its sides to restore original length:
Where:
\(\ell_0\) is original string length;
\(\ell\) is its current length;
\(\Delta\ell\) is string extension;
\(k\) is a material’s property, called spring constant with units N/m.
We can expand \(k\) as a property of string’s cross-sectional area, multiplied by Youngs’ modulus:
Where:
\(E\) is Young’s modulus, a material’s resistance to axial strain, in N/m²;
\(A_0\) is string’s cross-sectional area, in m².
The product \(EA_0\) is the axial rigidity, in newtons.
The expanded equation for tension:
The function \(\max\) chooses the larger argument, so negative strain produces zero tension.
Real material becomes thinner when pulled, but in this equation \(A_0\) is constant, which suitable only for small extension. And that’s what allows us to think about it as a spring.
7.3. Contour length and slack#
A 1 m piece can have its ends 10 cm apart without losing 90 cm of material. Its contour length is the distance along the string, following its bends, rather than the straight distance between its ends.
Our discrete model represents that contour by many short links. Each link may stretch, but we impose a lower bound at its manufactured length:
Bringing the supports closer must therefore bend the chain instead of shortening its links to nothing. This is an additional modelling assumption, not a consequence of Hooke’s law. It treats each resolved link as unable to buckle internally; folds must be represented by several links.
The tensile spring only pulls. The length constraint can also supply a separating reaction when a link would become too short. That reaction is not an elastic compression spring: no compression energy is stored and returned. Both its mathematical rule and the numerical corrections that approximate it are given below.
7.4. Internal friction#
As a string deforms, some of its mechanical energy becomes heat. This internal friction makes oscillations gradually die down. Elastic tension and bending rigidity alone would return the energy they store.
We can model resistance to changing length by adding a force proportional to how fast the length changes. For a taut piece, the pulling force needed to stretch it becomes
Here \(c\geq0\) is a damping coefficient in N s/m, and \(d\ell/dt\) is the rate of length change in m/s. For a uniform piece, \(c=\eta A_0/\ell_0\), where \(\eta\) is the material’s viscosity in N s/m². Stretching increases the required pull; shortening reduces it. For a tension-only string, clamp the total to zero if this formula gives a negative value. The spring and damping terms acting together are called the Kelvin–Voigt model.
This friction depends on deformation: moving the whole string without changing its shape or keeping it at a fixed stretch produces no internal damping.
7.5. Bending rigidity#
Bend a string without stretching it: the material on the outside of the bend extends, while the material on the inside shortens. Their resistance creates a turning effect, called a bending moment, that tries to straighten the string. A thin thread bends easily; a thick cord made of the same material is harder to bend.
For a string that is straight when unloaded, a simple model is
Where:
\(M_b\) is the bending moment needed to hold the bend, in N m;
\(\kappa\) is curvature, in 1/m: a circular bend of radius \(R\) has curvature \(1/R\);
\(B\) is bending rigidity, in N m²; bending rigidity stores elastic energy, just as stretching does;
\(E\) is Young’s modulus, the same material stiffness used for tension, in N/m²;
\(I\) describes how the cross-sectional area is distributed around the bending axis, in m⁴.
For a solid circular string of radius \(R\), \(I=\pi R^4/4\). Doubling its radius therefore multiplies bending rigidity by sixteen. This assumes a solid elastic material; fibers sliding inside a braided cord can make it bend more easily.
Friction can also resist changes of bending, even when the string’s length barely changes. For bending in a plane, add a term proportional to the rate of curvature change:
Here \(\kappa\) is signed curvature, distinguishing the two bending directions, and its time derivative follows the same piece of material. The coefficient \(D\geq0\) has units N m² s. For the same uniform Kelvin–Voigt material used above, \(D=\eta I\). The extra moment opposes changing the bend; holding a fixed bend produces no damping. Trubat and Molins: Rheological Damping of Slender Rods derives both the extension and bending forms.
Extension damping and bending damping describe different changes of shape. Either can be used on its own as a simplified model, or both can act together. The point-chain simulation below includes both elastic effects and both kinds of damping.
7.6. Try the string#
One end is fixed. Grab the orange end, pull it sideways, up, or down, then release it. The hand stops supporting that end, and the string falls and swings under gravity. The fixed end remains attached. Use Reset string to start again or Let go to release without dragging. Mouse and touch both work. With the scene focused, Space grabs or releases, arrow keys move the held end, and Escape releases it. The demonstration string is 4 m long. Keep holding the orange end and move it up and down to send a disturbance toward the fixed end. Slow motion runs at quarter speed so the motion along the string is easier to follow. Normal speed is the default. Switch it off for normal-speed playback; this changes playback speed, not the physical forces.
You can also open the string experiment on its own.
The demo and Python example below use the same 80-link, 4 m model and the same update: a semi-implicit Euler update followed by contour-length corrections, with a step of 0.0001 s. The left end is fixed; the other is prescribed only while held. Release sets that endpoint’s velocity to zero; the rest of the string keeps its current motion. The hand approaches the pointer smoothly at up to 6 m/s of physical time and reaches at most 4.6 m from the fixed end. Rendering is separate from physical ticks; slow devices may run behind real time.
The forces and constraint reactions are specified below and included in both implementations. There is no air resistance, ground contact, or self-contact.
7.7. A continuous mathematical model#
Label each material point by its distance \(s\) along the original string, with \(0\leq s\leq L_0\). Here \(L_0\) is the total manufactured length in metres. This label stays with the material as it moves. Let \(\vec{r}(s,t)\) be its current position in metres at time \(t\) in seconds. The string can move in two or three dimensions.
The derivative \(\partial\vec{r}/\partial s\) describes how position changes between neighboring material labels. Its length is the local ratio of current length to original length. Define that ratio as \(\lambda\), and the unit vector along the string as \(\vec{\tau}\):
Both \(\lambda\) and \(\vec{\tau}\) are dimensionless. The direction requires \(\lambda>0\); a collapsed piece has no defined tangent. The tension includes extension damping, with zero force in slack material. The time derivative of \(\lambda\) is taken at a fixed material label \(s\). Let \(\mu\) be mass per original length in kg/m, so total mass is \(\mu L_0\).
The contour-length assumption has an ideal continuous form too. Let \(R\) be the magnitude of its separating axial reaction, in newtons:
A stretched piece has \(R=0\). At its minimum length, a reaction can prevent further shortening. Since \(T\) pulls neighboring pieces together and \(R\) separates them, the signed axial force is \(T-R\).
Consider a short piece between labels \(s\) and \(s+\Delta s\). Its right neighbor exerts \(((T-R)\vec{\tau})(s+\Delta s,t)\); its left neighbor exerts \(-((T-R)\vec{\tau})(s,t)\). Subtract these forces, add gravity and the force from bending, and apply Newton’s law. Dividing by \(\Delta s\) and taking the limit gives
Here \(\vec{g}\) is gravitational acceleration in m/s². The vector \(\vec{f}_b\) is force per original length from bending rigidity and bending damping, in N/m. We compute its discrete counterpart from changes of angle below. Every term has units N/m. The second time derivative is acceleration. The spatial derivative is the difference of neighboring tension vectors per original length. Even equal tension magnitudes can give a net force when their directions differ.
7.8. Boundary and initial conditions#
The equation alone does not specify an experiment. Let \(\vec{a},\vec{b}\) be the fixed support positions in metres. While the hand holds the right end, our boundary conditions are
The support reactions supply whatever forces are necessary to keep those positions fixed. Our supports fix position but allow the string to rotate, so they impose no bending moment. A free end would instead have zero total internal force and zero bending moment, unless an external end force or moment is applied. Releasing the hand switches the right end to this free condition. The demo sets that endpoint’s velocity to zero at release, so release does not throw it; all other nodes retain their velocities. Gravity continues to act.
We must also prescribe an initial shape \(\vec{r}(s,0)\) and initial velocity \(\partial\vec{r}(s,0)/\partial t\) consistent with the supports. Putting supports farther apart than \(L_0\) creates initial stretch, or pretension. Putting them closer together allows a sagging shape without initially stretching every piece. Gravity then changes the shape and the local tensions together.
7.9. Replacing the string with points#
The computer cannot store a position for every real value of \(s\). Instead, split the original string into \(N\) equal segments with rest length \(\ell_0=L_0/N\), joined by \(N+1\) points called nodes. Node \(i\) has position \(\vec{r}_i\), velocity \(\vec{v}_i\), and mass \(m_i\).
Give half of each segment’s mass to each of its endpoints. Interior nodes receive two halves, and the end nodes receive one:
These masses add to \(\mu L_0\). For segment \(i\), joining nodes \(i\) and \(i+1\), compute
The dot product measures the relative velocity along the segment; \(\dot\ell_i\) is its rate of length change in m/s. For a stretched segment, its elastic and damping contributions combine into
The vector \(\vec{q}_i\), in newtons, is the segment’s force on its left node. Add \(+\vec{q}_i\) to node \(i\) and \(-\vec{q}_i\) to node \(i+1\). Each internal interaction is equal and opposite, so internal forces cannot accelerate the whole string’s center of mass. Slack segments exert zero force. At exactly zero length, set their force to zero without dividing by length.
The force evaluated before constraint reactions is, for an interior node,
Here \(\vec{F}_{b,i}\) is the total bending force on that node, calculated next. The length-constraint reactions are applied separately as impulses after the Euler prediction, so they must not be counted again in \(\vec{F}_i\).
This chain is a spatial approximation of our single string, not a collection of unrelated strings. When doubling the segment count, halve segment masses and double segment spring constants. Double the segment damping coefficients too, since \(c=\eta A_0/\ell_0\). Keeping the old spring constants would make a different, softer string.
7.10. Bending between segments#
In two dimensions, three neighboring nodes form a bend. Let \(\theta_i\) be the signed turn from segment \(\vec{d}_{i-1}\) to segment \(\vec{d}_i\), in radians. It is zero when the segments continue in a straight line. For small extension and well-resolved bends, approximate curvature by \(\theta_i/\ell_0\). For a string whose original shape is straight, choose the bending energy
Each term represents one bend; \(B/\ell_0\) has units N m, so the energy is in joules. To calculate forces, we need to know how moving a node changes that energy. The intermediate step is how moving the node changes the angle.
For bend \(i\), define
This angle gradient is a two-component vector: its components say how much \(\theta_i\) changes per metre of horizontal or vertical movement of node \(j\). For a small displacement \(\delta\vec{r}_j\), with the other nodes held still,
We can calculate these vectors directly from the two segment vectors. For a planar vector \(\vec{d}=(d_x,d_y)\), define
This is the gradient of the segment’s direction angle with respect to its components. Since the turning angle is the outgoing direction minus the incoming direction,
The same angle-change rule gives the turning rate when all three nodes move:
Now use the gradient to differentiate the bend’s energy with respect to a node’s position:
This is the elastic force on that node from this bend. Adding bending damping gives the moment and its corresponding nodal force:
Thus each bend supplies three force vectors, one for each of its nodes. A node can belong to several bends; add their contributions to obtain the force used in the rate function:
In the Python implementation below, gradients contains these three
\(\vec{G}_{i,j}\) vectors, theta_rate is their velocity dot-product sum,
and moment is \(m_{b,i}\).
The line force[i - 1:i + 2] -= moment * gradients adds the bend’s three
contributions to the force array. Dividing that accumulated array by node masses
then gives the bending contribution to acceleration.
The damping contribution removes mechanical energy at rate \(D\dot\theta_i^2/\ell_0\).
The three angle gradients add to zero, so each bend’s forces sum to zero. Rigid translation or rotation changes neither the turning angles nor their rates. Bending therefore acts on shape changes, including in slack material. When refining the chain, keep \(B\) and \(D\) fixed and recompute the coefficients \(B/\ell_0\) and \(D/\ell_0\). The formulas require nonzero segment lengths and a resolved bend rather than an exact 180-degree reversal.
7.11. State and rate function#
In Basics of simulation software, the state combined one body’s position and velocity. For the string, collect the positions and velocities of all nodes:
The vector \(\vec{s}\) is the simulation state; the scalar \(s\) used earlier labels material along the string. Masses, rest lengths, and material constants are fixed parameters, so they do not need to be part of the evolving state.
The rate function \(f(\vec{s})\) computes the derivatives from the supplied state. For each free node,
Compute tension, bending, and damping forces from the supplied positions and velocities, add gravity, then divide each node’s total force by its mass. For the permanently fixed left end, return zero position and velocity derivatives. While the right end is held, its position derivative is the prescribed hand velocity and its acceleration is zero for that tick. Once released, it uses the same free-node rates as the interior nodes. The held/free flag and hand target are boundary data supplied to the rate function and tick.
Evaluating \(f\) does not advance time or change its input. It only answers: how fast would this state change?
7.12. Advancing the simulation#
We still use a first-order Euler method, but evaluate internal friction at the new velocity. Strong internal friction makes a fully explicit step require impractically short time steps. Solving the friction implicitly avoids that restriction without adding a new damping force. This changes how we advance the state, not the rate function or its physical forces.
Separate the state into all positions \(\vec{r}\) and all velocities \(\vec{v}\), and write \(\vec{a}(\vec{r},\vec{v})\) for the acceleration part of \(f\). At each tick solve
The segment directions and angle gradients come from the old positions. The viscous forces use the solved velocities, and the positions advance with those velocities. This is a semi-implicit, first-order Euler update rather than the fully explicit version in the basics chapter. The prescribed support velocities are boundary conditions on the solve. The complete tick then applies the contour correction described next:
7.13. Solving for velocity#
Every material contribution has the form
Here \(\vec{G}_a\) contains a deformation gradient at the affected nodes and zeros elsewhere. For extension it contains \((-\vec{n}_i,\vec{n}_i)\), with \(e_a=EA_0(\ell_i-\ell_0)/\ell_0\) and \(c_a=\eta A_0/\ell_0\). For bending it contains the three angle gradients already derived, with \(e_a=B\theta_i/\ell_0\) and \(c_a=D/\ell_0\). The extension contribution is included only for a stretched link whose total tension is positive.
Let \(\mathsf M\) be the diagonal mass matrix, with each node’s mass repeated for its two coordinates. For a fixed choice of taut links, the velocity equation becomes a linear system:
Solve this system, then check the tensions with the new velocities. If a total tension would be negative, omit that link’s contribution and solve again; if a previously omitted stretched link now has positive tension, include it again. The implementation uses Newton iterations with backtracking when the taut/slack choice changes. Bending contributions remain active in either direction.
The matrix is symmetric and positive definite on the free coordinates.
Each term couples at most three neighboring nodes, so only five diagonals below the main diagonal need storage.
solve_band uses this structure for a Cholesky solve without building a dense matrix.
The Python solve_material_velocity function and JavaScript solver perform the same update.
The JavaScript version reuses its working arrays and specializes the Cholesky
calculation for this fixed bandwidth. Bending contributions are assembled once
per substep and reused while the solver adjusts which tensile terms are active.
Implicit integration for stiff internal forces is also discussed in
Baraff and Witkin: Large Steps in Cloth Simulation.
The nominal tick is at most 0.0001 s. Fast motion subdivides it further. For maximum node speed \(v_{\max}\), propose a step no larger than \(0.9\,(0.01\ell_0)/v_{\max}\); at zero speed use the remaining tick duration. The 10% margin leaves room for acceleration and contour corrections, avoiding repeated solves of steps that would predictably exceed the movement allowance. Retry with half that step if any node, including its contour correction, moves more than \(0.01\ell_0\). This limits geometric changes during rapid drags and tight folds. It is temporal refinement, not a force or velocity damping term.
The interactive demo gives physics about 8 ms per animation frame, checking the budget after each accepted substep. An ongoing solve finishes before yielding. If rapid motion requires more work, the demo discards overdue wall time and runs more slowly while keeping the same integration steps and material coefficients. This prevents extra substeps from blocking drawing and pointer input.
The integration and contour corrections can introduce energy error. Halving the step and disabling material friction are useful separate checks: a fast decay should come from the chosen friction coefficients, rather than disappear when the numerical step is refined.
7.14. Correcting shortened links#
For a predicted link, compute its current length \(\ell_i\) and direction \(\vec{n}_i=(\vec{r}_{i+1}-\vec{r}_i)/\ell_i\). The required separation increase is
Let \(\omega_i=1/m_i\) be inverse mass for a free node, in 1/kg. Set \(\omega_i=0\) at a fixed or held node, so a constraint cannot move a prescribed point. If \(c_i>0\) and the two inverse masses do not both vanish, apply
These corrections increase the pair’s distance by \(c_i\). Equal masses move equally; a lighter node moves farther; a prescribed node does not move. For two free nodes, the mass-weighted corrections sum to zero.
Repair the positions separately from the velocities. Turning a position repair into \(\Delta\vec{r}/h\) can add kinetic energy and make compressed links rebound. Instead, for a link at its minimum length, measure its relative closing velocity:
Here \(J_i\) is a separating impulse in N s. Apply
This cancels closing motion without creating a separating rebound. For two free nodes the impulse conserves momentum and cannot increase kinetic energy. A moving prescribed endpoint can still do work on the string. This inelastic constraint rule is a numerical modelling choice, separate from material viscosity. It stores no elastic compression energy.
A correction can shorten a neighboring link, so visit all links twice: first from the hand toward the fixed end, then in the opposite direction. Recompute lengths and directions after each correction. The code treats links within \(10^{-6}\ell_0\) above their minimum as touching, repairs any shortage, and cancels closing velocity for those pairs. Two sweeps leave a finite residual; this is an approximate constraint solve, not an exact projection of the whole chain. Position repair can still change potential energy, so check convergence by reducing the time step.
For the demonstration, \(B=10^{-4}\) N m² and \(D=10^{-5}\) N m² s. Here bending friction is an independent effective coefficient, rather than assuming a single uniform Kelvin–Voigt material. Its relaxation time \(D/B=0.1\) s and the axial value \(\eta A_0/(EA_0)=0.3\) s are chosen independently. This damps small bending ripples while leaving the resistance to a static bend unchanged. The smaller bending stiffness lets a lifted hanging end form a local fold; a much stiffer strand resists this and transmits a stronger push upward. The 80 links (81 nodes) give a rest spacing of 5 cm. This reduces solver work and the number of adaptive substeps needed during fast motion. Tight folds have less geometric detail than with 160 links: in a narrow sag, neighboring segments can turn by about 34 degrees rather than 17 degrees. The material coefficients stay the same; masses and discrete forces scale with the rest spacing. An exactly straight planar chain also has no preferred buckling direction; small existing bends select the side on which it folds. Finite bending stiffness and changes in tension can still disturb the upper string.
7.15. Moving and releasing the hand#
While the right end is held, limit the target to 4.6 m from the fixed support. Let \(\vec{b}_{\mathrm{target}}\) be that bounded target, and \(\vec{d}_H=\vec{b}_{\mathrm{target}}-\vec{r}_N\) its displacement from the held endpoint. Pointer events arrive as separate position samples. Reaching each sample immediately would repeatedly start and stop the hand, sending short pulses into the string. Instead, use a critically damped controller to approach the target smoothly. Before evaluating \(f\), calculate
Limit the acceleration magnitude to \(a_{\max}=100\) m/s², then update the prescribed velocity with Euler and limit its speed:
Leave a zero vector at zero when applying either magnitude limit. This controller describes the imposed hand motion, not a force or damping term in the string. The Euler prediction moves the held point at the resulting prescribed velocity; constraint corrections leave it alone.
On release, set only \(\vec{v}_N\) to zero and remove its prescribed status. It then has its ordinary endpoint mass and responds to gravity, internal forces, damping, and constraint impulses. The permanent support remains fixed.
7.16. Starting with a held string#
Reset starts with a straight, unstretched string and zero velocities. Both supports are initially at the same height, a distance \(L_0\) apart. Gravity pulls the interior nodes down, producing extension, tension, and a sagging shape. Internal extension and bending friction dissipate the resulting oscillations. No equilibrium shape is precomputed or imposed.
The demonstration uses illustrative coefficients \(\eta A_0=30\) N s and \(D=10^{-5}\) N m² s so the straight-start motion damps visibly within a few simulated seconds. These are effective material parameters, not an external drag term or measured properties of a particular cord. They also damp small, rapid ripples more strongly than broad waves.
7.17. Python implementation#
This is the same 4 m, 80-link model used by the demo: tension, extension damping, angle-based bending and bending damping, and the contour-length constraint. The coefficients are illustrative; extension and bending use independent effective damping coefficients.
"""The article's Python version of the interactive string's physical tick."""
import numpy as np
N = 80
L0, mu0, EA = 4.0, 0.02, 100.0 # m, kg/m, N
etaA, h = 30.0, 0.0001 # N s, s
B = 0.0001 # N m^2
D = 0.00001 # N m^2 s: effective bending friction
ell0 = L0 / N
g = np.array([0.0, -9.81])
anchor = np.array([-0.4 * L0, 0.65 * L0])
hand_reach, hand_speed = 1.15 * L0, 1.5 * L0
hand_response, hand_acceleration = 0.015, 100.0 # s, m/s^2
mass = np.full(N + 1, mu0 * ell0)
mass[[0, -1]] *= 0.5
def initial_state():
r = np.column_stack((anchor[0] + L0 * np.arange(N + 1) / N,
np.full(N + 1, anchor[1])))
return np.array([r, np.zeros_like(r)])
def f(state, held=True):
r, v = state
force = mass[:, None] * g
for i in range(N):
d = r[i + 1] - r[i]
length = np.linalg.norm(d)
if length > ell0:
direction = d / length
length_rate = np.dot(v[i + 1] - v[i], direction)
tension = max((EA * (length - ell0)
+ etaA * length_rate) / ell0, 0.0)
q = tension * direction
force[i] += q
force[i + 1] -= q
for i in range(1, N):
a, b = r[i] - r[i - 1], r[i + 1] - r[i]
aa, bb = np.dot(a, a), np.dot(b, b)
if min(aa, bb) < 1e-24:
continue # No resolved direction at exactly zero length.
theta = np.arctan2(a[0]*b[1] - a[1]*b[0], np.dot(a, b))
pa = np.array([-a[1], a[0]]) / aa
pb = np.array([-b[1], b[0]]) / bb
gradients = np.array([pa, -pa - pb, pb])
theta_rate = np.sum(gradients * v[i - 1:i + 2])
moment = (B * theta + D * theta_rate) / ell0
force[i - 1:i + 2] -= moment * gradients
dr = v.copy()
dv = force / mass[:, None]
dr[0] = dv[0] = 0.0
if held:
dv[-1] = 0.0 # dr[-1] remains the prescribed hand velocity.
return np.array([dr, dv])
def material_terms(r):
"""Frozen geometric gradients, elastic forces, and viscous coefficients."""
terms = []
for i in range(N):
d = r[i + 1] - r[i]
length = np.linalg.norm(d)
if length > ell0:
direction = d / length
terms.append((2*i, np.concatenate((-direction, direction)),
EA*(length - ell0)/ell0, etaA/ell0, True))
for i in range(1, N):
a, b = r[i] - r[i - 1], r[i + 1] - r[i]
aa, bb = np.dot(a, a), np.dot(b, b)
if min(aa, bb) < 1e-24:
continue
theta = np.arctan2(a[0]*b[1] - a[1]*b[0], np.dot(a, b))
pa = np.array([-a[1], a[0]]) / aa
pb = np.array([-b[1], b[0]]) / bb
gradients = np.concatenate((pa, -pa - pb, pb))
terms.append((2*(i - 1), gradients, B*theta/ell0, D/ell0, False))
return terms
def solve_band(matrix, rhs):
"""Cholesky solve; matrix[i, k] stores the entry (i, i-k)."""
size, stride = matrix.shape
width = stride - 1
for i in range(size):
for j in range(max(0, i - width), i + 1):
value = matrix[i, i - j]
for k in range(max(0, i - width), j):
value -= matrix[i, i - k] * matrix[j, j - k]
matrix[i, i - j] = (np.sqrt(value) if i == j
else value / matrix[j, 0])
result = rhs.copy()
for i in range(size):
for j in range(max(0, i - width), i):
result[i] -= matrix[i, i - j] * result[j]
result[i] /= matrix[i, 0]
for i in range(size - 1, -1, -1):
for j in range(i + 1, min(size, i + width + 1)):
result[i] -= matrix[j, j - i] * result[j]
result[i] /= matrix[i, 0]
return result
def solve_material_velocity(state, h, held):
"""Solve v1 = v0 + h*a(r0, v1), with the supports prescribed."""
terms = material_terms(state[0])
old = state[1].ravel().copy()
old[:2] = 0
size = old.size
masses = np.repeat(mass, 2)
known = np.zeros(size, dtype=bool)
known[:2] = True
if held:
known[-2:] = True
def rate(term, v):
start, gradient, _, _, _ = term
return np.dot(gradient, v[start:start + gradient.size])
def active(term, v):
_, _, elastic, damping, unilateral = term
return not unilateral or damping == 0 or elastic + damping*rate(term, v) > 0
def objective(v):
value = np.sum(0.5*masses*(v - old)**2 - h*masses*np.tile(g, N + 1)*v)
for term in terms:
_, _, elastic, damping, unilateral = term
speed = rate(term, v)
value += h * (max(0.0, elastic + damping*speed)**2/(2*damping)
if unilateral and damping > 0
else elastic*speed + 0.5*damping*speed**2)
return value
current = old.copy()
for _ in range(40):
matrix = np.zeros((size, 6))
matrix[:, 0] = masses
rhs = masses * (old + h*np.tile(g, N + 1))
flags = [active(term, current) for term in terms]
for enabled, term in zip(flags, terms):
if not enabled:
continue
start, gradient, elastic, damping, _ = term
for a, ga in enumerate(gradient):
i = start + a
rhs[i] -= h*elastic*ga
for b in range(a + 1):
matrix[i, a - b] += h*damping*ga*gradient[b]
for i in range(size):
for j in range(max(0, i - 5), i):
if known[i]:
rhs[j] -= matrix[i, i - j]*old[i]
matrix[i, i - j] = 0
elif known[j]:
rhs[i] -= matrix[i, i - j]*old[j]
matrix[i, i - j] = 0
matrix[known, 0] = 1
rhs[known] = old[known]
candidate = solve_band(matrix, rhs)
if all(flag == active(term, candidate) for flag, term in zip(flags, terms)):
return candidate.reshape(N + 1, 2)
cost, fraction = objective(current), 1.0
for _ in range(24):
trial = current + fraction*(candidate - current)
if objective(trial) <= cost + 1e-15:
break
fraction *= 0.5
if np.max(np.abs(trial - current)) < 1e-10:
return trial.reshape(N + 1, 2)
current = trial
raise RuntimeError("Internal friction solve did not converge")
def project_contour(state, held):
r, v = state
inverse_mass = 1 / mass
inverse_mass[0] = 0.0
if held:
inverse_mass[-1] = 0.0
for indices in (range(N - 1, -1, -1), range(N)):
for i in indices:
d = r[i + 1] - r[i]
length = np.linalg.norm(d)
if length > ell0 * (1 + 1e-6):
continue
wi, wj = inverse_mass[i:i + 2]
if wi + wj == 0:
continue
direction = d / length if length > 1e-12 else np.array([1.0, 0.0])
correction = max(0.0, ell0 - length) / (wi + wj)
left = -wi * correction * direction
right = wj * correction * direction
r[i] += left
r[i + 1] += right
# Cancel closing motion; position repair does not supply a rebound.
closing = np.dot(v[i + 1] - v[i], direction)
impulse = max(0.0, -closing) / (wi + wj)
v[i] -= wi * impulse * direction
v[i + 1] += wj * impulse * direction
def substep(state, h, held=True, target=None):
start = state.copy()
if held:
point = start[0, -1].copy() if target is None else np.array(target)
offset = point - anchor
point = anchor + offset * min(1.0, hand_reach / max(np.linalg.norm(offset), 1e-12))
delta = point - start[0, -1]
omega = 1 / hand_response
acceleration = omega**2 * delta - 2 * omega * start[1, -1]
acceleration *= min(1.0, hand_acceleration / max(np.linalg.norm(acceleration), 1e-12))
velocity = start[1, -1] + h * acceleration
start[1, -1] = velocity * min(1.0, hand_speed / max(np.linalg.norm(velocity), 1e-12))
prediction = start.copy()
prediction[1] = solve_material_velocity(start, h, held)
prediction[0] += h * prediction[1]
project_contour(prediction, held)
return prediction
def tick(state, h, held=True, target=None):
result = state.copy()
if held and target is None:
target = result[0, -1].copy()
limit, remaining = 0.01 * ell0, h
while remaining > h * 1e-12:
speed = np.max(np.linalg.norm(result[1], axis=1))
# Leave room for acceleration and contour corrections before retrying.
part = min(remaining, 0.9 * limit / max(speed, 1e-12))
for attempt in range(21):
prediction = substep(result, part, held, target)
movement = np.max(np.linalg.norm(prediction[0] - result[0], axis=1))
if movement <= limit * (1 + 1e-6):
break
part *= 0.5
else:
raise RuntimeError("Could not resolve the string motion")
result = prediction
remaining -= part
return result
def release(state):
released = state.copy()
released[1, -1] = 0.0
return released
if __name__ == "__main__":
state = release(initial_state())
for _ in range(round(0.1 / h)): # 0.1 physical seconds after letting go.
state = tick(state, h, held=False)
r, v = state
print("Free endpoint after 0.1 s:", r[-1])
Here state[0] holds all positions and state[1] all velocities.
The rate function returns an array of the same shape containing velocities and accelerations.
tick copies its input, prescribes hand velocity if held, solves the Euler velocity equation, advances positions, and corrects the contour.
release changes only the endpoint velocity; pass held=False on subsequent ticks.
The runnable example releases the initially straight string immediately.
To hold or move the endpoint, call tick(state, h, held=True, target=point);
point is a two-component position in metres. For a fixed hold, pass the same target on every tick; omitting it uses the current endpoint as the next target.
The accepted r values are the points to draw as a connected line.
Rendering and quarter-speed playback do not change the physical tick,
as discussed in Basics of simulation software.
7.18. Choosing a time step and checking the result#
A stiff string has fast vibrations, including lengthwise motion that may be almost invisible. For the uniform, taut chain with small extension, its longitudinal wave speed in original-length coordinates is approximately
The travel time across one segment, \(\ell_0/c_{\mathrm{axial}}\), gives a useful scale for choosing the time step. A fully explicit Euler step can add numerical energy to vibrations, so resolving their period alone is not a stability guarantee. The demo instead solves the velocity-dependent friction implicitly and advances positions with the new velocities. Its geometry is still evaluated at the old positions, so elastic motion and changes of constraint still require temporal resolution. For our numbers, \(c_{\mathrm{axial}}\approx70.7\) m/s and \(\ell_0/c_{\mathrm{axial}}\approx0.000354\) s. The example uses \(h=0.0001\) s and runs for 0.1 s. Check convergence over the intended duration by reducing the time step; this short example does not establish stability for arbitrarily long playback. Refining the chain also requires checking the time step again.
Bending introduces another time scale, roughly \(\ell_0^2\sqrt{\mu/B}\). Shorter bends can vibrate faster, so halving segment lengths reduces this scale by a factor of four. Extension damping adds a time scale of order \(\mu\ell_0^2/(\eta A_0)\), which is about 0.42 microseconds for these coefficients. Increasing viscosity can therefore require a smaller fully explicit step, even though it dissipates physical energy. The implicit friction solve avoids this particular restriction.
Strong bending damping would also restrict a fully explicit time step; its short-wavelength scale is roughly \(\mu\ell_0^4/D\). These are scales rather than universal stability bounds.
This speed differs from the small transverse-wave speed \(c_{\mathrm{transverse}}=\sqrt{T/\mu}\) in OpenStax: Wave Speed on a Stretched String. That formula uses tension and mass per current length, and describes small sideways motion about a taut state. Our extensible model also permits axial vibrations; their faster speed can control the software time step. For nearly inextensible strings, an implicit integrator or a length-constraint solver can be more practical.
Several simple experiments help check the implementation:
With gravity disabled and zero initial velocity, a straight, uniformly stretched string between fixed supports stays at rest. Neighboring tensions cancel at every interior node.
With free endpoints and no external forces, total momentum stays constant up to numerical roundoff.
With fixed supports and internal friction, a displaced string should lose mechanical energy overall. Extension damping dissipates energy only while segment lengths change. Bending damping dissipates energy while turning angles change. The discrete integrator can introduce energy error, so compare runs with smaller time steps.
Bring the supports closer and check that the contour does not disappear into shortened links. Compare the shortest link and the distribution of turning angles as the bottom rounds out.
With gravity disabled, a bent string has restoring forces even if all segments are at their rest lengths. Set \(B=D=0\) to remove those bending forces and compare the motion.
Halve the time step, then refine the segments with their masses and stiffnesses rescaled. Compare sag and oscillation timing at matching physical times.
7.19. From one string to cloth#
A string stores energy when it stretches; a flexible sheet also resists changes of angle between material directions. Crossing many strings suggests warp and weft, but their intersections need a rule for shear, and real cloth needs some bending resistance too. The next chapter, Elastic cloth model, adds these ingredients. The software pattern remains familiar: represent material with finitely many points, compute forces from deformation, apply supports, and advance time.