iconEuler Examples

Computing the Curve Length

This notebook computes the curve length using Maxima and Euler. For
an example, we start with the logarithmic spiral.
>a=0.1; plot2d("exp(a*x)*cos(x)","exp(a*x)*sin(x)",line=1,r=2,xmin=0,xmax=2*pi);  ...
>  insimg;

Curve Length

>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

>&df
                                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.
>res(a=0.1)
8.78817491636
The next example is the parabola.
>plot2d("x^2",xmin=-1,xmax=1,r=1); insimg;

Curve Length

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

2.95788571509
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);  ...
>  plot2d(x,y,points=1,style="[]",add=1); insimg;

Curve Length

The result is approximately the same as the continous result.
>i=1:cols(x)-1; sum(sqrt((x[i+1]-x[i])^2+(y[i+1]-y[i])^2))
2.95191957027
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.
>plot2d(eq,r=2,niveau=0); insimg;

Curve Length

We are interested in the part in the positive quadrant.
>plot2d(eq,a=0,b=2,c=0,d=2,niveau=0); insimg;

Curve Length

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]
                               3
                              l  + 1

And define a function in Maxima.
>function f(l) &= rhs(sol[1])
                                 3 l
                                ------
                                 3
                                l  + 1

Now we use the function to plot the curve in Euler.
>plot2d(&f(x),&x*f(x),line=1,xmin=-0.5,xmax=2,a=0,b=2,c=0,d=2,r=2); insimg;

Curve Length

Compute the curve element ds.
>&ratsimp(sqrt(diff(f(l),l)^2+diff(l*f(l),l)^2))
                    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.
>&integrate(%,l,0,1)
          1
         /          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)
          0

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.
>2*romberg(&(sqrt(diff(f(x),x)^2+diff(x*f(x),x)^2)),0,1)
4.91748872168
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) ...
ds=mxm("sqrt(diff(@fx,x)^2+diff(@fy,x)^2)");
return romberg(ds,a,b);
endfunction
>curvelength("x","x^2",-1,1)
2.95788571509
>2*curvelength(mxm("f(x)"),mxm("x*f(x)"),0,1)
4.91748872168
Finally the archimedian spiral.
>curvelength("x*cos(x)","x*sin(x)",0,2*pi)
21.2562941482
>plot2d("x*cos(x)","x*sin(x)",xmin=0,xmax=2*pi,square=1); insimg;

Curve Length

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), %()
                                   2
                             sqrt(x  + 1)


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

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

>plot2d("3*x^2-1","3*x^3-1",xmin=-1/sqrt(3),xmax=1/sqrt(3),square=1); insimg;

Curve Length

>&integrate(sol,x)
                                2     3/2
                            (9 x  + 4)
                            -------------
                                  9

We compute the length of the curve defined by a point on a rolling
wheel.
>ex &= x-sin(x); ey &= 1-cos(x);
First we plot the curve.
>plot2d(ex,ey,xmin=0,xmax=2pi,square=1); insimg;

Curve Length

>&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))

>&integrate(f,x,0,2*%pi)
                                  8

Of course, numerical intergration works too.
>romberg(mxm("f"),0,2*pi)
8

Examples Homepage