iconEuler Examples

Curve Length and Curvature

This notebook computes the curve length using Maxima and Euler. For an example, we start with the logarithmic spiral.

>a=0.1;  ...

Curve Length and Curvature

>function fx(t) &= exp(a*t)*cos(t)
                              a t
                             E    cos(t)

>function fy(t) &= exp(a*t)*sin(t)
                              a t
                             E    sin(t)

Compute the elements of curve length ds, and simplify the result.

>&sqrt(diff(fx(t),t)^2+diff(fy(t),t)^2), df &= trigreduce(radcan(%))
                a t           a t        2
       sqrt((a E    sin(t) + E    cos(t))
                                           a t           a t        2
                                     + (a E    cos(t) - E    sin(t)) )

                                2       a t
                          sqrt(a  + 1) E

                                2       a t
                          sqrt(a  + 1) E

And integrate to find the curve length.

>res &= integrate(df,t,0,2*%pi)
                                     2 pi a
                            2       E         1
                      sqrt(a  + 1) (------- - -)
                                       a      a

Now we can compute the length of our spiral.


The next example is the parabola.


Curve Length and Curvature

>&integrate(sqrt(1+diff(x^2,x)^2),x,-1,1), %()
                         asinh(2) + 2 sqrt(5)


Let us compute the length of a line segment with points on the curve.

>x=-1:0.2:1; y=x^2; plot2d(x,y,r=1);  ...

Curve Length and Curvature

The result is approximately the same as the continuous result.

>i=1:cols(x)-1; sum(sqrt((x[i+1]-x[i])^2+(y[i+1]-y[i])^2))

Cartesian Curve

The next example is called the Cartesian curve. It is an implicit function.

>eq &= x^3+y^3-3*x*y
                            3            3
                           y  - 3 x y + x

First, we plot it.


Curve Length and Curvature

We are interested in the part in the positive quadrant.


Curve Length and Curvature

We solve it for x.

>&eq with y=l*x, sol &= solve(%,x)
                          3  3    3        2
                         l  x  + x  - 3 l x

                               3 l
                         [x = ------, x = 0]
                              l  + 1

And define a function in Maxima.

>function f(l) &= rhs(sol[1])
                                 3 l
                                l  + 1

Now we use the function to plot the curve in Euler.


Curve Length and Curvature

Compute the curve element ds.

                    8       6       5       3       2
            sqrt(9 l  + 36 l  - 36 l  - 36 l  + 36 l  + 9)
                        12      9      6      3
                  sqrt(l   + 4 l  + 6 l  + 4 l  + 1)

Maxima cannot integrate this.

         /          8       6       5       3       2
         [  sqrt(9 l  + 36 l  - 36 l  - 36 l  + 36 l  + 9)
         I  ---------------------------------------------- dl
         ]              12      9      6      3
         /        sqrt(l   + 4 l  + 6 l  + 4 l  + 1)

So we use fast numerical integration in Euler. To get the part in the first quadrant, we integrate from 0 to 1, and double the result.


We can put the computation of the curve length into a function in Euler. The function will use Maxima each time it is called to compute the derivative.

>function curvelength (fx,fy,a,b) ...
 return romberg(ds,a,b);

Finally the Archimedean spiral.


Curve Length and Curvature

The following function computes the line element for given expressions fx and fy in Maxima only.

>function ds(fx,fy) &&= sqrt(diff(fx,x)^2+diff(fy,x)^2)
                           2              2
                  sqrt(diff (fy, x) + diff (fx, x))

>sol &= ds(x*cos(x),x*sin(x))
                                  2                      2
          sqrt((cos(x) - x sin(x))  + (sin(x) + x cos(x)) )

>&sol | trigreduce | expand, &integrate(%,x,0,2*pi), %()
                             sqrt(x  + 1)

                  asinh(2 pi) + 2 pi sqrt(4 pi  + 1)


Another example.

>sol &= radcan(ds(3*x^2-1,3*x^3-1))
                          3 x sqrt(9 x  + 4)


Curve Length and Curvature

                                2     3/2
                            (9 x  + 4)

The Cycloid

We compute the length of the curve defined by a point on a rolling wheel. The wheel is rolling to the left, so that its center is in

Curve Length and Curvature

at time t. The wheel has radius 1 for simplicity, so that it has turned to the left once if t=2*pi. The point, where the wheel touches the ground at time t=0 follows the curve

Curve Length and Curvature

To check this, set t=0. Then observe that (sin(t),cos(t)) does indeed turn to the left.

>ex &= x-sin(x); ey &= 1-cos(x);

We plot the curve and the state of the wheel at t=2.

>plot2d(ex,ey,xmin=0,xmax=2pi,square=1); ...
   plot2d("2+cos(x)","1+sin(x)",xmin=0,xmax=2pi,>add,color=red); ...
   plot2d([2,ex(2)],[1,ey(2)],color=red,>add); ...

Curve Length and Curvature

>&sqrt(diff(x-sin(x),x)^2+diff(1-cos(x),x)^2) | radcan, f &= trigsimp(%)
                        2         2
                sqrt(sin (x) + cos (x) - 2 cos(x) + 1)

                          sqrt(2 - 2 cos(x))


Of course, numerical integration works too.



Now 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.

Curve Length and Curvature

It is the norm of the derivative vector, or speaking in terms 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 length 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.

Curve Length and 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)


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);  ...

Now, we setup a function for the curvature of a curve [fx,fy]. We can only motivate the formula. The vector

Curve Length and Curvature

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

Curve Length and 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); ...
                                  3  2
                                om  r

The speed is the following.

>&assume(om>0); &trigsimp(dsf(gx,gy,t))
                                 om r

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

Curve Length and Curvature

We take this as a definition of curvature for general curves, since we can assume the acceleration and the speed to be asymptotically constant locally.

>function kr(fx,fy,t) &&=  dsf(fx,fy,t)^3 ...
                               dsf (fx, fy, t)
        mabs(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))

We get the well known formula for the curvature of the parabola.

                                2     3/2
                            (4 t  + 1)

The general formula for the curvature of the graph of a function

Curve Length and Curvature

it the following.

>&remvalue(f); &depends(f,t); $kr(t,f,t)

Curve Length and Curvature

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

The following factor is common to both coordinates.

>function kra(fx,fy,t) &&=  dsf(fx,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.

                                [0, 0]

For the parabola, we get a special curve.

                                3  6 x  + 1
                          [- 4 x , --------]

Let us plot this in Euler. First the parabola.


Then add the curve of the midpoints.


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); ...

Curve Length and 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)", ...

Curve Length and Curvature

To see everything, we need a less complicated plot.

>a=1; k=0.15;

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)", ...

Curve Length and Curvature

What about the Archimedean 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 the curve converges to a circle.

>plot2d(mxm("ux(x*cos(a*x),x*sin(a*x),x)"), ...

Curve Length and Curvature