The problem is to analyze the behavior of an involute of the cubical parabola near the singular points and at infinity.

To do this efficiently, it is best to work with partially evaluated expressions.

We want to investigate the singular points of an involute of the cubical parabola y = x^3.

If the curve is given parametrically by "conjugate(r)(t)", an involute is given by

"conjugate(E)(t)=conjugate(r)(t)-(conjugate(r)'(t))/(|conjugate(r)'(t)|)*(∫)[a]^(t)|conjugate(r)'(tau)| ⅆtau"

What we want to do is to leave the integral unevaluated and compute only its derivatives. The simplest way is to define S := proc (a, t) options operator, arrow; Int(s(tau), tau = a .. t) end proc. Both diff and evalf will be handled automatically by Maple. We'll do it in a slightly more complicated way, leaving S(a, t) completely inert with the exception of the condition S(a, a)=0 and implementing a more efficient numerical evaluator.

forget(diff); forget(series); r := proc (t) options operator, arrow; [t, t^3] end proc; s := unapply(sqrt((diff(r(t)[1], t))^2+(diff(r(t)[2], t))^2), t); S := proc (a, t) options operator, arrow; `if`(t = a, 0, 'S(a, t)') end proc; `diff/S` := proc (at, ft, t) options operator, arrow; s(ft)*(diff(ft, t))-s(at)*(diff(at, t)) end proc; E := unapply(r(t)-`~`[`*`](diff(r(t), t), S(a, t)/s(t)), [a, t])

forget(evalf); `evalf/S` := (proc () local acur, Ssol; acur := undefined; Ssol := subs(dsolve({diff(f(t), t) = s(t), f(a) = 0}, f(t), numeric, parameters = [a], output = listprocedure), f(t)); proc (a, t) if [a, t]::(list(numeric)) then if a <> acur then Ssol(parameters = [a]); acur := a end if; Ssol(t) else 'S(a, t)' end if end proc end proc)()

plot([[op(r(t)), t = -1 .. 3], [op(E(1, t)), t = -1.5 .. 4]], view = -1 .. 3, scaling = constrained)


simplify(diff(E(a, t), t))

[18*S(a, t)*t^3/(9*t^4+1)^(3/2), -6*t*S(a, t)/(9*t^4+1)^(3/2)]


The singular points are easily determined now: we need either S(a, t)=0, which implies t=a, or t=0. We assume a>0.

Compute the Taylor series around t=a first.

`~`[series](E(a, t+a), t = 0, 4); ser := collect(convert(%, polynom), t, normal)

[-18*a^2*(6*a^4-1)*t^3/(9*a^4+1)^2+9*a^3*t^2/(9*a^4+1)+a, 2*(36*a^4-1)*t^3/(9*a^4+1)^2-3*a*t^2/(9*a^4+1)+a^3]


The center of the expansion E(a, a) is the point [a, a^3] on the cubical parabola.

There are two quadratic terms, which give us the slope of the tangent line. The tangent coincides with the normal to the cubical parabola at this point.

Rotate the axes to make the tangent vertical.

rot := proc (e, a) options operator, arrow; convert(Student:-LinearAlgebra:-RotationMatrix(a).Vector(e), list) end proc

rot(ser-E(a, a), (`@`(arctan, op))(`~`[coeff](ser, t, 2))); ser := `assuming`([collect(%, t, normal)], [a > 0])

[-12*a^2*t^3/(9*a^4+1)^(3/2), -2*(162*a^8+9*a^4-1)*t^3/(9*a^4+1)^(5/2)+3*a*t^2/(9*a^4+1)^(1/2)]


In the new coordinates, the leading terms are t^3 and t^2, and the involute has the shape of a semicubical parabola. When t increases, the involute moves from the first to the second quadrant. In terms of x and y, `~`[x]*alpha*y^(3/2), and the coefficient alpha is found from the two leading terms. For t>0:

`assuming`([simplify(coeff(ser[1], t, 3)/coeff(ser[2], t, 2)^(3/2))], [a > 0])



Now consider the second singular point t=0 in the original coordinate system.

ser := convert(`~`[series](E(a, t), t = 0, 6), polynom)

[-S(a, 0)+(9/2)*S(a, 0)*t^4+(18/5)*t^5, -3*S(a, 0)*t^2-2*t^3]


The center is the point [-S(a, 0), 0], which corresponds to the point [0, 0] on the cubical parabola. The leading terms t^4 and t^2 indicate that the tangent line is vertical and that this is a return point, with both branches lying in the second quadrant in the coordinate system shifted by [-S(a, 0), 0].

In terms of x and y, `~`[x]*alpha*y^2, and to determine the next term, we need to analyze the two series expansions. t^4 gives a t^4 term and fixes alpha. t^2*(t^3) gives a t^5 term, but the coefficient doesn't match the coefficient at t^5 in x. So the next term has to be y^(5/2) in order to match t^5. For t>0, (b*t^3+a*t^2)^(5/2) = t^5*(b*t+a)^(5/2), which gives a^(5/2)*t^5 plus higher order terms. No other terms contribute, and we can equate the coefficients at t^4 and t^5 in x and y:

collect(alpha*ser[2]^2+beta*coeff(ser[2], t, 2)^(5/2)*t^5, t); (`@`(simplify, solve))(`~`[coeff](%-ser[1], t, [4, 5]), {alpha, beta})

4*alpha*t^6+(beta*(-3*S(a, 0))^(5/2)+12*alpha*S(a, 0))*t^5+9*alpha*S(a, 0)^2*t^4


{alpha = (1/2)/S(a, 0), beta = -(4/45)*3^(1/2)/(-S(a, 0))^(5/2)}


We have obtained the coefficients in `~`[x]*alpha*y^2+beta*y^(5/2). alpha and beta are negative, while for t<0, beta has the opposite sign, and the branch corresponding to t<0 is the one which is closer to the vertical tangent.

This could have been done by using Groebner:-Basis to eliminate t from the two series and then using algcurves:-puiseux. But it is still necessary to determine the truncation order.

There is also a separate case a=0.

ser := `~`[series](E(0, t), t = 0, 6)

[series((18/5)*t^5+O(t^6),t,6), series(-2*t^3+O(t^7),t,7)]


Now the singularity is an inflection point, and the involute has a vertical tangent and moves from the second to the fourth quadrant. In the fourth quadrant, `~`[x]*alpha*(-y)^(5/3), and we find alpha in the same way as before:

coeff(ser[1], t, 5)/(-coeff(ser[2], t, 3))^(5/3)



Finally let's consider the asymptotic behavior of an involute at infinity. Now we do have to investigate S(a, t). One possible way is to apply the binomial formula to s(t) = 3*t^2*sqrt(1+1/(9*t^4))  and integrate term by term:

`assuming`([int(3*tau^2*binomial(1/2, k)*(9*tau^4)^(-k), tau = a .. t)], [0 < a and a < t])

3*binomial(1/2, k)*(9^(-k)*a^(-4*k+3)-9^(-k)*t^(-4*k+3))/(4*k-3)


This is the power series expansion of S(a, t) for large t, valid for 9*a^4 > 1. Substituting it into E(a, t), we obtain the required asymptotics.

E(a, t)

[-S(a, t)/(9*t^4+1)^(1/2)+t, -3*t^2*S(a, t)/(9*t^4+1)^(1/2)+t^3]


In the x component, we get -(1/3)*t+o(1)+t. In the y component, we get -t^3-C+o(1)+t^3. Thus the involute approaches a horizontal asymptote when t goes to +infinity, and the constant term -C gives the position of the asymptote:

`C__+` := `assuming`([unapply(sum(3*binomial(1/2, k)*9^(-k)*a^(-4*k+3)/(4*k-3), k = 0 .. infinity), a)], [a > 1/sqrt(3)])

proc (a) options operator, arrow; -a^3*hypergeom([-3/4, -1/2], [1/4], -(1/9)/a^4) end proc


That in fact extends to all positive a. For the asymptotics at -infinity, first we compute the integral over [a, -a], this time using the power series for small t. For the integral from -a to -t, S(-a, -t) = -S(a, t), and we need to subtract the value found above:

`C__-` := `assuming`([unapply(sum(int(binomial(1/2, k)*(9*tau^4)^k, tau = a .. -a), k = 0 .. infinity)-`C__+`(a), a)], [0 < a and a < 1/sqrt(3)])

proc (a) options operator, arrow; -2*a*hypergeom([-1/2, 1/4], [5/4], -9*a^4)+a^3*hypergeom([-3/4, -1/2], [1/4], -(1/9)/a^4) end proc


For negative a, `C__+` and `C__-` are swapped.

This could have been done by converting S(a, t) into a difference of two integrals of hypergeometric type, for which Maple is able to find the asymptotics automatically.

inv := proc (a, T) plots:-display(plot([op(r(t)), t = -1 .. 2], color = black), plot([op(E(a, t)), t = -1.5 .. T], color = red), plottools:-line(r(T), E(a, T), linestyle = dot), map(proc (y) options operator, arrow; plottools:-line([-1, y], [3, y], linestyle = spacedot) end proc, [-`C__+`(a), -`C__-`(a)]), scaling = constrained, view = [-1 .. 3, -1 .. 3]) end proc

Explore(inv(a, t), a = -.3 .. 1.3, t = -1. .. 4., initialvalues = [a = 1., t = 3.])



Please Wait...