# Rocket Propulsion, Part 2

See Part 1, where I introduced the problem. Constant acceleration of an object (such as a spacecraft) produces a linearly increasing velocity and therefore a quadratically increasing kinetic energy (since EK = ½mv²). If it took a constant energy supply to produce the constant acceleration, there would seem to be a mismatch between the amount of energy consumed for propulsion and the amount of kinetic energy gained. Momentum must be conserved, so for the craft to increase velocity forward, it must eject some mass backward. I explored this using the case of a thruster that ejected a small amount of mass per second, and it turns out that while the kinetic energy of the craft increases quadratically, there is far more kinetic energy given to the expelled propellant, and the total kinetic energy appeared to increase linearly, at least in a spreadsheet that evaluated the first 60 seconds.

I would like to make this approach a little more rigorous. In particular, I assumed a drive that ejected a constant amount of mass each second, and looked at the behavior over the first minute. I am hoping that this approximates the continuous case, but I’ll want to verify this by using very small time intervals, if needed. And I also want to explore the behavior of the rocket over a much longer period of time. For this, I will need more than a spreadsheet — I will write a program in Python. Also, I assumed that the drive ejected the mass at a constant velocity relative to the craft. There is no guarantee that this will take the same amount of energy over time, since the mass of the craft changes. Let’s work it out more rigorously.

To make the math simpler, for each time interval, let’s choose a reference frame moving with the craft, so the initial velocity = 0. The craft will consume a certain amount of energy E to push apart the reaction mass and craft. Let’s assume that the drive is perfectly efficient and that all the energy consumed goes to increasing the kinetic energy of the objects. Let m1 be the mass of the craft (after the mass is ejected), and m2 be the reaction mass ejected. The final velocities will be v1 and v2, respectively. Momentum (p=mv) must be conserved, so the initial and final momentum must be the same:

$p_i = p_f$

$(m_1 + m_2) \cdot 0 = m_1 v_1 + m_2 v_2$

$m_1 v_1 + m_2 v_2 = 0$

Kinetic energy (EK = ½mv²) is not conserved; we are adding an amount of energy E each time interval.

$E_{Ki} + E = E_{Kf}$

$\frac{1}{2}(m_1 + m_2) \cdot 0^2 + E = \frac{1}{2}m_1 {v_1}^2 + \frac{1}{2}m_2 {v_2}^2$

$\frac{1}{2}m_1 {v_1}^2 + \frac{1}{2}m_2 {v_2}^2 = E$

I’m primarily interested in v1, the final velocity of the craft, so I’ll solve the momentum equation above for v2

$m_1 v_1 + m_2 v_2 = 0$

$m_2 v_2 = - m_1 v_1$

$v_2 = -\frac{m_1 v_1}{m_2}$

and substitute it into the kinetic energy equation:

$\frac{1}{2}m_1 {v_1}^2 + \frac{1}{2}m_2 {v_2}^2 = E$

$\frac{1}{2}m_1 {v_1}^2 + \frac{1}{2}m_2 (-\frac{m_1 v_1}{m_2})^2 = E$

$m_1 {v_1}^2 + m_2 (-\frac{m_1 v_1}{m_2})^2 = 2E$

$m_1 {v_1}^2 + \frac{m_2 {m_1}^2 {v_1}^2}{{m_2}^2} = 2E$

$m_1 {v_1}^2 + \frac{ {m_1}^2 {v_1}^2}{m_2} = 2E$

${v_1}^2 (m_1 + \frac{{m_1}^2}{m_2} ) = 2E$

${v_1}^2 (\frac{m_1 m_2}{m_2} + \frac{{m_1}^2}{m_2} ) = 2E$

${v_1}^2 (\frac{m_1 m_2 + {m_1}^2}{m_2} ) = 2E$

${v_1}^2 [\frac{m_1 (m_1+m_2)}{m_2}] = 2E$

${v_1}^2 = 2E \frac{m_2}{m_1 (m_1 + m_2) }$

$v_1 = \pm \sqrt{2E \frac{m_2}{m_1 (m_1 + m_2) }}$

I am specifying that the propellant is ejected in the negative direction (v2 < 0 ), so logically and by the momentum equation above, the craft will move in the positive direction (v1 > 0); therefore, we can discard the negative solution.

$v_1 = \sqrt{2E \frac{m_2}{m_1 (m_1 + m_2) }}$

This answer looks reasonable. Dividing energy by mass (m1), doubling, and taking the square root are just the inverses of the kinetic energy equation. The $\frac{m_2}{m_1 + m_2}$ is the fraction of the original mass that was ejected.

At last, I’m ready to code this into my program. The code will reevaluate the variables each timestep Δt. Given a predefined rate of propellant use (in mass/time), we know how much is used each interval: multiply the rate of propellant use by the timestep to get m2. We can subtract that from the ship’s prior mass to get its new mass, m1. We’ll use the velocity equation derived just above to derive its velocity increase (recall that I defined those equations in a reference frame where the craft was starting at rest), so we’ll add that to the craft’s prior velocity to get its new velocity. The energy E added each step is the predefined rate of energy production multiplied by the timestep.

Let’s use initial conditions similar to the example in the last post. I’ll start with a 1,000-kg craft. It will eject 1 g/s, and it will be powered by a 500-W drive (since in the last post, the kinetic energy increased by 5,000,005 mJ / 10 s = 5,000.005 J / 10 s = 500.0005 W). Running the simulation for 60 seconds yields a velocity of 0.0600 m/s and a kinetic energy of 1.80 J, just as in the prior example. The approach seems sound. Let’s run it for longer. I’ll say that 10% of the mass of the craft is propellant. Let’s see what happens:

The drive will work for 100,000 s (which makes sense since we eject 1 g/s, and we started with a reaction mass of 100,000 g — 10% of 1,000,000 g), or 27.8 h. It will reach a velocity of 105 m/s (235 mph, 378 km/h). The acceleration looks constant. Is it? I’ll add in a line in the program to calculate the acceleration, by taking the new velocity minus the old velocity, and dividing by the timestep: a = (vfvi) / Δt.

Interesting. I actually thought the acceleration would be roughly constant over this period. It starts off with an acceleration of 1.00×10−3 m/s², and ends with an acceleration of 1.11×10−3 m/s². I assume that this rate of increase is because the craft is 10% lighter by the end, so a constant thrust will accelerate it more. We can calculate the force by taking the acceleration and multiplying by the mass at that time (F=m1a).

Well, it certainly looks constant. That matches what I’ve read about ion thrusters. It starts and ends at 1.00 N of thrust. That means that the increasing acceleration is due to the progressive lightening of the craft. Our 500-W drive is able to produce 1 N of thrust when ejecting 1 g/s. That’s interesting. I wonder how it depends on the parameters (it doesn’t seem to depend much on the mass of the ship), and if we can derive a mathematical relationship — I’ll investigate that later. Finally, let’s explore the issue that sparked my interest in this problem in the first place, the kinetic energy. It just takes another line of code to calculate the kinetic energy as EK = ½m1v1²

As expected, it looks parabolic. Of course, the velocity is not increasing linearly, as we know the acceleration is actually increasing over time. But the mass is decreasing, too.  It reaches a kinetic energy of 5.00 MJ (5.00×106 J). We already know that this must only be a small portion of the total kinetic energy. Let’s find the kinetic energy of the propellant. This is slightly more complex to code, since each at each time step we need to add in the (unchanging) kinetic energy of the new bit of propellant that’s been ejected, rather than looking at the changing velocity and mass of the craft at that instant. And I’ll track the sum of the two kinetic energies as well.

I notice several interesting things about this graph. As expected, the kinetic energy of the exhaust is much more than the kinetic energy of the spacecraft (note the 10-fold increase in scale from the previous graph). The total kinetic energy does appear to be increasing linearly: the kinetic energy of the exhaust ends at 4.50×107 J, and after adding in the 0.500×107 J of spacecraft kinetic energy, we get a total of 5.00×107 J. In a way, this is a relief, because we’ve run our 500-W power source for 1.00×105 s, and multiplying the two gives 5.00×107 J. This gives me some confidence that the coding is correct. I also note that the kinetic energy of the exhaust increases, but that rate of increase slows over time (it curves away from the total line). Of course, this is to be expected, because the spacecraft’s kinetic energy is increasing at an increasing rate (and there is a fixed amount of total energy), but I’m still interested in understanding why this is happening. Let’s look at the exhaust velocity (v2 in our derivation above).

The exhaust velocity starts and ends at -1000 m/s (for graphing purposes, I flipped the sign to make it positive). This is also a relief, because our engine power was based on the earlier model which used a  craft that ejected propellant at -1000 m/s. If the exhaust velocity is constant, and we’re ejecting a constant amount of mass each second, why is the rate of kinetic energy dropping off? It’s because the constant exhaust velocity is relative to the spacecraft. But the velocity is less, relative to our rest observational frame. To begin with, if the craft is motionless and ejecting propellant at -1000 m/s, we’ll also measure that velocity. But if the craft were moving at, say, 10 m/s, then the propellant would be ejected at -990 m/s relative to us. It would have less kinetic energy than the propellant that was previously ejected. (From the observer’s point of view, the propellant had been moving forward at 10 m/s, so even though it was pushed out with the same energy, it didn’t achieve as high a speed moving backward.)

If the spacecraft’s kinetic energy were truly parabolic, it would eventually cross any arbitrary line, including the line of total kinetic energy. Since this presumably can’t happen, I’ll want to explore what actually does occur — in the next post.

# Rocket Propulsion, Part 1

The kinetic energy of an object (the energy it has due to its motion) is equal to ½mv2; that is, half the product of its mass and the square of its velocity. But a question I saw recently pointed out some interesting consequences of this, especially in regards to travel in space.

We often think about constant forces accelerating objects. If the force is constant, so is the acceleration (by F=ma). And if the acceleration is constant, then the velocity increases at a constant rate. It would seem that it would take a constant supply of energy to supply a constant force (that is, the total energy used would increase linearly). But the kinetic energy of the object would increase quadratically. If it took x amount of energy to increase the velocity by a certain amount, it would take 4x to double that and 9x to triple that. Where is this extra energy coming from?

Let’s construct a scenario and put in some test numbers. I’m going to stick to velocities well below the speed of light and analyze the situation according to classical mechanics, not taking relativity into account. We’ll start with summarizing the basic physics involved. If you’re already familiar with the basics of momentum, feel free to skip ahead.

Of course, forces can’t be applied arbitrarily. There has to be a source. Even if a craft has plenty of power, conservation of momentum affects how it can move. Recall that (linear) momentum is the product of mass and velocity (p = mv); each object has a momentum. The sum of the momenta of a closed system (m1v1 + m2v2 + m3v3 + …) is conserved and cannot change; the total momentum of all the objects in the system before an event must be equal to the total momentum afterwards. For an isolated spacecraft, mivimfvf. If a spacecraft is initially at rest, then vi = 0 and so momentum p = 0. If it is to start moving with a positive velocity vf, then it will have a positive momentum mvf. Of course, 0≠mvf. (Or even if it’s not at rest, in order to accelerate, the velocity on the right will be higher than on the left, and therefore, so will the momentum). We either need to add a positive mv term to the left side of the equation or a negative mv term to the right (and since mass is positive, it will have to be an object with positive velocity before or with negative velocity after).

A positive term on the left would represent a second object with positive velocity, flying toward the spacecraft. This could be a laser beam on Earth fired towards the craft, particles from the solar wind pushing a solar sail, and so on. However, it’s not the type of propulsion we are looking for, since it has to be produced by an external source and can only push the craft away from the source.

The other way for momentum to be conserved is to add a negative mv term on the right side of the equation. That is, something must be pushed backwards. On Earth, this could be air (for an airplane), or the Earth itself (for a car, or for people walking). In the near vacuum of space, there is not much to push against, and so any mass pushed backward must come from the spacecraft itself. This is called reaction mass. Let’s designate the craft mass 1, and the propellant mass 2. The momentum equation would thus be (m1 + m2)vim1v1f + m2v2f.

In this example, the questioner suggested the object be a spacecraft with an ion thruster, which uses an electric field to accelerate ions out the rear of the craft, pushing the craft forward. These use only small amounts of mass, and generate a roughly constant thrust. This sounds perfect for our initial exploration, since we’ll get a constant force and acceleration and the mass of the craft shouldn’t change significantly.

Finally, it’s time to put in some test numbers. We’ll start with a spacecraft mass of 1000 kg to start — of course, it will get lighter over time. It will start at rest. Let’s assume it can eject 1 g (=0.001 kg) of mass each second at -1000 m/s. To make it easier, we’ll break this down into 1-second intervals. At t=0, the craft is 1000 kg and is moving at 0 m/s. After the first second, it ejects 1 g of mass at -1000 m/s. This is a momentum of -1 kg·m/s. The craft now still has a mass of 1000 kg (technically 999.999 kg), and since the initial momentum was 0, it must have a momentum of 1 kg·m/s. Dividing by the new mass gives a velocity of 0.001 m/s, and a kinetic energy of 0.0005 J = 0.500 mJ.

We’ll use the momentum equation from above. We know the mass of the craft and the propellant, and the initial velocity. The propellant velocity vpf (the propellant velocity) is vi – 1000 (starting at the velocity of the craft, and being ejected at a relative -1000 m/s). If mc is the mass of the craft, mp is the mass of the propellant, and mc’ = mcmp is the new mass of the craft, then

$m_c v_i = m_{c^\prime}v_f + m_p v_{pf}$

$m_c v_i - m_p v_{pf} = m_{c^\prime} v_f$

$\frac{m_c v_i - m_p v_{pf}}{m_{c^\prime}} = v_f$

Or put simply, we need to find the initial momentum, subtract the final momentum of the propellant to get the final momentum of the craft, and divide by the new mass to find the final velocity (at each step).

Continuing in this fashion, at t = 2, the mass is down to 999.998 kg, velocity is 0.002 m/s, and kinetic energy is 2.00 mJ. Let’s look at the first minute:

Time / s Mass / kg Velocity / (m/s) Acceleration / (m/s²) Kinetic energy / mJ
0 1000 0 0 0
10 1000 0.0100 0.00100 50
20 1000 0.0200 0.00100 200
30 1000 0.0300 0.00100 450
40 1000 0.0400 0.00100 800
50 1000 0.0500 0.00100 1250
60 1000 0.0600 0.00100 1800

Over this minute, we can see that the craft stays roughly the same mass (it drops to 999.94 kg) and has roughly constant acceleration (it actually increases from 0.001000000 to 0.00100006 m/s²). But the kinetic energy is going up quadratically, as one would expect with constant acceleration and subsequent linearly increasing velocity. What are we missing? Of course, we’ve ignored the propellant that’s been ejected. It may have very low mass compared to the craft, but it has proportionally high velocity, and remember that kinetic energy is ½mv2. In that first second when the 1,000,000-g craft reached 0.001 m/s and a kinetic energy of 0.500 mJ, we ejected 1 g of reaction mass at -1000 m/s. With this high velocity, this mass has a kinetic energy of 500,000 mJ! That is one million times the kinetic energy of the craft. The combined kinetic energy is 500,000.5 mJ. Let’s look at what happens over the first minute, remembering that we have to add up the kinetic energy of each block of mass that’s been ejected up to that point.

Time / s Kinetic energy / mJ
Spacecraft Propellant Total
0 0 0 0
10 50 4,999,955 5,000,005
20 200 9,999,810 10,000,010
30 450 14,999,565 15,000,015
40 800 19,999,220 20,000,020
50 1250 24,998,775 25,000,025
60 1800 29,998,230 30,000,030

This, then, is the resolution to the paradox. The vast majority of the kinetic energy goes to the propellant, not to the spacecraft. The spacecraft does accelerate and its kinetic energy does go up quadratically, but it is only a tiny fraction of the total kinetic energy (which in this example, goes up by 5,000,005 mJ every 10 seconds).

Of course, we used a model where the amount of mass ejected was very small relative to the craft. What will happen over longer periods, or if we use more mass? Any quadratically increasing function will eventually pass a linear one. In order to explore these issues, we’ll have to make our model a bit more rigorous (in particular, we have to clarify how our drive works). I’ll look at these in the next post.

# Motion Along a Curve, Part 3

See Part 1 and Part 2, where I discuss trying to find equations of motion for a ball rolling along a track defined by an arbitrary function y(x). So far, I’ve worked out my general approach, and tested it on the very simple case of an inclined line. I do want to tackle some complex functions, but first I want to summarize the method, and to incorporate some changes I learned as I tried my first solution.

I had initially stated I wanted to start with x(0) = 0, to keep it simple. But I didn’t end up needing this restriction. My equation for velocity used the initial yi, which we got from plugging in the initial value of x. Also, we used the x(0) = 0 to find the constant after integrating. But these happened later in the process.

I also had started with an initial velocity of zero. That did make a difference. But looking back on part 1, I don’t think it would complicate the equation too much, and if it’s zero it will be an extra term to just drop out. Let’s go back to the conservation of energy equation, and keep vi in this time.

$\displaystyle U_i + K_i = U + K$

$\displaystyle mgy_i + \frac{1}{2}m{v_i}^2 = mgy + \frac{1}{2}mv^2$

$\displaystyle gy_i - gy = \frac{1}{2}(v^2 - {v_i}^2)$

$\displaystyle 2g(y_i - y) = v^2 - {v_i}^2$

$\displaystyle 2g(y_i - y) + {v_i}^2 = v^2$

$\displaystyle v = \sqrt{2g(y_i-y)+{v_i}^2} \blacktriangleleft$

I’m leaving the ± in this time. Strictly speaking, I’m not treating this as the magnitude of the vector, since magnitudes must be positive. Rather, I want to consider a velocity vector that can point ether forwards or backwards along the direction of the curve. I’m going to allow motion in both directions, not just forward.

Recall the graph showing components of the velocity vector:

Now I want to find the x-component. As I discovered last time, I don’t need to bother with the y-component — once I find an equation for x(t), I can use that directly to obtain y(t). The x-component will be

$\displaystyle v_x = v \frac{1}{\sqrt{(y')^2+1}}$

$\displaystyle v_x = \pm \sqrt{\frac{2g(y_i - y) + {v_i}^2}{(y')^2+1}}$

$\displaystyle \frac{dx}{dt} = \pm \sqrt{\frac{2g(y_i - y) + {v_i}^2}{(y')^2+1}} \blacktriangleleft$

Since I will have y and y′ in terms of x, I would need to rearrange to solve the differential equation. Let’s see how far I can take the general case:

$\displaystyle \pm \sqrt{\frac{(y')^2+1}{2g(y_i-y)+{v_i}^2}} \, dx = dt$

$\displaystyle \pm \int \sqrt{\frac{(y')^2+1}{2g(y_i-y)+{v_i}^2}} \, dx = \int dt$

$\displaystyle \pm \int \sqrt{\frac{(y')^2+1}{2g(y_i-y)+{v_i}^2}} \, dx = t \blacktriangleleft$

where I did not include a constant of integration on the right side, since it can be absorbed into the constant that the left integral will produce.

So, the general approach should be as follows: Given our equation y(x), find y′(x). Plug in those expressions, plug in the initial velocity, and plug in the initial height y[x(0)]. Integrate, and solve for x in terms of t to get x(t), then plug that into y(x) to get y(t). I’ll test if this approach can actually work in the next post.

# Motion Along a Curve, Part 2

I would like to find equations for the position of a small ball rolling along a curve given by an equation, rather than a flat incline. In Part 1, I introduced the problem and set up a test case. Now let’s try to actually make some headway.

Last time, we showed that

$\displaystyle v = \sqrt{2g(y_i - y)}$

where v is the magnitude of the velocity vector v, g is the acceleration due to gravity, and y is the height. I note that if y>yi, then the quantity inside the radical is negative and v will not be real — which makes sense, since the ball should never be able to roll higher than its starting position. (Recall that I am starting with an initial velocity of zero). Since we’re assuming that the ball rolls along the function, its direction will be given by the derivative of y with respect to x. I do see a potential problem with notation arising. The value of y depends on x (by definition, since that’s how I’m defining the curve). But x and y both depend on t if we’re thinking about the motion of the ball. I want to define our variables carefully. We should think of x as a function x(t), and y as a function of x, so y(x) = y[x(t)]. I am going to reserve the prime notation y′ specifically to refer to the derivative with respect to x.

$\displaystyle y' = \frac{dy}{dx}$

We know that

$\displaystyle \mathbf{v} = \langle r, \angle\theta \rangle$

The magnitude will be v. If we take the tangent line to the curve and sketch a triangle, we can see that for a change in x of 1, y changes by y′, so tan θ = y′/1. I illustrated this below with an arbitrary curve.

$\displaystyle \mathbf{v} = \langle v, \angle \arctan y' \rangle$

Let’s now find the component vectors vx and vy, where of course vx = r cos θ and vy = r sin θ.

$\displaystyle \mathbf{v} = \langle v \cos \arctan y', v \sin \arctan y' \rangle$

There’s actually no need for these arctangents. If I go back to the earlier sketched triangle, I can see that the hypotenuse is $\sqrt{(y')^2 +1}$ and so I can read the cosine and sine directly off that triangle (1 over the hypotenuse and y’ over the hypotenuse, respectively):

$\displaystyle \mathbf{v} = \left\langle v \frac{1}{\sqrt{(y')^2 + 1}}, v \frac{y'}{\sqrt{(y')^2 +1}} \right\rangle$

$\displaystyle \mathbf{v} = \frac{v}{\sqrt{(y')^2+1}}\langle 1, y' \rangle$

Let’s substitute in our expression for the magnitude, v:

$\displaystyle \mathbf{v} = \sqrt{\frac{2g(y_i-y)}{(y')^2+1}} \langle 1, y' \rangle \blacktriangleleft$

I believe that this should be the general solution, given the constraints (ball must roll along the curve, uniform gravity, no friction, starting x and v of 0). Let’s pause to make some observations. The expression in the numerator is the magnitude of the velocity which we found in the last post. The potential/kinetic energy relationship is obvious, with the product of g and the change in height from potential energy, and the times two and the square root being inverses from kinetic energy. The <1, y′> ensures the proper direction, and the denominator adjusts for the length of that direction vector (it basically makes it a unit vector).

This equation is so far pretty straightforward. However, in order to find equations for x and y in terms of time t, we’d have to plug in the expression for y and y′, then integrate the velocity. Let’s try this with what should be a very simple example, the incline (line) from the prior entry:

$\displaystyle y = -\frac{x\sqrt{3}}{3} + 1$

$\displaystyle y' = -\frac{\sqrt{3}}{3}$

$\displaystyle y_i = y(0) = 1$

Plugging these into our equation for v gives:

$\displaystyle \mathbf{v} = \sqrt{\frac{2g(y_i-y)}{(y')^2+1}} \langle 1, y' \rangle$

$\displaystyle \mathbf{v} = \sqrt{\frac{2g[1-(-\frac{x\sqrt{3}}{3}+1)]}{(-\frac{\sqrt{3}}{3})^2+1}} \left\langle 1, -\frac{\sqrt{3}}{3} \right\rangle$

I am starting to regret setting my initial angle to π/6. There was no benefit to choosing a “nice” value for the angle; it would have been better to choose a slope with an integer value. However, it’s too late to go back and change it now. Let me see if I can simplify this.

$\displaystyle \mathbf{v} = \sqrt{\frac{2g\frac{x\sqrt{3}}{3}}{\frac{1}{3}+1}} \left\langle 1, -\frac{\sqrt{3}}{3} \right\rangle$

$\displaystyle \mathbf{v} = \sqrt{\frac{2gx\sqrt{3}}{\frac{4}{3}\cdot 3}} \left\langle 1, -\frac{\sqrt{3}}{3} \right\rangle$

$\displaystyle \mathbf{v} = \sqrt{\frac{gx\sqrt{3}}{2}} \left\langle 1, -\frac{\sqrt{3}}{3} \right\rangle$

Let’s split this into its component vectors, and switch to magnitudes:

$\displaystyle \mathbf{v} = \mathbf{v_x} + \mathbf{v_y}$

$\displaystyle \mathbf{v} = v_x \boldsymbol{\hat{\i}} + v_y \boldsymbol{\hat{\j}}$

$\displaystyle v_x = \sqrt{\frac{gx\sqrt{3}}{2}}$

and

$\displaystyle v_y = -\sqrt{\frac{gx\sqrt{3}}{2}} \cdot\frac{\sqrt{3}}{3} = -\sqrt{\frac{gx\sqrt{3}}{6}}$

Let’s now try integrating the first one:

$\displaystyle \frac{dx}{dt} = \sqrt{\frac{gx\sqrt{3}}{2}}$

$\displaystyle \frac{1}{\sqrt{x}} \, dx = \sqrt{\frac{g\sqrt{3}}{2}} \, dt$

$\displaystyle \int x^{-\frac{1}{2}} dx = \int \sqrt{\frac{g\sqrt{3}}{2}} \, dt$

$\displaystyle \frac{x^{\frac{1}{2}}}{\frac{1}{2}} + C_1 = t \sqrt{\frac{g\sqrt{3}}{2}} + C_2$

$\displaystyle 2\sqrt{x} = t \sqrt{\frac{g\sqrt{3}}{2}} + C_3$

$\displaystyle \sqrt{x} = t \sqrt{\frac{g\sqrt{3}}{8}} + C_3$

This is looking promising. Before I square both sides, let’s deal with the constant. Since we’re starting at a position x=0 at t=0,

$\displaystyle 0 = 0 + C_3$

so clearly the constant is zero. (This won’t be so simple for the y-component.)

$\displaystyle \sqrt{x} = t \sqrt{\frac{g\sqrt{3}}{8}}$

I will square both sides, realizing that I might be introducing extraneous solutions. In fact, t cannot be negative, since the square root of x is nonnegative.

$\displaystyle x(t) = \frac{g\sqrt{3}}{8} t^2, \, t\geq 0 \blacktriangleleft$

Before tackling the y-direction, let’s test this answer. Let’s plug in our value from the last post for t, for the time when it reaches the x-axis.

$\displaystyle x\left(\sqrt{\frac{8}{g}}\right) = \frac{g\sqrt{3}}{8} \left(\sqrt{\frac{8}{g}}\right)^2 = \frac{g\sqrt{3}}{8} \cdot \frac{8}{g} = \sqrt{3}$

This is indeed the value of x at the base of the triangle (recall height 1, width √3, hypotenuse 2 for a 30-60-90 triangle)! Now let’s find y. Let’s start with our equation for the y-component of the velocity:

$\displaystyle v_y = -\sqrt{\frac{gx\sqrt{3}}{6}}$

This one is going to be a little more complicated.

$\displaystyle \frac{dy}{dt} = -\sqrt{\frac{gx\sqrt{3}}{6}}$

I already see a complication. My equation has x, not y. We of course have our starting equation for y in terms of x, and it would be easy to solve that for x and substitute it in. However, if we want to try this on other functions, especially those that aren’t injective (that is, where multiple values of x might result in the same value of y — in a parabola, or sine curve for instance), then solving will be difficult. I think it would be better to try to use the chain rule.

$\displaystyle \frac{dy}{dx} \frac{dx}{dt} = -\sqrt{\frac{gx\sqrt{3}}{6}}$

$\displaystyle -\frac{\sqrt{3}}{3} \frac{dx}{dt} = -\sqrt{\frac{gx\sqrt{3}}{6}}$

I just realized something. We already have $\frac{dx}{dt}$. It was the vx we started with to find x(t) above!

Wait. That expression is in terms of x, and we have no y’s left. Substituting that in will just result in an identity.

$\displaystyle -\frac{\sqrt{3}}{3} \sqrt{\frac{gx\sqrt{3}}{2}} = -\sqrt{\frac{gx\sqrt{3}}{6}}$

$\displaystyle -\sqrt{\frac{gx\sqrt{3}}{6}} = -\sqrt{\frac{gx\sqrt{3}}{6}}$

This is bad news, because I really don’t think that solving the starter equation for y will be a practical solution for most functions. But if we don’t have a y or a dy, then we’re not going to get an equation for y.

Scratch that. There’s no need to do any of this! We just derived an equation x(t). We already know y(x). So we can find y[x(t)]!

$\displaystyle y(x) = -\frac{x\sqrt{3}}{3} + 1$

$\displaystyle x(t) = \frac{g\sqrt{3}}{8}t^2$

So, plugging in the equation for x in,

$\displaystyle y(t) = -\frac{\sqrt{3}}{3} \left(\frac{g\sqrt{3}}{8}\right) t^2 + 1$

$\displaystyle y(t) = -\frac{g}{8}t^2 + 1, \, t\geq 0 \blacktriangleleft$

This looks right. It’s of a similar form to our equation for x(t), but with a factor of -1 instead √3 (which follows from our 30-60-90 triangle, and y is decreasing. And it has a +1 to raise the initial height. Just to be sure, let’s check it with our test value for t:

$\displaystyle y\left(\sqrt{\frac{8}{g}}\right) = -\frac{g}{8}\left(\sqrt{\frac{8}{g}}\right)^2 + 1 = -1 + 1 = 0$

And that’s correct, because we specifically picked the test case for when the ball reaches the x-intercept. So our final solution is

$\displaystyle x(t) = \frac{g\sqrt{3}}{8} t^2, \,t\geq0 \blacktriangleleft$

$\displaystyle y(t) = -\frac{g}{8}t^2 + 1, \, t\geq 0 \blacktriangleleft$

The method works! In future explorations, it would be interesting to try this for a more complicated starter equation.

The restriction t≥0 is interesting. If it weren’t there, the equations for the motion of the ball would still be valid, since it would represent the ball being launched backwards up the slope, coming to rest for an instant at (0,1) at t=0, then rolling back down the slope. I’m not surprised that it is there, though, because my method explicitly assumed that the ball’s velocity vector points in the direction <1, y′>, which points in the forward x direction. Let’s trace back the restriction on t. It comes from where it times a constant equals √x, which must be positive. Why did we have √x, as opposed to ±√x? That came from the equation for the magnitude of the velocity v. That in turn came from the fact that K = ½mv2. Should we have kept the negative solution? No, because v represents the magnitude, and the magnitude is nonnegative. Still, a negative v would correspond to a vector pointing in the opposite direction, when put into our component equation. I believe that those solutions should be allowed. I could justify it by saying that v was not the magnitude, but rather ±magnitude, to allow the vector to point in either direction. Now that I think about it, this will occur for every problem we attempt with this method. Since the laws of physics and therefore the equations of motion should be symmetric over time, this means that if the ball is starting at rest at time = 0, we could “run the film in reverse” to see what path it took to get to the starting point. This is an interesting development, and I wonder if expanding the method in this manner would allow the ball to roll backwards at other times, too — say if it were to roll down a valley, then up the other side. If we allow these “negative magnitudes,” might we be able to get an equation which expresses the motion as it rolls back down into the valley and back up the original side?

# Motion Along a Curve, Part 1

Earlier this week, I saw an intriguing problem suggested, one involving motion along a curve. The idea is this: in elementary mechanics, we learn about how an object will move when placed on an inclined surface. If we restrict motion to two dimensions, this surface can be represented by a line. With gravity as the propelling force, can we find a general approach to deriving its motion along an arbitrary curve? In other words, given a function f(x), can we find functions describing the position (and therefore velocity and acceleration) over time?

Here are my initial thoughts. One, I am sure that this type of problem has been analyzed before. That doesn’t matter; I am going to try figuring it out myself. Two, I think we can use conservation of energy to determine the object’s velocity at any point. However, that will depend on knowing its height, and I don’t know how this approach would handle the object flying up off a short hill, or dropping away from a cliff. So three, I am going to restrict the object’s motion to the curve — consider it a track, rather than a 1D surface. That means that four, I can use the derivative to find the direction, but that means the object can only move forward. Five, I have to ignore friction. My conservation of energy approach will require that all energy be potential or kinetic; I won’t know hot to deal with losses to friction. Six, I am using constant downward gravitation. Seven, I am starting with initial position x=0 and initial velocity of 0. I don’t believe this method requires it, but I am already concerned about the complexity of the math. And seven, I am planning to obtain a velocity equation and integrate to find position, but I don’t know if the equation will have an analytic solution (or if I will know how to integrate it).

These are some significant limitations, of course. If this approach works, I suspect that it could be expanded to deal with several of these. For instance, we could keep in terms for initial position and velocity. Gravity or whichever force could be represented by a vector field rather than a constant force; this would make the potential energy term much more complex. I also wonder if I could incorporate friction by adding in an extra term, but I think this would turn the equation into a more complex differential equation.

Let’s start with a simple case that we can solve using conventional means, so we’ll have an answer to check later. I can also use this to check my conservation of energy approach. Let’s have the object rolling down an incline of 30° (π/6), starting at a height of 1. The ball will roll down to the right (I find it more intuitive to imagine a rolling ball rather than a frictionless object sliding, especially if the surface is curved). The equation for this surface will be

$\displaystyle y(x) = -\frac{x\sqrt{3}}{3}+1 \, \bullet$

where the slope is -tan(π/6).

You can see that it would form a right triangle with height 1, length √3, and hypotenuse 2. The slope is therefore 1/√3, or √3/3. Perhaps it would have been better to use an incline of 60° (π/3) so that the slope would be -√3, but I like this one better. Let’s try the conservation of energy approach to see what the velocity would be when the ball reaches the bottom (that is, y = 0). The total energy E is the sum of the potential energy U and the kinetic energy K. This should remain constant, so

$\displaystyle E_i = E_f$

$\displaystyle U_i + K_i = U_f + K_f$

We know that potential energy is given by

$\displaystyle U = mgh$

assuming constant gravitation g (the acceleration due to gravity), with mass m and height h; and kinetic energy is given by

$\displaystyle K = \frac{1}{2}mv^2$

where v is the magnitude of the velocity vector v. I confess I am not very proficient with vectors, but clearly the ball is moving in two dimensions and we will need both components in the x and y directions. I am going to try to not be sloppy but to carefully think about what we mean by v or v every time I write it.

We assume that the initial velocity, and therefore kinetic energy, is zero.

$\displaystyle U_i + K_i = U_f + K_f$

$\displaystyle mgh_i + \frac{1}{2}m\cdot0^2 = mgh_f + \frac{1}{2}mv_f^2$

The height is y, and we can divide both sides by the nonzero mass m:

$\displaystyle gy_i = gy_f + \frac{1}{2}v_f^2$

Let’s solve for v:

$\displaystyle 2g(y_i-y_f)=v_f^2$

$\displaystyle v_f = \sqrt{2g(y_i-y_f)} \blacktriangleleft$

This should be the general case, so if this checks out, we can use this equation to develop the method. As a check, if the units of g are m/s² and y is in m, the units will be m/s, appropriate for velocity. Also note that this is simply the magnitude of the velocity vector (the speed); the actual vector could point in any direction if all we know is conservation of energy. That’s why I earlier constrained the ball to move along the curve — that will give us our direction. Plugging in our initial and final values of y, we get

$\displaystyle v_f = \sqrt{2g(1-0)} = \sqrt{2g}$

Next, to verify this solution, let’s solve using the traditional free-body approach. Gravity applies a force downwards of W=mg. This can be split into two components, one perpendicular to the surface mg cos(θ), and one along the surface mg sin(θ). The perpendicular force will be balanced out by the equal and opposite normal force, leaving a force of mg sin(θ). This is simple motion in one dimension, and let me introduce a variable s to represent its position along the surface (so I don’t generate confusion with our x’s and y’s). Think of it as placing a tape measure from the bottom of the incline up to the starting point. We know that for a given constant acceleration a,

$\displaystyle v = \int \! a \, dt = at + v_0$

$\displaystyle s = \int \! v \, dt = \int \! (at+v_0) \, dt = \frac{1}{2}at^2 + v_0 t + s_0$

By Newton’s second law, F = ma, so a = g sin(θ) and will be negative, since it’s pointing towards decreasing s (since this is motion in one dimension, we can use the sign to expression direction rather than need vectors). Initial velocity is zero, and initial position is 2 (the distance along the incline, if the bottom is at s=0).

$\displaystyle 0 = -\frac{1}{2}g \sin(\frac{\pi}{6}) \cdot t^2 + 0t + 2$

$\displaystyle \frac{1}{2} g \left( \frac{1}{2}\right)t^2 = 2$

$\displaystyle \frac{gt^2}{4} = 2$

$\displaystyle t^2 = \frac{8}{g}$

$\displaystyle t = \pm \sqrt{\frac{8}{g}}$

I can discard the negative solution, since I’m interested in the behavior starting from time = 0, not how it might have been launched before to come to a temporary stop at time = 0. Let’s plug this time to get our velocity from v = at (since the initial velocity is zero):

$\displaystyle v = at = -g\left(\frac{1}{2}\right)\sqrt{\frac{8}{g}} = - \sqrt{\frac{g^2}{4}}\sqrt{\frac{8}{g}} = -\sqrt{2g}$

This velocity is negative since it points downwards along the s direction, along the incline. Its magnitude will be the absolute value which is what we obtained earlier using the conservation of energy method. Armed with this success, and with some test values, we’ll be ready to work on the actual problem in the next post.

# Simple Harmonic Motion

I came across an interesting approach to a problem involving a spring the other day, and it made me wonder if I could derive the equations for simple harmonic motion. The basic principle of a spring is Hooke’s law, that the spring exerts a restoring force proportional to its displacement from equilibrium position: F = –kx, where k is the spring constant representing the stiffness of the spring. I wondered if I could derive the formulas for its motion, and where the sinusoidal components would appear from, since clearly there are no trigonometric functions in this basic equation. To begin with, from Newton’s second law, F = ma:

$\displaystyle ma = -kx$

$\displaystyle a = -\frac{kx}{m}$

We know that acceleration is the derivative of velocity and velocity is the derivative of position:

$\displaystyle a = \frac{dv}{dt}$

and

$\displaystyle v = \frac{dx}{dt} .$

We can rewrite these as a dt = dv and v dt = dx, and integrate. However, our acceleration equation has a in terms of position x, not time t, so I would prefer to integrate with respect to x (as at this point I don’t know how x depends on t — that’s what we’re trying to derive!).

$\displaystyle a \, dt = dv$

$\displaystyle a \left( \frac{dx}{v} \right) = dv$

$\displaystyle a \, dx = v \, dv$

This equation is well-known and could have been used as a starting point, but it’s nice to derive it from the most basic principles. We can now substitute our expression in for a and integrate both sides.

$\displaystyle -\frac{kx}{m} \, dx = v \, dv$

$\displaystyle \int \! \left( -\frac{kx}{m} \right) \, dx = \int \! v \, dv$

$\displaystyle -\frac{k}{m} \int \! x \, dx = \int \! v \, dv$

$\displaystyle -\frac{k}{m} \left( \frac{x^2}{2} + C_1 \right) = \frac{v^2}{2} + C_2$

Let’s solve this for v:

$\displaystyle \frac{v^2}{2} + C_2 = -\frac{k}{m} \left( \frac{x^2}{2} + C_1 \right)$

$\displaystyle \frac{v^2}{2} + C_2 = -\frac{kx^2}{2m} + C_3$

$\displaystyle \frac{v^2}{2} = -\frac{kx^2}{2m} + C_4$

$\displaystyle v^2 = -\frac{kx^2}{m} + C_5$

$\displaystyle v = \pm \sqrt { -\frac{kx^2}{m} + C_5 }$

There’s still no sign of sine appearing anywhere (I know; I couldn’t resist). I’m also a bit uncomfortable with the plus/minus sign, but velocity can be both positive and negative, so I don’t think that discarding the negative solution is justified. I want to integrate again to remove the velocity, but we face a similar problem as before: we have velocity in terms of position x, not time t. We can rearrange the differential again:

$\displaystyle v = \frac{dx}{dt}$

to get v on the dx side:

$\displaystyle dt = \frac{dx}{v} .$

Let’s substitute in our expression for v, and integrate:

$\displaystyle dt = \frac{dx}{\pm \sqrt { -\frac{kx^2}{m} + C_5 } }$

$\displaystyle \int dt = \int \! \frac{dx}{\pm \sqrt { -\frac{kx^2}{m} + C_5 } }$

This is not pretty, and I’m not sure how to integrate the right side. I believe I can pull out the plus/minus as a constant.

$\displaystyle t + C_6 = \pm \int \! \frac{dx}{\sqrt { -\frac{k}{m} x^2 + C_5 } }$

I had to consult a table of integrals for assistance with this integration. Of course, it can be solved with trignometric substitution! All of a sudden, we have a hint of where a sine function will appear. We need to get the radical in the form $\sqrt{a^2 - b^2 x^2 }$:

$\displaystyle \sqrt{C_5 - \frac{k}{m} x^2 }$

$\displaystyle a = \sqrt{C_5} = C_7$

$\displaystyle b= \sqrt{\frac{k}{m}}$

Now, we set

$\displaystyle x = \frac{a}{b} \sin \theta :$

$\displaystyle x = C_7\sqrt{\frac{m}{k}}\sin \theta$

Let’s now take our original integral back, and start substituting.

$\displaystyle t + C_6 = \pm \int \! \frac{dx}{\sqrt { C_5 -\frac{k}{m} x^2 } }$

$\displaystyle t + C_6 = \pm \int \! \frac{dx}{\sqrt { {C_7}^2 -\frac{k}{m} x^2 } }$

$\displaystyle t + C_6 = \pm \int \! \frac{dx}{\sqrt { {C_7}^2 -\frac{k}{m} \left( C_7\sqrt{\frac{m}{k}}\sin \theta \right) ^2 } }$

This expression is complicated, but I can already see how it will simplify out.

$\displaystyle t + C_5 = \pm \int \! \frac{dx}{\sqrt{ {C_7}^2 - \frac{k}{m} \cdot {C_7}^2 \cdot \frac{m}{k} \cdot \sin^2 \theta} }$

$\displaystyle t + C_5 = \pm \int \! \frac{dx}{\sqrt{ {C_7}^2 - {C_7}^2 \sin^2 \theta } }$

$\displaystyle t + C_5 = \pm \int \! \frac{dx}{\sqrt{{C_7}^2 (1 - \sin^2 \theta)} }$

$\displaystyle t + C_5 = \pm \frac{1}{C_7} \int \! \frac{dx}{\sqrt{1 - \sin^2 \theta} }$

Let’s convert the dx to a dθ:

$\displaystyle x = C_7\sqrt{\frac{m}{k}}\sin \theta$

$\displaystyle \frac{dx}{d\theta} = C_7\sqrt{\frac{m}{k}}\cos \theta$

$\displaystyle dx = C_7\sqrt{\frac{m}{k}}\cos \theta \, d\theta$

Substituting this back in yields:

$\displaystyle t + C_5 = \pm \frac{1}{C_7} \int \! \frac{C_7\sqrt{\frac{m}{k}}\cos \theta \, d\theta }{\sqrt{1 - \sin^2 \theta} }$

$\displaystyle t + C_5 = \pm \frac{C_7}{C_7}\sqrt{\frac{m}{k}}\int \! \frac {\cos \theta}{\sqrt{1 - \sin^2 \theta}} d\theta$

$\displaystyle t + C_5 = \pm \sqrt{\frac{m}{k}}\int \! \frac{\cos \theta}{\sqrt{\cos^2 \theta}} d\theta$

$\displaystyle t + C_5 = \pm \sqrt{\frac{m}{k}}\int \! \frac{\cos \theta}{\cos \theta} d\theta$

$\displaystyle t + C_5 = \pm \sqrt{\frac{m}{k}}\int d\theta$

$\displaystyle t + C_5 = \pm \theta \sqrt{\frac{m}{k}} + C_8$

$\displaystyle t + C_9 = \pm \theta \sqrt{\frac{m}{k}}$

Remembering our equation for x in terms of θ:

$\displaystyle C_7\sqrt{\frac{m}{k}}\sin \theta = x$

$\displaystyle \sqrt{\frac{m}{k}}\sin \theta = C_{10}x$

$\displaystyle \sin \theta = C_{10}x\sqrt{\frac{k}{m}}$

$\displaystyle \theta = \arcsin \left(C_{10}x\sqrt{\frac{k}{m}} \right)$

Finally, we can substitute this back in for θ and solve for x:

$\displaystyle t + C_9 = \pm \left[ \arcsin \left(C_{10}x\sqrt{\frac{k}{m}} \right) \right] \sqrt{\frac{m}{k}}$

$\displaystyle \pm t\sqrt{\frac{k}{m}} + C_{11} = \arcsin \left(C_{10}x\sqrt{\frac{k}{m}} \right)$

The plus-minus is driving me crazy, but I can’t justify dropping it — even though I want to.

$\displaystyle \sin \left( \pm t\sqrt{\frac{k}{m}}+ C_{11} \right) = C_{10}x\sqrt{\frac{k}{m}}$

Actually, I think I can get rid of it now! Since sin(-x) = -sin(x) = sin(x+π), I can incorporate the +π into the constant!

$\displaystyle C_{10}x\sqrt{\frac{k}{m}} = \sin \left( t\sqrt{\frac{k}{m}} + C_{12} \right)$

$\displaystyle x = C_{13} \sin \left( t\sqrt{\frac{k}{m}} + C_{12} \right)$

Let’s rewrite this equation using the conventional symbols for the constants.

$\displaystyle x = A \sin \left( t\sqrt{\frac{k}{m}} - \varphi \right)$

There it is. The basic sine wave equation is readily apparent, with the constants representing the amplitude and the phase shift, and the radical representing the angular frequency. The fourth possible constant (a “+ D”) is absent, and I wondered if I had accidentally dropped a constant somewhere — since that entire family of curves should be solutions. That is, x should be able to oscillate around 1, -5, whatever. But the reason for its absence is our starting point of Hooke’s law, F=-kx, which assumes that the equilibrium point is x=0. We would have had to start with F=-k(xxeq) to have had that last constant.

It’s pretty neat to see the equation emerge from just four basic equations: Hooke’s law F = –kx, Newton’s second law F=ma, and the definitions of velocity and acceleration v=x′ and a=v′!

# Snail Problem, Part 5

I hope no one is tired of snails yet. In the last several posts, I tackled the following problem: a snail crawls along a twig of length L0 at 1 cm/d, and the twig grows (along its entire length) at 2 cm/d. Will the snail reach the end, and when? In parts 1, 2, and 3, I solved this problem, and then in part 4 I turned to the general case, leaving the snail’s velocity and twig growth velocity unspecified.

We found that the snail’s position, x, can be given by the following equation, where g is the growth rate of the twig, s is the (inherent) velocity of the snail, and t is the time:

$\displaystyle x = \frac{s}{g}(gt+L_0)\ln\frac{gt+L_0}{L_0}.$

We want to find the time t where the snail reaches the end of the twig. Since the length of the twig is gt+L0, we set x equal to this.

$\displaystyle x = gt+L_0$

Substituting our equation for x gives

$\displaystyle \frac{s}{g}(gt+L_0)\ln\frac{gt+L_0}{L_0}=gt+L_0$.

To solve this for t, we start by dividing both sides by gt+L0 (note that this is the twig length, and is only applicable for nonzero twig lengths).

$\displaystyle \frac{s}{g}\ln\frac{gt+L_0}{L_0}=1$

We can multiply both sides by g/s (note that the snail velocity cannot be zero).

$\displaystyle \ln\frac{gt+L_0}{L_0}=\frac{g}{s}$

Exponentiating both sides gives

$\displaystyle e^{\ln\frac{gt+L_0}{L_0}}=e^{\frac{g}{s}}$

$\displaystyle \frac{gt+L_0}{L_0}=e^{\frac{g}{s}}$

$\displaystyle gt+L_0 = L_0 e^{\frac{g}{s}}$

$\displaystyle gt = L_0 e^{\frac{g}{s}}-L_0$

$\displaystyle gt = L_0(e^{\frac{g}{s}}-1)$

And finally, the time is

$\displaystyle t_f = \frac{L_0}{g}\left(e^{g/s}-1\right) \blacksquare$

So the time it takes is the initial length of the twig divided by its growth rate (which incidentally is the amount of time it would have taken to grow to its start position) times the quantity e to the power of the ratio of the twig and snail rates minus 1.

Interestingly, this will be finite no matter how small s is. If the snail is extremely slow, then the exponent g/s will be very large, but as long as s is positive, the snail will eventually reach the end of the twig — no matter how fast the twig grows!

Let’s also see how long the twig will be when the snail reaches the end. Recall that the length of the twig is

$\displaystyle L(t) = gt+L_0$

We can substitute our time above in for t:

$\displaystyle L(t_f) = g\left[\frac{L_0}{g}\left(e^{g/s}-1\right)\right] +L_0$

$\displaystyle L(t_f) = L_0\left(e^{g/s}-1\right)+L_0$

$\displaystyle L(t_f) = L_0\left[(e^{g/s}-1)+1\right]$

And finally, the length is

$\displaystyle L(t_f) = L_0 e^{g/s}\, \blacksquare$