I warmly recommend using
wxMaxima, a free computer algebra system for Windows, Mac OS, Linux, and other operating systems.
I'm very familiar with numerical integration of particle trajectories. I used to do molecular dynamics simulator development using non-quantum mechanical classical potential models and force field models, which model the total force acting on each atom in a simulation, then use discrete time steps and predictor-corrector models to integrate the particle trajectories at very high precision over time.
I like to use \$t\$ for time, \$x\$ for position, \$v\$ for velocity, \$a\$ for acceleration, \$j\$ for jerk or jolt, and \$s\$ for snap.
This is one of my favourite examples:
What kind of acceleration/velocity should a train or elevator have, so that it would feel the smoothest to us humans?
We want everything to be standstill (zero) at time \$t = 0\$, and maximum velocity reached at time \$t = T\$ with zero acceleration left, and maximum acceleration \$A\$ at midpoint during acceleration, \$t = T/2\$. Let's start with a generic cubic polynomial for the jerk:
$$j(t) = J_3 t^3 + J_2 t^2 + J_1 t + J_0$$
This also gives us snap (the rate of change in jerk):
$$s(t) = \frac{d j(t)}{d t} = 3 J_3 t^2 + 2 J_2 t + J_1$$
and acceleration (jerk being the rate of change in acceleration):
$$a(t) = \int_{0}^{t} j(\tau) d\tau = \frac{J_3}{4} t^4 + \frac{J_2}{3} t^3 + \frac{J_1}{2} t^2 + J_0 t$$
Note that because these are definite integrals (integrals over a specific interval), there is no integration constants added.
The requirements we have are
$$\left\lbrace \begin{aligned}
j(0) &= 0 \\
j(T) &= 0 \\
a(0) &= 0 \\
a(T/2) &= A \\
a(T) &= 0 \\
\end{aligned} \right . \quad \iff \quad \left\lbrace \begin{aligned}
J_0 &= 0 \\
J_3 T^3 + J_2 T^2 + J_1 T + J_0 &= 0 \\
0 &= 0 \\
\frac{J_3}{64} T^4 + \frac{J_2}{24} T^3 + \frac{J_1}{8} T^2 + \frac{J_0}{2} T &= A \\
\frac{J_3}{4} T^4 + \frac{J_2}{3} T^3 + \frac{J_1}{2} T^2 + J_0 T &= 0 \\
\end{aligned} \right .$$
This is a system of equations we can solve. (If there was no solution, we'd need to bump \$j(t)\$ one degree higher, until a solution exists.)
The solution is \$J_0 = 0\$, \$J_1 = 32 A / T^2\$, \$J2 = -96 A / T^3\$, and \$J3 = 64 A / T^4\$:
$$\begin{aligned}
s(t) &= \frac{192 A}{T^4} t^2 - \frac{192 A}{T^3} t + \frac{32 A}{T^2} \\
j(t) &= \frac{64 A}{T^4} t^3 - \frac{96 A}{T^3} t^2 + \frac{32 A}{T^2} t \\
a(t) &= \frac{16 A}{T^4} t^4 - \frac{32 A}{T^3} t^3 + \frac{16 A}{T^2} t^2 \\
v(t) &= \frac{16 A}{5 T^4} t^5 - \frac{8 A}{T^3} t^4 + \frac{16 A}{3 T^2} t^3 \\
x(t) &= \frac{8 A}{15 T^4} t^6 - \frac{8 A}{5 T^3} t^5 + \frac{4 A}{3 T^2} t^4 \\
\end{aligned}$$
Let's say that the final velocity (at time \$t = T\$) the velocity is \$V\$:
$$v(T) = V = \frac{8 A T}{15} \quad \iff \quad A = \frac{15 V}{8 T}$$
The distance traveled during acceleration is \$x(T) = 4 A T^2 / 15 = T V / 2\$.
Maximum acceleration or deceleration occurs when its derivative (jerk) is zero; at \$t = T/2\$, when \$a(t) = a(T/2) = A\$. (The two other extrema are \$t = 0\$ and \$t = T\$, but \$a(0) = a(T) = 0\$.)
People are most likely to stumble when the magnitude of jerk is maximum. The extrema of a continuous differentiable function occur when its derivative is zero, i.e. when snap is zero, \$s(t) = 0\$. \$s(t)\$ is a second degree function, so the roots are easily calculated, and are \$t = T/2 \pm T/\sqrt{12}\$, or \$t \approx 0.2113 T\$ and \$t \approx 0.7887 T\$. There, \$j(t) = \pm16 A / (\sqrt{27} T)\$, i.e. maximum magnitude of jerk is \$16 A / (\sqrt{27} T) \approx 3.0792 A/T \approx 5.7735 V/T^2\$.
While the velocity is a fifth-degree polynomial and position six-degree polynomial, finding say "time \$t\$ when the elevator/train has moved distance \$L\$" is relatively simple to do
numerically, when \$T\$ and \$A\$ (or \$V\$) are known: the key is realizing that position \$x(t)\$ and velocity \$v(t)\$ are monotonically increasing functions for \$0 \le t \le T\$, so both
Newton's method and the
bisection method starting with \$t=0\$ and \$t=T\$
will find a single solution \$x(0) \le L \le x(T)\$ and/or \$v(0) \le V \le v(T)\$.
Algebraic solutions may or may not exist, but even when they do, they tend to be annoyingly complicated to calculate.
For example, if \$T = 1\$, \$A = 15/4 = 3.75\$, then \$V = 2\$, and \$x(t) = 2 t^6 - 6 t^5 + 5 t^4\$. \$x(0) = 0\$ and \$x(T) = x(1) = 1\$. Then, \$x(1/2) = 5/32\$, \$x(1/4) = 29/2048\$, \$x(3/4) = 1053\$.
By Newton's method, \$f(t) = x(t) - L = t^4 ( 2 t^2 - 6 t + 5 ) - L\$, \$d f(t) / d t = f'(t) = v(t) = 12 t^5 - 30 t^4 + 20 t^3\$. Iteration is
$$t_{n+1} = t_{n} - \frac{f(t_{n})}{f'(t_{n})} = (10 t^6 - 24 t^5 + 15 t^4 + L) / (12 t^5 - 30 t^4 + 20 t^3)$$ and \$t_0 = 1/2\$ is a good starting point for any valid \$L\$ (\$0 \le L \le 1\$).
For example, for \$L = 0.5\$, \$t_1 = 0.84375\$, \$t_2 \approx 0.745861\$, \$t_3 \approx 0.7420736\$, and \$x(0.7420736) \approx 0.500016\$, so one only needs a few iterations to get a pretty good numerical solution. If you clamp each \$t_n\$ to the valid range (\$0 \le t_n \le T\$), it will converge, because both position and velocity are monotonically increasing functions within this range.
Home
Help
Search
Login
Register
