Curvature of 2 Dimensional Curves

The best means of visualizing the curvature of a curve in 2 dimensional space is to animate the osculating circle as the parameter increases. Since the curvature is the reciprocal of the radius of the osculating circle, the curvature increases when the osculating circle gets smaller and the curvature decreases when the osculating circle gets larger. It is this idea that we want to visualize.

First, we define a vector-valued function r , and then we compute its velocity v and its acceleration a .

> r:=[4*cos(t),3*sin(t)]; #Enter a vector function of t
v:=diff(r,t);
a:=diff(v,t);

Now let's compute the speed and the unit tangent vector. We use the simplify and expand commands to make the expressions more tractable.

> speed_raw:=sqrt(innerprod(v,v)); #speed computed but not simplified
speed:=simplify(speed_raw);
T:=simplify(expand((v/speed)));

Next we find the curvature and the unit normal. We begin by differentiating the unit tangent with respect to time.

> dT_over_dt:=simplify(diff(T,t));

The magnitude of the derivative of the unit tangent is the rate of change of the angle of inclination of the velocity vector. That is,

Subsequently, the curvature , which is called kappa , is the angular velocity over the speed, and the normal N is the unit vector in the direction of the derivative of T.

> dtheta_over_dt:=sqrt(innerprod(dT_over_dt,dT_over_dt));
kappa:=simplify(dtheta_over_dt/speed);
N:=simplify(expand(dT_over_dt/dtheta_over_dt));

We are now ready to find the center and radius of the osculating circle. In particular, the radius of the osculating circle is 1/ and the center of the osculating circle is a distance of 1/kappa from the curve along the unit normal.

Thus, the center of the osculating circle is the sum of the position vector r and the normal vector scaled by 1/kappa.

> osc_center:=simplify(expand(r+N/kappa));

Let's draw the unit tangent, the unit normal, and the osculating circle at a given point on the curve. The curve is graphed on the range

t_min .. t_max

and the point on the curve is determined by the value of t_now.

> t_min:=0;
t_max:=2*Pi;
t_now:=Pi/6;

# plot the curve, the unit tangent, and the unit normal at this point
p1:=plot([r[1],r[2],t=t_min..t_max],color=black,thickness=2):
r_point:=evalf(subs(t=t_now,r)):
T_point:=evalf(subs(t=t_now,T)):
N_point:=evalf(subs(t=t_now,N)):
p2:=arrow(vector(r_point),vector(T_point),0.1,0.2,0.1,color=red):
p3:=arrow(vector(r_point),vector(N_point),0.1,0.2,0.1,color=green):
# compute the osculating circle center and radius
kappa_now:=evalf(subs(t=t_now,kappa)):
center_now:=evalf(subs(t=t_now,osc_center)):
p4:=circle(center_now,1/kappa_now,color=blue):
display(p1,p2,p3,p4,scaling=constrained);

In order to more clearly visualize the relationship between the osculating circle and the curve, let's redraw the osculating circle continuously as we traverse the curve. We do so by creating a series of plots like the one above which are then shown in sequence.

> t_min:=0:
t_max:=2*Pi:
num_frames:=10: #Number of frames in the animation

t_delt:=(t_max-t_min)/num_frames:
plot_array:=array(1..num_frames):

#First, we draw a static image of the curve
static_plot:=plot([r[1],r[2],t=t_min..t_max],color=grey):

#Next, we create the sequence of images making up the animation
for i from 1 to num_frames do
t_now:=t_min+i*t_delt:

# plot the curve, the unit tangent, and the unit normal at this point
p1:=plot([r[1],r[2],t=t_min..t_now],color=black):
r_point:=evalf(subs(t=t_now,r)):
T_point:=evalf(subs(t=t_now,T)):
N_point:=evalf(subs(t=t_now,N)):
p2:=arrow(vector(r_point),vector(T_point),0.1,0.2,0.1,color=red):
p3:=arrow(vector(r_point),vector(N_point),0.1,0.2,0.1,color=green):

# compute the osculating circle center and radius
kappa_now:=evalf(subs(t=t_now,kappa)):
center_now:=evalf(subs(t=t_now,osc_center)):

# plot osculating circle if possible
if (kappa_now = 0) then
plot_array[i]:=display(p1,p2,p3):
else
p4:=circle(center_now,1/kappa_now,color=blue):
plot_array[i]:=display(p1,p2,p3,p4):
end if:
end do:
plot_list:=convert(plot_array,'list'):
ani_plot:=display(plot_list,insequence=true):
display(static_plot,ani_plot,scaling=constrained);