iconEuler Examples

Curvature

by R. Grothmann

In this notebook, we use Maxima to compute the curvature, and the
evolute of curves.

First, we define a function for the curve segment ds of a curve f(t) =
[fx(t),fy(t)] in Maxima.

Curvature

It is the norm of the derivative vector, or speaking in tems of
physics, the speed at the time t.
>function dsf(fx,fy,t) &&= sqrt(diff(fx,t)^2+diff(fy,t)^2)
                           2              2
                  sqrt(diff (fy, t) + diff (fx, t))

Of course, this is the legngth of the circumference for a circle, if
we take one turn in one second.
>&assume(r>0); &trigsimp(dsf(r*cos(2*pi*x),r*sin(2*pi*x),x))
                                2 pi r

The length of the curve is the integral over the speed. This is
obvious in physics. With the right definition for the curve length, it
can be proved exactly in math.

Curvature

Let us compute the length of the parabola from -1 to 1.
>&integrate(dsf(x,x^2,x),x,-1,1), %()
                         asinh(2) + 2 sqrt(5)
                         --------------------
                                  2

2.95788571509
We can check this numerically in Euler. We add the lengths of curve
segments for this.
>t=-1:0.01:1; s=t^2; n=cols(t);  ...
>  sum(sqrt((differences(t))^2+(differences(s))^2))
2.95787080795
Now, we setup a function for the curvature of a curve [fx,fy]. We can
only motivate the formula. The vector

Curvature

is perpendicular to f'(t) with the same length. Thus the scalar
product

Curvature

is the part of the acceleration perpendicular to the path times the
speed. Let us compute this for a general path along a circle with
constant speed.
>gx &= r*cos(om*t); gy &= r*sin(om*t); ...
>   &trigsimp(diff(gx,t)*diff(gy,t,2)-diff(gy,t)*diff(gx,t,2))
                                  3  2
                                om  r

The speed is the following.
>&trigsimp(dsf(gx,gy,t))
                              abs(om) r

If we divide the speed to the third power by the expression above, we
get the radius for this specific case.

Curvature

We take this as a definition of curvature for general curves, since we
can assume the accelaration and the speed to be asymptotically
constant locally.
>function kr(fx,fy,t) &&=  dsf(fx,fy,t)^3 ...
>  /abs(diff(fx,t)*diff(fy,t,2)-diff(fx,t,2)*diff(fy,t))
                                 3
                              dsf (fx, fy, t)
        ------------------------------------------------------------
        abs(diff(fx, t) diff(fy, t, 2) - diff(fx, t, 2) diff(fy, t))

Let us check with a circle of radius r.
>&assume(r>0); &trigsimp(kr(r*cos(om*x),r*sin(om*x),x))
                                  r

Heuer us the well know formula for the curvature of the parabola.
>&ratsimp(kr(t,t^2,t))
                                2     3/2
                            (4 t  + 1)
                            -------------
                                  2

The general formula for the graph of a function

Curvature

it the following.
>&depends(f,t); $kr(t,f,t)

Curvature

Now we want to cumpute the center of the curvature at any point
[fx,fy] depending on t. Again, we cannot privide an explanation for
the formula here.

The following factor is common to both coordinates.
>function kra(fx,fy,t) &&=  dsf(fx,fy,t)^2 ...
>  /(diff(fx,t)*diff(fy,t,2)-diff(fx,t,2)*diff(fy,t))
                               2
                            dsf (fx, fy, t)
        -------------------------------------------------------
        diff(fx, t) diff(fy, t, 2) - diff(fx, t, 2) diff(fy, t)

Now the x and y coordinates depending on x.
>function ux(fx,fy,t) &&= fx-diff(fy,t)*kra(fx,fy,t)
                   fx - kra(fx, fy, t) diff(fy, t)

>function uy(fx,fy,t) &&= fy+diff(fx,t)*kra(fx,fy,t)
                   fy + kra(fx, fy, t) diff(fx, t)

Combine to a curve.
>function krm(fx,fy,t) &&= [ux(fx,fy,t),uy(fx,fy,t)]
                    [ux(fx, fy, t), uy(fx, fy, t)]

Let us test with the circle of radius r.
>&trigsimp(krm(r*cos(x),r*sin(x),x))
                                [0, 0]

For the parabola, we get a special curve.
>&ratsimp(krm(x,x^2,x))
                                      2
                                3  6 x  + 1
                          [- 4 x , --------]
                                      2

Let us plot this in Euler. First the parabola.
>plot2d("x","x^2",xmin=-1,xmax=1,a=-1.5,b=1.5,c=-1,d=2);
Then add the midpoint curve.
>plot2d(mxm("ux(x,x^2,x)"),mxm("uy(x,x^2,x)"),xmin=-1,xmax=1,color=2,add=1);
Now compute the radius and the center at a specific point.
>xt=0; m=mxmeval("krm(x,x^2,x)",x=xt); r=mxmeval("kr(x,x^2,x)",x=xt);
>plot2d("m[1]+cos(x)*r","m[2]+sin(x)*r",xmin=0,xmax=2*pi,add=1,color=2);
>plot2d(m[1],m[2],points=1,style="o",add=1); insimg;

Curvature

Next, let us investigate the logarithmic spiral.
>a=1; k=0.15; plot2d("a*exp(k*x)*cos(x)","a*exp(k*x)*sin(x)", ...
>  xmin=-4*pi,xmax=4*pi,r=7,n=500); insimg;

Curvature

To see everything, we need a less complicated plot.
>a=1; k=0.15;
>plot2d("a*exp(k*x)*cos(x)","a*exp(k*x)*sin(x)",xmin=0,xmax=2*pi,r=3,n=500);
Let us compute the centers of curvature.
>&assume(a>0); &ratsimp(krm(a*exp(k*t)*cos(t),a*exp(k*t)*sin(t),t))
                         k t              k t
                 [- a k E    sin(t), a k E    cos(t)]

This is again a logarithmic spiral. We add it to the plot.
>plot2d("-a*k*exp(k*x)*sin(x)","a*k*exp(k*x)*cos(x)", ...
>  color=2,xmin=0,xmax=2*pi,add=1,n=500); insimg;

Curvature

What about the Archimedian spiral?
>a=2; plot2d("x*cos(a*x)","x*sin(a*x)",xmin=0,xmax=2*pi,r=8,n=500);
This spiral adds the same amount to the radius at each winding.
Let us compute the centers of curvature.
>&assume(a>0); &ratsimp(krm(x*cos(a*x),x*sin(a*x),x))
                          2  2
         a x cos(a x) - (a  x  + 1) sin(a x)
        [-----------------------------------, 
                      3  2
                     a  x  + 2 a
                                                   2  2
                                  a x sin(a x) + (a  x  + 1) cos(a x)
                                  -----------------------------------]
                                               3  2
                                              a  x  + 2 a

And add to the plot. We see that their curve converges to a circle.
>plot2d(mxm("ux(x*cos(a*x),x*sin(a*x),x)"), ...
>  mxm("uy(x*cos(a*x),x*sin(a*x),x)"),xmin=0,xmax=2*pi,add=1,color=2); insimg;

Curvature

>

Examples Homepage