Wilbert's website at SocSci

> Physics> Smooth steps

physics/smooth.html 2020-10-31

Smooth steps and minimal jerk

scaled to 0,0 - 1,1:
        x = 10 * t**3 -  15 * t**4  +  6 * t**5
        v = 30 * t**2 -  60 * t**3  + 30 * t**4
        a = 60 * t    - 180 * t**2 + 120 * t**3
        j = 60 * (1  - 6 * t + 6 * t**2)

to get physics units:
t /= duration D
x *= amplitude A
v */ A/D
a */ A/D^2

zero jerk: t = 
        p(x):= 10 * x**3 -  15 * x**4  +  6 * x**5;
        v(x):= 30 * x**2 -  60 * x**3  + 30 * x**4;
        a(x):= 60 * x    - 180 * x**2 + 120 * x**3;
        j(x):= 60 * (1  - 6 * x + 6 * x**2);
solve(j(x), x);
                    3 +/-sqrt(3) 
              x = - -----------
                         6

peak velocity:
v(0.5) = 1.875 A/D
peak acceleration:
a((3-sqrt(3))/6) = 10/sqrt(3) A/D^2 \approx 5.77 A/D^2
peak jerk:
j(0) = 60 A/D^3
j(0.5) = -30 A/D^3
j(1) = 60 A/D^3

lowest order polynomial with finite jerk (sometimes called minimum jerk)
x(t) := a + b*t + c*t**2 + d*t**3 + e*t**4 + f*t**5;
define (v(t), (diff(x(t), t)));
define (a(t), (diff(v(t), t)));
define (j(t), (diff(a(t), t)));
solve([x(0)=0, x(1)=1, v(0)=0, v(1)=0, a(0)=0, a(1)=0], [a, b, c, d, e, f]);
/* [[a = 0, b = 0, c = 0, d = 10, e = - 15, f = 6]] */
p(x) = 10*x**3 + -15*x**4 + 6*x**5


lowest order polynomial with finite jerk, more general case with non-zero v0 and a0
x(t) := a + b*t + c*t**2 + d*t**3 + e*t**4 + f*t**5;
define (v(t), (diff(x(t), t)));
define (a(t), (diff(v(t), t)));
define (j(t), (diff(a(t), t)));
solve([x(0)=0, x(1)=1, v(0)=v0, v(1)=0, a(0)=a0, a(1)=0], [a, b, c, d, e, f]);
/* [[a = 0, b = v0, c = a0/2, d = - 6 v0 -3/2 a0 + 10), e = 8 v0 + 3/2 a0 - 15), f = -3 v0 - 1/2 a0 + 6)]] */
p(x, v0, a0) = v0*x + 0.5*a0*x**2 - (6*v0+1.5*a0-10)*x**3 + (8*v0+1.5*a0-15)*x**4 - (3*v0+0.5*a0-6)*x**5
v(x, v0, a0) = v0 + a0*x - (18*v0+4.5*a0-30)*x**2 + (32*v0+6*a0-60)*x**3 - (15*v0+2.5*a0-30)*x**4
a(x, v0, a0) = a0 - (36*v0+9*a0-60)*x + (96*v0+18*a0-180)*x**2 - (60*v0+10*a0-120)*x**3
j(x, v0, a0) = -36*v0-9*a0+60 + (192*v0+36*a0-360)*x - (180*v0+30*a0-360)*x**2

lowest order polynomial with finite jerk, general case with non-zero v0, v1, a0 and a1
x(t) := a + b*t + c*t**2 + d*t**3 + e*t**4 + f*t**5;
define (v(t), (diff(x(t), t)));
define (a(t), (diff(v(t), t)));
define (j(t), (diff(a(t), t)));
solve([x(0)=0, x(1)=1, v(0)=v0, v(1)=v1, a(0)=a0, a(1)=a1], [a, b, c, d, e, f]);
p=v0*x + 0.5*a0*x**2 - (4*v1+6*v0-0.5*a1+1.5*a0-10)*x**3 + (7*v1+8*v0-a1+1.5*a0-15)x**4 - (3*v1+3*v0-0.5*a1+0.5*a0-6)*x**5
(t,v0,a0,v1,a1) = 
x = (t^4*(14*v1+16*v0-2*a1+3*a0-30))/2-(t^3*(8*v1+12*v0-a1+3*a0-20))/2 -(t^5*(6*v1+6*v0-a1+a0-12))/2+t*v0 + (a0*t^2)/2
v = 2*t^3*(14*v1+16*v0-2*a1+3*a0-30)-(3*t^2*(8*v1+12*v0-a1+3*a0-20))/2 -(5*t^4*(6*v1+6*v0-a1+a0-12))/2+v0+a0*t
a = 6*t^2*(14*v1+16*v0-2*a1+3*a0-30)-3*t*(8*v1+12*v0-a1+3*a0-20) -10*t^3*(6*v1+6*v0-a1+a0-12)+a0
j = 12*t*(14*v1+16*v0-2*a1+3*a0-30)-3*(8*v1+12*v0-a1+3*a0-20) -30*t^2*(6*v1+6*v0-a1+a0-12)

lowest order polynomial with finite acceleration step for start and finite jerk for finish (ease out)
x(t) := a + b*t + c*t**2 + d*t**3 + e*t**4;
define (v(t), (diff(x(t), t)));
define (a(t), (diff(v(t), t)));
define (j(t), (diff(a(t), t)));
solve([x(0)=0, x(1)=1, v(0)=v0, v(1)=v1, a(1)=a1], [a, b, c, d, e]);
x(t) = t*v0 + (-(t^2*(6*v1+6*v0-a1-12))/2)+t^3*(5*v1+3*v0-a1-8)-(t^4*(4*v1+2*v0-a1-6))/2
v(t) = (-t*(6*v1+6*v0-a1-12))+3*t^2*(5*v1+3*v0-a1-8)-2*t^3*(4*v1+2*v0-a1-6)+v0
a(t) = 6*t*(5*v1+3*v0-a1-8)-6*t^2*(4*v1+2*v0-a1-6)-6*v1-6*v0+a1+12
j(t) = 6*(5*v1+3*v0-a1-8)-12*t*(4*v1+2*v0-a1-6)

lowest order polynomial with finite acceleration steps (linear)
x(t) := a + b*t + c*t**2 + d*t**3;
define (v(t), (diff(x(t), t)));
define (a(t), (diff(v(t), t)));
define (j(t), (diff(a(t), t)));
solve([x(0)=0, x(1)=1, v(0)=v0, v(1)=v1], [a, b, c, d]);
x(t) = t*v0 + t^2*(-v1-2*v0+3) + t^3*(v1+v0-2)
v(t) = 3*t^2*(v1+v0-2)+2*t*((-v1)-2*v0+3)+v0
a(t) = 6*t*(v1+v0-2) + 2*(-v1-2*v0+3)
j(t) = 6*(v1+v0-2)

attempt to find higher order polynomial with peak jerk less than 'minimum jerk'
x(t) := a + b*t + c*t**2 + d*t**3 + e*t**4 + f*t**5 + g*t**6;
define (v(t), (diff(x(t), t)));
define (a(t), (diff(v(t), t)));
define (j(t), (diff(a(t), t)));
solve([x(0)=0, x(1)=1, v(0)=0, v(1)=0, a(0)=0, a(1)=0], [a, b, c, d, e, f, g]);
/* [a = 0, b = 0, c = 0, d = 10 - g, e = 3 g - 15, f = 6 - 3 g                                                                   ]]
there is no value for g that reduces the maximum jerk:
j(x, g) = 120*g*x**3 + 60*(6-3*g)*x**2 + 24*(3*g-15) * x + 6*(10-g), j(0,g) = 60-6*g, j(1, g) = 60 + 6*g
*/

even higher order polynomial with peak jerk less than 'minimum jerk' 
x(t) := a + b*t + c*t**2 + d*t**3 + e*t**4 + f*t**5 + g*t**6 + h*t**7;
define (v(t), (diff(x(t), t)));
define (a(t), (diff(v(t), t)));
define (j(t), (diff(a(t), t)));
solve([x(0)=0, x(1)=1, v(0)=0, v(1)=0, a(0)=0, a(1)=0], [a, b, c, d, e, f, g, h]);
[[a = 0, b = 0, c = 0, d = (- g) - 3 h + 10, e = 3 g + 8 h - 15, f = (- 3 g) - 6 h + 6]]
p(x, g, h) = (- g - 3*h + 10)*x**3 + (3*g + 8*h - 15)*x**4 + (- 3*g - 6*h + 6)*x**5 + g*x**6 + h*x**7
j(x, g, h) = 210*h*x**4 + 120*g*x**3 + 60*(-6*h - 3*g + 6)*x**2 + 24*(15*h + 3*g - 15) * x + 6*(- 10*h - g + 10)
j(0, g , h) = 60 - 60*h - 6*g # 10h + g > 0
j(1, g , h) = 60 + (150)*h +6*g # 25h+g < 0
# at t=0 and t=1 this is met by: -10 h < g < -25 h
# for instance g=5, h=-.27 has a peak abs(jerk) < 50 A/D^3 in the entire domain making it less jerky than 'minimum jerk'

#gnuplot
v0=1
a0=1
v1=0
a1=0

set key top left
set yrange [-4:6]

set multiplot layout 2,2 rowsfirst

set label 1 'infinite jerk' at graph 0.8,0.9 font ',8'
p1(x) = x*v0 + x**2*(-v1-2*v0+3) + x**3*(v1+v0-2)
v1(x) = 3*x**2*(v1+v0-2)+2*x*((-v1)-2*v0+3)+v0
a1(x) = 6*x*(v1+v0-2) + 2*(-v1-2*v0+3)
j1(x) = 6*(v1+v0-2)
plot [0:1] p1(x) lw 2, v1(x), a1(x)


set label 1 'ease out' at graph 0.8,0.9 font ',8'
p2(x) = x*v0 + (-(x**2*(6*v1+6*v0-a1-12))/2)+x**3*(5*v1+3*v0-a1-8)-(x**4*(4*v1+2*v0-a1-6))/2
v2(x) = (-x*(6*v1+6*v0-a1-12))+3*x**2*(5*v1+3*v0-a1-8)-2*x**3*(4*v1+2*v0-a1-6)+v0
a2(x) = 6*x*(5*v1+3*v0-a1-8)-6*x**2*(4*v1+2*v0-a1-6)-6*v1-6*v0+a1+12
j2(x) = 6*(5*v1+3*v0-a1-8)-12*x*(4*v1+2*v0-a1-6)
plot [0:1] p2(x) lw 2, v2(x), a2(x)


set label 1 'x' at graph 0.8,0.9 font ',8'
p3(x) = (x**4*(14*v1+16*v0-2*a1+3*a0-30))/2-(x**3*(8*v1+12*v0-a1+3*a0-20))/2 -(x**5*(6*v1+6*v0-a1+a0-12))/2+x*v0 + (a0*x**2)/2
v3(x) = 2*x**3*(14*v1+16*v0-2*a1+3*a0-30)-(3*x**2*(8*v1+12*v0-a1+3*a0-20))/2 -(5*x**4*(6*v1+6*v0-a1+a0-12))/2+v0+a0*x
a3(x) = 6*x**2*(14*v1+16*v0-2*a1+3*a0-30)-3*x*(8*v1+12*v0-a1+3*a0-20) -10*x**3*(6*v1+6*v0-a1+a0-12)+a0
j3(x) = 12*x*(14*v1+16*v0-2*a1+3*a0-30)-3*(8*v1+12*v0-a1+3*a0-20) -30*x**2*(6*v1+6*v0-a1+a0-12)
plot [0:1] p3(x) lw 2, v3(x), a3(x)

plot [0:1][0:1] p1(x) lw 2, p2(x) lw 2, p3(x) lw 2


unset multiplot