Cycloid

While working through the MIT ocw lectures on multi-variable calculus (Prof. Denis Auroux, here), I particularly enjoyed his discussion about the cycloid. Above is a graphic from the wikipedia article (actually the graphic is an animated gif, but I grabbed one of the frames). The red curve is generated by the motion of a point on the edge of a rolling circle.

In addition to the beauty of the curve, it turns out that the length and area under the curve have simple values that are relatively easy to calculate. See wikipedia for the details.

One thing the article doesn't explain is how to get the "parametrization" for the curve. This looks hard, but is made easy by using vectors. It's explained in the second half of Auroux's fifth lecture.

Another thing the article doesn't explain is how to integrate

√(2 - 2 cos t)

Start from the double angle formula:

cos 2s = cos2s - sin2s cos 2s = 1 - 2 sin2s 2 sin2s = (1 - cos 2s) 2 sin2(t/2) = (1 - cos t)

It's straightforward from there.

The Mathworld article is also quite nice, and references a famous challenge in history, the one which led to this quote (in reference to Newton):

"Ah, I know the lion by his paw!"

Getting a little help from Newton

I've posted here on a number of topics and derivations in fundamental calculus. This is about finding the derivative for powers of x with negative or fractional exponents. It seems a trivial proof, until I got stopped in my tracks by the cofactors of the binomial distribution.

The idea is to find the slope of the curve at x by adding a little bit to x, (x + h), calculating f(x + h), then subtract f(x), and finally divide by h and finding the limit as h gets very small. It works great for y = x². But consider

y = x1/2

We want to do:

1/h [(x + h)1/2 - x1/2]

So thinking about the binomial expansion of the x + h term, I'm trying to figure what is the expansion for the 1/2 power? I mean, what can "n choose k" mean when n is 1/2?

The other thing that is a good approach here is to bring x out of the sum as follows

(x + h)1/2 = x1/2(1 + h/x)1/2

To expand this, we just remember Newton's super-duper binomial expansion from here. If the expression is

(1 + Q)m/n

the expansion is

1 + (m/n)Q + (m/n)((m-n)/2n)Q2 + (m/n)((m-n)/2n)((m-2n)/3n)Q3 + ..

We only need the first two terms because Q² is (h/x)² which we will ignore. So we have for the expansion

(1 + h/x)1/2 = 1 + 1/2(h/x)

and the derivative is the limit as h => 0 of

1/h [x1/2(1 + 1/2(h/x)) - x1/2] = 1/h (1/2 h) 1/x1/2 = 1/2 x-1/2

Similarly for m = -2 and n = 1 the first two terms are 1 and -2 and the expansion we need will be

(1 + (h/x))-2 = = 1 - 2(h/x)

The derivative is the limit as h => 0 of

1/h [x-2 (1 - 2(h/x)) - x-2] = = 1/h x-2 (-2)(h/x) = -2 x-3

Easy when you know how.

The slope of the sine curve

Just to finish with the subject from last time (here), let's show that the slope of the sine curve at x is cos(x) and for the cosine curve at x it is -sin(x).

We start from x, and then move a little bit h. Using the rule for sum of sines (here):

sin(x + h) = sin(x) cos(h) + cos(x) sin(h) f(x + h) - f(x) = = sin(x) cos(h) + cos(x) sin(h) - sin(x) d/dx sin(x) = lim(h -> 0) [f(x + h) - f(x)]/h = 1/h [sin(x) cos(h) + cos(x) sin(h) - sin(x)] = -(1/h) sin(x) (1 - cos(h)) + (1/h) cos(x) sin(h)

The first term is -sin(x) times (1/h)(1 - cos(h)); last time we showed that (1/h)(1 - cos(h)) equals zero in the limit as h -> 0.

The second term is cos(x) times (1/h) sin(h); we showed that (1/h) sin(h) approaches 1 as h -> 0. Thus,

d/dx sin(x) = cos(x)

Since sine and cosine are periodic with cosine "ahead"

sin(x + π/2) = cos(x) cos(x + π/2) = -sin(x)

If

y = cos(x) = = sin(x + π/2)

Let

u = x + π/2

Then

d/dx cos(x) = = d/dx sin(x + π/2) = d/du sin(u) d/dx (x + π/2) = d/du sin(u) = cos(u) = cos(x + π/2) = -sin(x) d/dx cos(x) = -sin(x)

Know your limits

To find the derivative of the sine function, we need to determine two limits. The first is:

lim (h -> 0) sin(h)/h

In the figure, h is the measure of the angle and also the length of the arc of the circle between A and D (in radians). (I'm using h instead of x because ultimately this will be a small change to x, as in x + h).



Kline sets up this inequality:

area OAB < area OAD < area OED area OAB = (1/2) OB AB = (1/2) cos(h) sin(h) area OAD = h/(2π) π 12 = (1/2) h area OED = (1/2) OD ED = (1/2) sin(h)/cos(h)

The last result follows since OAB and OED are similar triangles so that ED/OD = AB/OB, and OD = 1. Thus, the inequality becomes:

(1/2) cos(h) sin(h) < (1/2) h < (1/2) sin(h)/cos(h)

As long as we do not actually reach h = 0, we are allowed to divide by sin(h):

cos(h) < h/sin(h) < 1/cos(h)

As h approaches 0, cos(h) and its inverse 1/cos(h) both approach 1 so the ratios h/sin(h) and sin(h)/h both approach 1 as well.

Strang does the first part by extending the segment AB (dotted line) and noting that

2 AB < 2 h 2 sin(h) < 2h sin(h)/h < 1

The other thing we need is:

lim (h -> 0) (1 - cos(h))/h

Following Strang:

sin2(h) + cos2(h) = 1 sin2(h) = 1 - cos2(h) = (1 + cos(h)) (1 - cos(h))

Since sin(h) < h (see the figure again):

(1 + cos(h)) (1 - cos(h)) = sin2(h) < h2

Divide by h and also by 1 + cos(h):

(1 - cos(h))/h < h/(1 + cos(h))

Strang says, note that all terms are positive, so that:

0 < (1 - cos(h))/h < h/(1 + cos(h))

Now, as as h -> 0 the right side must go to 0/2 (that is, 0) and our ratio is "caught in the middle."

Here is a simple Python script to calculate the same values. Notice that although 1-cos(h) and h are both approaching 0, h is approaching more slowly, so the ratio tends to the 0 in the limit. First the output:

h 1.0000000000, sin(h) 0.8414709848, h/sin(h) 0.8414709848 h 0.5000000000, sin(h) 0.4794255386, h/sin(h) 0.9588510772 h 0.1666666667, sin(h) 0.1658961327, h/sin(h) 0.9953767962 h 0.0416666667, sin(h) 0.0416546114, h/sin(h) 0.9997106733 h 0.0083333333, sin(h) 0.0083332369, h/sin(h) 0.9999884260 h 0.0013888889, sin(h) 0.0013888884, h/sin(h) 0.9999996785 h 0.0001984127, sin(h) 0.0001984127, h/sin(h) 0.9999999934 h 0.0000248016, sin(h) 0.0000248016, h/sin(h) 0.9999999999 h 0.0000027557, sin(h) 0.0000027557, h/sin(h) 1.0000000000 h 1.0000000000, cos(h) 0.5403023059, (1-cos(h))/h 0.4596976941 h 0.5000000000, cos(h) 0.8775825619, (1-cos(h))/h 0.2448348762 h 0.1666666667, cos(h) 0.9861432316, (1-cos(h))/h 0.0831406106 h 0.0416666667, cos(h) 0.9991320700, (1-cos(h))/h 0.0208303194 h 0.0083333333, cos(h) 0.9999652780, (1-cos(h))/h 0.0041666426 h 0.0013888889, cos(h) 0.9999990355, (1-cos(h))/h 0.0006944443 h 0.0001984127, cos(h) 0.9999999803, (1-cos(h))/h 0.0000992063 h 0.0000248016, cos(h) 0.9999999997, (1-cos(h))/h 0.0000124008 h 0.0000027557, cos(h) 1.0000000000, (1-cos(h))/h 0.0000013779 h 0.0000002756, cos(h) 1.0000000000, (1-cos(h))/h 0.0000001378 h 0.0000000251, cos(h) 1.0000000000, (1-cos(h))/h 0.0000000133 h 0.0000000021, cos(h) 1.0000000000, (1-cos(h))/h 0.0000000000 h 0.0000000002, cos(h) 1.0000000000, (1-cos(h))/h 0.0000000000 h 0.0000000000, cos(h) 1.0000000000, (1-cos(h))/h 0.0000000000

and the listing:

from __future__ import division from math import sin, cos h = 1 for i in range(1,10): h /= i y = sin(h) print 'h %3.10f, sin(h) %3.10f, h/sin(h) %3.10f' % (h, y, y/h) print h = 1 for i in range(1,15): h /= i y = cos(h) print 'h %3.10f, cos(h) %3.10f, (1-cos(h))/h %3.10f' % (h, y, (1-y)/h)

Area under the curve

I had a post the other day about the product rule in differential calculus, and how to extend it to find the derivative of a reciprocal and a quotient. The product rule also does heavy duty in the opposite direction, in the method called integration by parts. Strang has a wonderful diagram of that (shown above), which I've reproduced here using matplotlib.

The product rule is:

(uv)' = u v' + v u'

Integrating both sides and isolating u v' we get:

∫u v' = u v - ∫v u'

where for the definite integral we have to remember to evaluate u times v at both of the limits of integration, subtracting the value at the lower limit from the upper one.

∫u v' = u2 v2 - u1 v1 - ∫v u'

Can you see each of these four terms as areas in the diagram? It's very pretty.

from math import sqrt import numpy as np import matplotlib.pyplot as plt @np.vectorize def f(x): return 100*sqrt(x) + 0.2*x**3 # the curve dx = 0.01 X = np.arange(0,10,dx) plt.plot(X,f(X),'r',lw=6,zorder=1) # the points u1 = 2 u2 = 10 plt.scatter(u1,f(u1),s=250, color='k',zorder=2) plt.scatter(u2,f(u2),s=250, color='k',zorder=2) # u dv X2 = np.arange(u1,u2+dx,dx) plt.fill_between(X2,f(X2),0, color='0.85',zorder=0) # v du YMAX = f(10) plt.fill_between(X2,YMAX,f(X2), color='0.35',zorder=0) X3 = np.arange(0,u1+dx,dx) plt.fill_between(X3,YMAX,f(u1), color='0.35',zorder=0) plt.text(6,100,s='$\int$ $u$ $dv$', fontsize=24) plt.text(2,300,s='$\int$ $v$ $du$', fontsize=24,color='w') plt.text(0,525, s='$y = 100 \sqrt{x}$ + $0.2 {x^3}$', fontsize=24,color='r') plt.text(2.5,115,s='$u_1$,$v_1$', fontsize=18) plt.text(10.5,500,s='$u_2$,$v_2$', fontsize=18) ax = plt.axes() r1 = plt.Rectangle((0,0),width=u1, height=f(u1),facecolor='cyan') ax.add_patch(r1) ax.set_xlim(-1,12.2) ax.set_ylim(-10,600) plt.savefig('example.png')

Calculus in 5 minutes (last part)

I am unable to resist saying a word or two about how we get the standard method. I can't emphasize enough that you don't need to know this to use calculus! But it might be interesting. Here goes...

Remember:

y = xn y' = n xn-1

Where does this come from?

Well, what is the slope of the curve at any point x? Imagine that we increase x by a small amount, we'll call that amount h. The slope of the curve is the difference in y [ the value of the function f(x) ], between the two points x and x + h. How much does it change?

Suppose:

y = x3 slope = [ f(x + h) - f(x) ] / h f(x + h) = (x + h)3 = ?

Well,

(a + b)n is the binomial expansion of a + b (a + b)n = an + n an-1 b + stuff times b2 + more stuff

So...

(x + h)n = xn + n xn-1 h + stuff times h2 + more stuff

And there are three quick points:

• (1) The first term xⁿ is what we need to subtract, it is f(x)

• (2) The third and following terms involve h² (and higher powers). We ignore these terms!

Suppose h = 10-6, then h2 = 10-12 Not small enough? Suppose h = 10-60, then h2 = 10-120 Not small enough?…

The point is that h² goes away faster than h as h gets very small, and if you don't think it's getting small fast enough, just choose h to be even smaller...

• (3) What remains is the first term of the binomial expansion:

n xn-1 h

To get the slope, we need to divide by h:

slope = n xn-1 h / h

The h / h part is still equal to 1 even when h gets very small. So it goes away!

slope = n xn-1

There you have it.

As I said, this is for my son. Here is a picture from an earlier time:

Calculus in 5 minutes (part 4)

If you haven't read the earlier posts, they are here, here and here.

Now that was interesting, you say. But how about a little proof ? For a mathematician, the word "proof" is like waving the red handkerchief to a bull at Pamplona, but I'm not going there. However, it is always a good idea to check things a bit to see if they make sense. Consider the simple equation:

y = x2

We want to know the area under the curve.



I know how to do that, you say. Remembering that the curve itself is the derivative of the area function f(x), we have:

y' = x2 y = 1/3 x3 evaluated between x = 1 and x = 0: area = 1/3

That was way too easy. But consider this: how about the area above the curve?



Notice that we have

y = x2 x = √y = y1/2

I get a little freaked out seeing x as a function of y, so I am going to mentally turn the curve through 90 degrees and then make its mirror image reflection. I get this:



Now, our differentiation trick works for fractional exponents. So we have:

y' = x1/2 y = 2/3 x3/2 evaluated between x = 1 and x = 0: area = 2/3

And, as we might hope and expect, 2/3 + 1/3 does equal the area of the unit square. And note that it will work for any power of x. If, for example, we start with the curve y = x³, then the first area is 1/4 x⁴, and the second area is 3/4 x^4/3. It's easy to see that it works for any rational power of x. Suppose we have:

y' = xa/b the exponent is a/b, let a/b + 1 = c/b y = b/c xc/b x' = yb/a b/a + 1 = d/a x = a/d yd/a From above: a/b + 1 = c/b a + b = c Also from above: b/a + 1 = d/a a + b = d So: c = d a + b = c a/c + b/c = 1 a/d + b/c = 1

As we guessed, we see that b/c + a/d = 1.

(Of course, it helps that we chose the upper limit for measuring the area at x = 1, y = 1. Evaluation for some other bound would involve a bit more calculation).

Ladies and gentlemen, will you please…give it up for Mr. Leibnitz and Mr. Newton.

[Update: blogger messed with my code. Think I fixed it:]

R code:

f <- function(x) { return (x**2) } plot (f,0,1,cex=3) xvals = seq(0,1,length=1000) yvals = f(xvals) x = c(xvals,rev(xvals)) y = c(rep(0,1000),rev(yvals)) polygon(x,y,col='gray') f <- function(x) { return (x**2) } plot (f,0,1,cex=3) xvals = seq(0,1,length=1000) yvals = f(xvals) x = c(xvals,rev(xvals)) y = c(rep(1,1000),rev(yvals)) polygon(x,y,col='gray') f <- function(x) { return (x**0.5) } plot (f,0,1,cex=3) xvals = seq(0,1,length=1000) yvals = f(xvals) x = c(xvals,rev(xvals)) y = c(rep(0,1000),rev(yvals)) polygon(x,y,col='gray')

Calculus in 5 minutes (part 3)

This is a series about the essence, the merest whiff of the perfume of calculus. The first and second posts are here and here. Differentiation is useful, but integration allows us to do things which are nothing short of miraculous. For example, consider this modestly complex equation, and its curve.

y = x3 - x2 - x + 150

The integral calculus (the second part) allows us to easily calculate the area under this curve (between the curve and the x-axis), as shown in gray.

Remember what we said before: the relationship between the function which we plot as a curve, and its slope at any point is found by differentiating f(x) to f '(x). The slope is f '(x). The secret of integral calculus (don't tell anybody) is that the same relationship holds between the area under a curve, and the curve itself, but in reverse. That is, if the curve is f '(x), the area is f(x). It's that simple.

If we know f '(x), how do we find f(x)? Sometimes it's easy and sometimes it's hard. Suppose that every time we have a function f(x) and we find f '(x), we write the result down and save it in our little black book. When we encounter a function in an integration problem, we look in the book.

For the example, we need a function f(x), which when differentiated gives:

y' = x3 - x2 - x + 150

We look in our book, and there it is:

y = x4/4 - x3/3 - x2/2 + 150 x

Of course, the area depends on which endpoints we choose. (Looking at the figure should make this clear). So evaluate the integrated function between the limits x = 7 and x = -5. We have:

74/4 - 73/3 - 72/2 + 7*150 minus (-5)4/4 - (-5)3/3 - (-5)2/2 + (-5)*150

I get:

2401/4 - 343/3 - 49/2 + 1050 - 625/4 + 125/3 - 25/2 - 750

I'm too lazy to finish it. But I'll check the method itself by an interesting approach in the next post.

Why do I say it's a miracle? Because we get the area from the integrated equation just by using the endpoints. We don't have to worry about what lies between. That is amazing.

R code: f <- function(x) { return (x**3 - x**2 - x + 150) } plot (f,-5,7,xlim=c(-6,8),cex=3) xvals = seq(-5,7,length=1000) yvals = f(xvals) x = c(xvals,rev(xvals)) y = c(rep(0,1000),rev(yvals)) polygon(x,y,col='gray')

Calculus in 5 minutes (part 2)

In the first post of this series I mentioned that a car's odometer and speedometer perform related functions :), and these functions are a good way to think about the fundamental operations of calculus. Suppose a car is going in a straight line (say, on the Bonneville salt flats). We need to deal only with the positive x-axis.

Distance (often labeled s) is a function of time:

s = f(t)

If the velocity is constant and we start from s = 0 then:

s = v * t

Distance equals velocity times time. (60 miles per hour for 2 hours = 120 miles).

Now, velocity might be changing, it might be a function of time. Let's say we'll figure out how it changes later, and just write in a general way:

v = f(t)

If the velocity is changing, we call the rate of change in velocity the acceleration. (It is the slope of the curve for v as a function of t). If the acceleration is constant (like gravity), then:

v = a * t

For a car, imagine increasing the pressure on the gas pedal steadily so that after 1 second you are going 10 mph, after 2 seconds 20 mph, after 3 seconds 30 mph. If we continue at the same rate, we'll accelerate from 0 to 60 mph in 6 seconds.

As I said before, there is a special trick symbolized with a prime mark placed next to the s (the physicists put a dot on top). The trick is called differentiation:

s' = ... something.

What it is depends on the exact form of s as a function of time. If the acceleration is constant:

s = 1/2 a t2

Do you remember the rule from last time?

y = c xn y' = n c xn-1

Or, using the time and distance symbols:

s = 1/2 a t2 s' = v = a t s'' = a

So that is where the factor of 1/2 comes from in the equation involving acceleration. It is needed to cancel the the 2 that comes out of the exponent when we differentiate. One way to think about this is to say that the velocity after a period of constant acceleration is the average of the initial acceleration and the final acceleration, times the time. But this differentiation mumbo jumbo will work even if the velocity and acceleration are more complicated functions of time.

Calculus in 5 minutes (part 1)

Slope

Here is a post I wrote for my son. I hope he reads it! The essential ideas of calculus are surprisingly simple. See if you don't agree.

We can talk about the slope of a curve at a single point. This seems odd because "slope" refers to a change in y divided by a change in x. We ought to need two points, not just one. But we don't.

A single point can have a slope for the same reason that the earth seems to be flat where I'm standing right now even though I know it is curved---it is flat, on the ocean or the Mohave desert anyway. The carpenters who built your house assumed that the earth's curvature is negligible, and we will too. We can mentally zoom in indefinitely (ask an ant if the earth is flat).



So for any given point on a curve, we can always find another point very close nearby which has the same slope, and then we know the slope between them is constant (within any limits you want to establish).

---

Like driving a car

You can understand the basic procedures in calculus by thinking about driving a car (Strang). Most drivers refer frequently to the speedometer (measures velocity), if they're measuring distance they refer to the odometer.

If we plot accumulated distance as a function of time, the slope of that curve is the velocity shown on the speedometer. [This analogy will become useful in part 2].

---

Simple tricks

There are tricks to use to find the slope from the equation for a curve. Where the tricks come from isn't hard to understand, but it's not necessary---just use the trick. Consider the equation:

(y-1) = (x-2)2

Multiply it out to make things a bit more interesting:

y = x2 - 4x + 5

The trick is that the slope of the curve (designated as y') is easily calculated:

y = x2 - 4x + 5 y' = 2x - 4

We used three rules together:

(i) if we have various parts being added together, we take the slope for each part separately and then add them up.

(ii) the slope of a constant term (like y = 3) is zero (obviously, a horizontal line is flat, and has slope = 0).

(iii) if (c stands for constant):

y = c xn; y' = c n xn-1

---

So what

One reason this is useful is to find maximum and minimum points of a curve. These are points at which the slope is zero. For our example

y' = 2x - 4 0 = 2x - 4 x = 2

This happens when

y = x2 - 4x + 5 y = 4 - 8 + 5 = 1

This is easy to see by going back to the first form of the equation. The minimum is when the term in parentheses on the right is zero.

(y-1) = (x-2)2

There are other simple tricks for terms like 2^x and sin x.

R code:

f <- function(x) { x**2 - 4*x + 5} plot(f,from=-3, to = 7) plot(f,from=-3, to = 7,col='red',lwd=3) lines(c(-3,7),c(1,1),lty=2) lines(c(2,2),c(0,25),lty=2)

Slope of the sine curve

In another post about representing the sine and cosine as infinite series, I tried to show how the series form easily demonstrates the correctness of the derivatives of these functions, that if

y = sin(x) dy/dx = cos(x) y = cos(x) dy/dx = -sin(x)

But unfortunately this argument is circular (I think), since the series comes from a Taylor series, which is generated using the first derivative (and the second, third and more derivatives as well).

In this post, I want to do two things. Strang has a nice picture which makes clear the relationship between position on the unit circle and speed. This is it:



The triangles for velocity and position are similar, just rotated by π / 2.

It is clear from the diagram that at any point, the y-component of the velocity is cos(t), while the y-position is sin(t). Thus, the rate of change of sin(t) is cos(t). This is the result we've been seeking. Similarly, the rate of change of the x-position, cos(t), is -sin(t).

Strang also derives this result more rigorously starting on p. 64. That derivation is a bit complicated, although not too bad, and I won't follow the whole thing here. It uses his standard approach as follows:

dy / dx = lim(h -> 0) [sin(x + h) - sin(x)] / h

Applying a result found using just the Pythagorean theorem earlier (p. 31) for sin (s + t):

sin(s + t) = sin(s) cos(t) + cos(s) sin(t) cos(s + t) = cos(s) cos(t) - sin(s) sin(t)

He comes up with this expression for:

Δy / Δx = sin x [cos(h) - 1 / h] + cos x (sin h / h)

The problem is then to determine what happens to these two expressions in the limit as h -> 0. The first one is more interesting. As h gets small, | cos(h)-1 | gets smaller like h², so the ratio goes to 0.

R can help us see better.

h = 1 for (i in 1:6) { h = h/10 print (h) print (cos(h)-1) }

[1] 0.1 [1] -0.004995835 [1] 0.01 [1] -4.999958e-05 [1] 0.001 [1] -5e-07 [1] 1e-04 [1] -5e-09 [1] 1e-05 [1] -5e-11 [1] 1e-06 [1] -5.000445e-13

Here is the plot for the second one, which converges to 1, leaving us with simply cos(x):

f<-function(x) { x } curve(sin(x),from=0,to=0.5, ylim=c(0,0.6), lwd=3,col='red') curve(f,from=0,to=0.5, lwd=3,col='blue',add=T)

Area of a circle: simple calculus

I posted some time ago about a formula to find the area of the circle (knowing its circumference). I think it's credited to Euclid, although Google isn't being my friend at the moment, so I'm not sure.

To brush up on calculus, let's explore methods to find the area (as described in Strang's Calculus. He says:

"The goal here is to take a first step away from rectangles...that is our first step toward freedom, away from rectangles to rings."

In the example on the left panel of the figure, imagine slicing up the circle in the way of any standard integration problem to find the area under a curve. Simplify by considering only the area above the dotted line (y = 0). If the radius is R, then we can write y as a function of x:

y = (R2 + x2)1/2

Immediately, we run into a problem. I don't know a function which gives this as its derivative, so I don't know how to integrate it. Let's leave it until the end.

In the second approach, we picture the circle as being built up out of rings. As the radius r varies from 0 to R, the circumference of the ring is 2πr and we need to integrate:

A = ∫ 2πr dr = πr2 + C

evaluated between r = 0 and r = R:

A = πR2

As Strang emphasizes, now we can see the geometrical reason why the the derivative of the area is the circumference! How fast does the area grow? It grows in proportion to the circumference.



Exactly the same argument explains why the surface area of a sphere (4πr²) is the derivative of the volume (4/3πr³).

The third approach is to visualize the circle being sliced like a pie. This is really the same method we used in the previous post, which employed calculus on the sly. Now we have a series of triangles. The angle at the vertex of each thin triangle is dθ. The height is just R, and the base of each triangle is proportional to R, it is R dθ. So the area of the thin triangle is:

dA = 1/2 R2 dθ

We integrate (R is a constant) and evaluate between θ = 2π and θ = 0:

A = 1/2 R2 θ A = 1/2 R2 2π A = π R2

Back to the first method. We will try numerical integration. We will compute the area of the quarter-circle between x = 0 and x = R, and then simplify even further by considering just the unit circle with R = 1. We do this in a naive way, dividing the figure into a large number of thin segments, and calculate:

y = (1 - x2)1/2

We get a better approximation if we divide the very first interval in half. In R:

y <- function(x) { sqrt(1 - x**2) } d = 0.00001 x = seq(0,1,by=d) f = y(x) S = sum(f[-1]) + 0.5*f[1] S = S*d

> S [1] 0.7853982 > S*4 [1] 3.141593

Volume of a sphere: calculus

This derivation is straight from the book Strang's Calculus. We want to calculate the volume of a hemi-sphere as shown in the figure. The cross-sections of this solid are semi-circles. The area of each semicircle is:

A = 1/2 πr2 = 1/2 π(R2-x2)

To find the volume, we add up (integrate) the areas of each little slice:

V = ∫ A(x) dx = ∫ 1/2 π(R2-x2) dx = 1/2 π [ R2x - 1/3 x3 ]

Evaluate the expression in brackets between x = -R and +R:

= [R3 - 1/3 R3] - [-R3 + 1/3 R3] = 2/3 R3 - - 2/3 R3 = 4/3 R3 V = 1/2 π 4/3 R3 = 2/3 π R3

The total volume is twice this, or

V = 4/3 π R3

The volume can be found in either cylindrical or spherical coordinates. For spherical coordinates we have:

x = r cosθ sinφ y = r sinθ sinφ z = r cosφ

That is, φ is the angle of the vector r with the z-axis, so the z-component is r cosφ. The component orthogonal to that is the projection of r on the x,y-plane, and that is r sinφ. Then, θ is the angle of that projection (and r) with respect to the x-axis. The x- component is r sinφ cosθ (the order is not important), while the y-component is r sinφ sinθ.

In this case, the integral is a triple integral:

2π π R V = ∫ ∫ ∫ r2 sinφ dr dφ dθ 0 0 0

From reading Chapter 14 on multiple integrals in Strang, it seems that what we do here is to first integrate with respect to r, holding the angles constant:

2π π V = 1/3 R3 ∫ ∫ sinφ dφ dθ 0 0

UPDATE: A sharp-eyed reader caught a silly error in the part that follows, so it's been edited. I had set up the problem in the first half (see the figure at the beginning of the post) to use the hemisphere. But in the second part we have two angles, θ and φ, with different limits of integration.

Continuing with the standard approach in two dimensions, we let θ go from 0 to 2π. Now φ will go from 0 to π. Imagine rotating the great circle in the xy-plane around the x-axis: we need only to go from 0 to π. We could do the hemisphere by having φ go from 0 to π/2, but it would be a weird volume, with two wedges joined at the x-axis.

So, we have the integral of sinφ dφ, which is - cosφ, evaluated between φ = 0 and φ = π =>

= -cos(pi) - -(cos(0)) = -(-1) + 1 = 2 2π V = 2/3 R3 ∫ dθ 0

But this is just θ evaluated between θ = 0 and θ = 2π => 2π, so finally we have:

V = 4/3 R3 π

(Figures are from the book)

UPDATE 2: Coming back a couple of years later, I notice Google has made changes to blogger (editing the html) that retroactively messed up the formatting on this post. The limits of integration aren't shown at the correct column offset any more. My apologies, but I don't see a simple fix at the moment.

Volume of a cone: simple calculus

I've been brushing up on my calculus skills, though it can be challenging for an old fart like me. (In one ear and out the other). I saw a book online (Strang's Calculus), which I like because it is succinct, and also the author has a very interesting, idiosyncratic style. I liked it so much I got a hard copy. Now that I've studied a bit, I realize that the book is more than succinct, it's dense, filled with math, and contains a lot more than just calculus. It's a great book.

In Chapter 8, Applications of the Integral, we encounter Example 11, to find the volume of a cone. I've redrawn the diagram from the book, below. The idea is to move vertically from the top to the bottom, letting x be equal to the radius at each point. At each value for x, we draw a shell (the outside surface of a cylinder) as shown in red. As we move from top to bottom, we accumulate a collection of these cylinders whose volumes are summed to get the volume of the cone.

The key is to recognize that the distance from the top of the cone to the top of each cylinder has the same relationship to x as b does to r (e.g. when x = r, then this distance = b). The smallest and largest rectangles in the figure are similar.

The height of the cylinder plus this distance equals b.



So the height of each cylinder is

h = b-bx/r

From this point it's easy. The area of the cross-section of each shell is its diameter (2*x), times &pi times the small width dx; the volume is this area times the height. I'm going to leave out the multiplication symbol I usually use (*) for clarity:

= 2πxhdx = 2πx(b-bx/r)dx

We sum (integrate) all these cylinders:

= ∫2πx(b-bx/r) dx = ∫2πxb dx - ∫2πxb(x/r) dx = πx2b - (2/3)πx3b/r

evaluated between x = 0 and x = r:

= πr2b - (2/3)πr2b = (1/3)πr2b