MatLab simulation of a system of masses connected by springs using Euler's method

324 Views Asked by At

I am trying to simulate this one dimensional system of 5 particles of equal mass M and equal diameter D, connected by 4 springs of equal length L and equal elastic constant, K, using the Euler's method. No other forces influence it.
If $0$ is the equilibrium point of the first particle, then $ \forall i = 2,...,5$ the equilibrium positions of the particles should be $x_i ^{eq} = L\cdot (i - 1).$
Now I apply a displacement of the positions (in MatLab, I simply add to each $x_i^{eq}$ a different random value $r_i \in (-L/2 + D/2, L/2 - D/2)$), obtaining new positions $x_1 ^{(1)},...,x_5^{(1)}$.
Then I use the Euler's algorithm applied to the Hooke's law. In Matlab it looks like this:
for i = 2:nrsteps
_F = force(X(:,i - 1),K); _ X(:,i) = X(:,i - 1) + (h * U); _ U = U + (h * F/M);
where h is the time step, U is the velocity vector and force is a function that computes $f_1, ..., f_5$ as follows: $$f_1 = -K\cdot (x_1 - x_2)\\ f_i = -K \cdot (2x_i - x_{i-1} - x_{i+1}), i = 2, 3, 4\\ f_5 = -K \cdot(x_5 - x_4)$$
I have chosen M = K = 1, D = 2 and L = 10 as parameters. In the resulting video though, particles seem to oscillate indipendently from each other (they overlap), just like if each of them was connected to a different spring that is fixed at some margin.
What am I doing wrong? I think that my mistake may be adding the random $r_i$ to each $x_i^{eq}$ without any further assumption other than the range $r_i$ lies in. Any suggestion?

1

There are 1 best solutions below

0
On

The deviation from the rest position of the first spring is $x_2-x_1-L$, not simply $x_2-x_1$. Similarly for the other springs. Thus the forces should be \begin{align} f_1 &= -K·(x_1+L-x_2)\\ f_j &= K·(x_{j-1}+L-x_j) - K·(x_j+L-x_{j+1})=-K·(2x_j-x_{j-1}-x_{j+1})\\ f_5 &= K·(x_4+L-x_5)=-K·(x_5-x_4-L) \end{align}

You should not expect realistic behavior over longer time periods, the error of the Euler method grows initially like $O(t·h)$, later exponentially, so a rather small $h$ and thus a high number of Euler steps are needed to get a useful result.