Plane Curves, Curvature, and Arclength
copyright © 2000 by Paul Green and Jonathan Rosenberg
In this notebook, we will use MATLAB to perform computations and obtain plots involving parametrized curves in two dimensions. We will illustrate the possibilites by considering the cycloid, which you have already seen in the text of Ellis and Gulick. It is parametrized by (t – sin(t), 1 – cos(t)). Let us enter this information into MATLAB. We enter the cycloid as a three dimensional curve in order to be able to use the cross-product later on.
[ t-sin(t), 1-cos(t), 0]
Let us compute the velocity and acceleration (i.e., the first two derivatives), the speed and the curvature for the cycloid. We will differentiate symbolically with respect to t; however, we do not need to specify t as an argument to diff because MATLAB will automatically differentiate with respect to the variable nearest to x alphabetically, and there is only one variable in this case. First we define a function veclength for computing Euclidean lengths of vectors. We don't want to call it length or norm since those names are already reserved by MATLAB for other purposes. (If you prefer, you can also use an m-file version of this function, available here.)
veclength=inline('sqrt(v*transpose(v))')
veclength =
Inline function:
veclength(v) =
sqrt(v*transpose(v))
[ 1-cos(t), sin(t), 0]
[ sin(t), cos(t), 0]
((1-cos(t))^2+sin(t)^2)^(1/2)
cyccurve=veclength(cross(cycvel,cycacc))/cycspeed^3
(((1-cos(t))*cos(t)-sin(t)^2)^2)^(1/2)/((1-cos(t))^2+sin(t)^2)^(3/2)
We can use the symbolic pretty command to improve the appearance and intelligibility of the formula for the curvature.
2 2 1/2
(((1 -
cos(t)) cos(t) - sin(t) ) )
-------------------------------------
2 2 3/2
((1 - cos(t)) + sin(t) )
Note that the numerator contains a square followed by a square root, in other words an absolute value. We can get things into better form using simplify:
-1/2*csgn(-1+cos(t))/(2-2*cos(t))^(1/2)
The csgn function here is native not to MATLAB but to Maple, and is basically the same thing as the sign function, except that it's defined for complex expressions also (hence the `c' in front of the name). Since the cosine is never bigger than 1, -1+cos(t) is never positive, so the csgn term is just –1 and:
cyccurve=(1/2)/(2-2*cos(t))^(1/2)
1/2/(2-2*cos(t))^(1/2)
Problem 1:
Find symbolic expressions for the velocity, speed, acceleration and curvature for the hypocycloid parametrized by
hypocyc=[2*cos(t)+cos(2*t),2*sin(t)-sin(2*t),0]
[ 2*cos(t)+cos(2*t), 2*sin(t)-sin(2*t), 0]
We continue our study of the cycloid. Let us plot it. We use ezplot, which will plot a parametrized plane curve provided the first two arguments are the coordinate functions. The last argument gives the range of the parameter for the plot. The command axis normal changes ezplot's default scaling, which would have set the scales on the two axes to be equal.
ezplot(cycloid(1),cycloid(2),[0,4*pi]);axis normal
We have shown two arches of the cycloid in order to show the cusp at (2p, 0). There are two natural questions. One is the length of each arch, and the other is the behavior of the curvature at the cusp. The length of the cycloid can be computed symbolically, by integrating the speed.
cyclength=int(cycspeed,0,2*pi)
4*4^(1/2)
Of course, 4^(1/2) is really 2, so this is the (exact) number 8. However, MATLAB does not simplify 4^(1/2) to 2 since it doesn't know for sure that we want the positive square root. Now let us look more carefully at the curvature
1/(2sqrt(2-2cos(t))),
which clearly becomes infinite as cos(t) approaches 1. This behaviour is not apparent from the plot, which makes it look as though the curve becomes straight near the cusp. However it does become apparent from a plot of the curvature as a function of t.
Let us see what we can learn by making a detailed plot near the cusp. We will use the cusp at the origin for the sake of simplicity.
ezplot(cycloid(1),cycloid(2),[-.1,.1]);axis normal
ezplot(cycloid(1),cycloid(2),[-.01,.01]);axis normal
ezplot(cycloid(1),cycloid(2),[-.001,.001]);axis normal
In order to interpret these plots, which appear to be quite similar,
we begin by recalling that the curvature of a plane curve can be identified
with |
|, where 
is the angle between
the unit tangent and any convenient reference ray. In this case we choose our
reference ray to be the positive y-axis. Then for small values, 
can be identified
with tan(
), which is the reciprocal (because our reference line is
vertical instead of horizontal) of the slope of the tangent line. Now 
, the change in 
across our plotting
interval, appears to be the same in all three plots, as does 
, the length of the plotted curve. However, if we look at the
scales on the axes, we observe that as the length of the plotting interval
decreases by a factor of 10, the y-scale decreases by a factor of 100 and the
x-scale decreases by a factor of 1000. Since the y-scale dominates, we can
identify 
with 
, so that 
also decreases by a
factor of 100. On the other hand, the ratio between the apparent angle with the
positive y-axis and the actual angle scales with the ratio of the x-scale to
the y-scale. This only decreases by a factor of 10. Thus as the length of the
plotting interval decreases by a factor of 10, the ratio
increases by a factor of 10, showing that the curvature
approaches 
as t approaches 0.
Problem 2: Plot the hypocycloid and compute its
length. Plot the curvature as a function of t, and make some detailed
plots near the cusp at 
. Since the hypocycloid is a closed curve, your original plot
will only need to be from 
to 
.
Except for the plots, everything we have done so far could also have been done by hand. Let us now illustrate the power of MATLAB by studying a generalization of the cycloid called the trochoid. A general trochoid (except for an overall scale factor) is parametrized by ( t – b sin(t), 1 – b cos(t)). The cycloid is the case b=1.
trochoid=[t-b*sin(t),1-b*cos(t),0]
[ t-b*sin(t), 1-b*cos(t), 0]
Let us compute the velocity and acceleration (the first two derivatives), the speed, and the curvature, for the whole family of trochoids. We can do this symbolically, by differentiating with respect to t and regarding b as a symbolic constant. Again, we do not need to specify t as an argument to diff because t is alphabetically closer to x than b is.
[ 1-b*cos(t), b*sin(t), 0]
[ b*sin(t), b*cos(t), 0]
trochspeed=veclength(trochvel)
((1-b*cos(t))^2+b^2*sin(t)^2)^(1/2)
trochcurve=veclength(cross(trochvel,trochacc))/trochspeed^3
(((1-b*cos(t))*b*cos(t)-b^2*sin(t)^2)^2)^(1/2)/((1-b*cos(t))^2+b^2*sin(t)^2)^(3/2)
trochcurve=simplify(trochcurve)
(b^2*(b-cos(t))^2)^(1/2)/(1-2*b*cos(t)+b^2)^(3/2)
2 2 1/2
(b (b - cos(t)) )
------------------------
2 3/2
(1 - 2 b cos(t) + b )
Let us now consider the trochoid with b=1/2.
[ t-1/2*sin(t), 1-1/2*cos(t), 0]
ezplot(trochhalf(1),trochhalf(2),[0,4*pi]);axis normal
trochhalfcurve=subs(trochcurve,b,.5)
1/4*4^(1/2)*((1/2-cos(t))^2)^(1/2)/(5/4-cos(t))^(3/2)
We can see in this case that the curvature does not become infinite, since the denominator cannot vanish. We can also plot the curvature as a function of t, using ezplot.
We can see both from the plot and from the formula that the curvature vanishes when cos(t)=sqrt(1/2), and is not differentiable at those points. These points are the inflection points on the trochoid.
We compute the length of this trochoid using quad8 (to see why, try computing it with double and int):
trochhalfspeed=subs(trochspeed,b,.5)
((1-1/2*cos(t))^2+1/4*sin(t)^2)^(1/2)
trochhalflength=quad8(inline(vectorize(trochhalfspeed)),0,2*pi)
6.6825
We can also see from the curvature plot that this trochoid has its maximum curvature at t=0. We can compute the curvature at this point and verify that its derivative with respect to t vanishes.
2
subs(diff(trochhalfcurve),t,0)
0
Problem 3: Find symbolic expressions for the velocity, speed, acceleration and curvature for the family of hypotrochoids defined as
hyptroch=[2*cos(t)+c*cos(2*t), 2*sin(t)-c*sin(2*t), 0]
[ 2*cos(t)+c*cos(2*t), 2*sin(t)-c*sin(2*t), 0]
Notice that the hypocycloid is the hypotrochoid corresponding to c=1. Make an analysis, similar to the one above, for the hypotrochoid corresponding to c=0.5.
Now let us consider the trochoid corresponding to b=2.
[ t-2*sin(t), 1-2*cos(t), 0]
troch2curve=subs(trochcurve,b,2)
4^(1/2)*((2-cos(t))^2)^(1/2)/(5-4*cos(t))^(3/2)
ezplot(troch2(1),troch2(2),[0,4*pi]);axis normal
We see two new features. The curvature no longer takes the value 0, corresponding to the absence of inflection points, and the trochoid intersects itself. We compute the length of this trochoid.
troch2speed=subs(trochspeed,b,2)
((1-2*cos(t))^2+4*sin(t)^2)^(1/2)
troch2length=quad8(inline(vectorize(troch2speed)),0,2*pi)
13.3649
Let us now see if we can locate the point of self-intersection on troch2. It is clear that there is an intersection point on the y-axis, and also on the line x=2n pi for every integer n. We will find the intersection on the y-axis by setting the first coordinate of troch2 to 0 and solving for t.
Inline function:
g(t) = t-2*sin(t)
Zero found in the interval: [1.7199, 4.1469].
t0 =
1.8955
0 1.6380
Problem 4:
Make an analysis of the hypotrochoid for c=2 similar to the above. Include a plot of the hypotrochoid, a plot of its curvature, a computation of the length of the hyptrochoid, and a determination of its points of self-intersection. A simultaneous plot of the hypotrochoid and the circle of radius 3 may help to determine the latter.
Although we cannot express the length of a period of the trochoid symbolically as a function of b, we can still plot it. The script m-file trochlength.m computes the length of the trochoid for a range of values of b, and plots the result.
Additional Problems:
1) Consider the plot of the length of the trochoid against the parameter b obtained above.
a) What is the exact value of the limit of the length as b approaches 0?
b) Notice that the plot appears to be asymptotic to a line for large values of b.
i) Plot one period of the trochoid for some large values of b.
ii) On the basis of your plots, guess the equation for your asymptote.
iii) Verify your guess by superimposing a plot of your conjectured asymptote on the length plot.
2) Study the maximum and minimum values of the curvature along the trochoid as functions of b.
a) Compute the first two derivatives of the curvature with respect to t.
b) Verify, using your computation in part a), that for all values of b, the curvature achieves a maximum at even multiples of pi and a minimum at odd multiples of pi.
c) Use the observation of part b) to express the minimum and maximum of the curvature as symbolic functions of b.
d) Use ezplot to plot the functions of part c) on a single diagram.
e) Relate the appearance of your plot, the symbolic form of the curvature bounds, and the trochoids you plotted in 1)b)i, to explain the behaviour of the curvature bounds for large b.
3) Observe
that any curve given in polar coordinates by 
can be parametrized
as (f(t)cos(t), f(t)sin(t)). Use this
observation to analyze each of the following curves. Plot each curve, compute
its length, plot the curvature as a function of t, determine the curvature
bounds, and identify any cusps or self-intersections.
a) The
cardioid 
.
b) The
limaçon 
.
c) The
lemniscate 
.
4) Observe
that the graph 
of a function 
can be parametrized
as (t, f(t)). Use this observation to analyze the graph of
each of the following functions. Determine the curvature bounds if they exist,
compute the length of each graph between x = 0 and x = 1, and
identify any cusps.
a) x2.
b) x3.
c) x4.
d) x3/2.
e) ex.