#  Arclength 
# 
# Thm.  If  f' is continuous on [a,b], then the graph
#  of f over the interval has length  
>  len = Int(sqrt(1+(diff(f(x),x))^2),x=a..b);

                               b
                               /
                              |  /    /  d      \2\1/2
                       len =  |  |1 + |---- f(x)| |    dx
                              |  \    \ dx      / /
                             /
                             a
# To see this, write a sum which approximates the length: 
> len, `approx=`,Sum(sqrt((x[i+1]-x[i])^2+(f(x[i+1])-f(x[i]))^2),i=1..n);

                     n
                   -----
                    \                      2                          2 1/2
     len, approx=,   )   ((x[i + 1] - x[i])  + (f(x[i + 1]) - f(x[i])) )
                    /
                   -----
                   i = 1
# Now rewrite this approximating sum in the form
>  Sum(sqrt((x[i+1]-x[i])^2+(f(x[i+1])-f(x[i]))^2),i=1..n)=Sum(sqrt((x[i+1]-x[i])^2+ `f'`(t[i])^2* (x[i+1]-x[i])^2),i=1..n);

          n
        -----
         \                      2                          2 1/2
          )   ((x[i + 1] - x[i])  + (f(x[i + 1]) - f(x[i])) )    =
         /
        -----
        i = 1

              n
            -----
             \                      2           2                  2 1/2
              )   ((x[i + 1] - x[i])  + f'(t[i])  (x[i + 1] - x[i]) )
             /
            -----
            i = 1
# where we use the mean value theorem on each little 
# subinterval [x[i],x[i+1]].  Write delta x[i] = x[i+1]-x[i]
# and factor it out from under the square root
> Sum(sqrt((x[i+1]-x[i])^2+ `f'`(t[i])^2* (x[i+1]-x[i])^2),i=1..n) = Sum(sqrt(1+`f'`(t[i])^2),i=1..n)* Delta*x[i];

           n
         -----
          \                      2           2                  2 1/2
           )   ((x[i + 1] - x[i])  + f'(t[i])  (x[i + 1] - x[i]) )    =
          /
         -----
         i = 1

             /  n                     \
             |-----                   |
             | \                 2 1/2|
             |  )   (1 + f'(t[i]) )   | Delta x[i]
             | /                      |
             |-----                   |
             \i = 1                   /
# Now this is a Riemann sum which converges to the integral
# as the mesh of the partition goes to 0.
# 
# Arclength integrals are often impossible (or at least
# very difficult) to do symbolically.  For example,
# the length of a parabolic arc, say y = x^2 from x=0
# to x=1 is fairly difficult to do by hand, although
# Maple makes short work of it.
# 
> L := Int(sqrt(1+diff(x^2,x)^2),x=0..1) ;value(L); evalf(");

                                   1
                                   /
                                  |          2 1/2
                            L :=  |  (1 + 4 x )    dx
                                  |
                                 /
                                 0

                                 1/2
                            1/2 5    + 1/4 arcsinh(2)

                                   1.478942858
# Note how close this number is to the straight line 
# distance from 0,0 to 1,1.
# 
# If we try to calculate the length of one arch of the
# sin function,  Maple returns the question. 'Evalfing'
# the integral says use a numerical method to compute. 
> int(sqrt(1+cos(x)^2),x=0..Pi); evalf(");

                              Pi
                              /
                             |              2 1/2
                             |   (1 + cos(x) )    dx
                             |
                            /
                            0

                                   3.820197789
# Problem.   Walk 10 units along the parabola y = x^2
# starting from x = 0.  Where are you at?
# 
# Solution:  Let x1 be the x-coordinate of your ending 
# point. First estimate x1.  Since we are walking up in
# addition to the right, x1 should be less than 4.  A 
# lower bound is more difficult to achieve.
# 
# We want to solve the equation    
> eqn := Int(sqrt(1+4*x^2),x=0..x1) - 10;

                                 x1
                                 /
                                |           2 1/2
                        eqn :=  |   (1 + 4 x )    dx - 10
                                |
                               /
                               0
# for x1.  Note that the unknown is a limit of integration.
>  eqn := value(eqn);

                                  2 1/2
              eqn := 1/2 (1 + 4 x1 )    x1 + 1/4 arcsinh(2 x1) - 10
#  solve gives no information.
> solve(eqn,x1);\

# So, we can use plot to get an approximate answer
> plot(eqn,x1=0..10);
# An inspection of the graph shows that x1 is between 2 and 
# 4.  The 'exact' answer is given by fsolve.
> fsolve(eqn,x1,0..10);

                                   3.041301697
#