Previous ToC Up Next

3. Second-Order Differential Equations

3.1. Formulating the Problem

Alice: Now that we know how to solve a first-order differential equation, we can extend our methods immediately to treat the case of a general second-order differential equation

Here the force per unit mass exerted on a particle depends explicitly on both the position and the velocity of that particle. An example of such a force is the motion of a mass point under the influence of friction. Indeed, our physical interpretation of the case of a first-order differential equation followed from the above form in the limit of infinitely strong friction, where the velocity-dependent term dominated completely. Another example would be the force that is felt by an electron moving in an electromagnetic field (where would have be interpreted as a three-dimensional vector, a point we will come back to later in this chapter).

Bob: Wouldn't it be simpler to restrict ourselves immediately to the type of equations we are dealing with in stellar dynamics, without any velocity dependence in the forces?

Alice: I would prefer to hold off just a bit, because the simplest way to treat second-order equations is by following the same recipe as we did before. In that case, velocity dependence does not pose any problems.

The second-order equation above can be rewritten as a system of two first-order differential equations:

In principle we can look at this as a single first-order differential equation for a two-component vector.

Bob: Ah, that is a nice short-cut. Let's do that, and then we should be able to use all the results from the previous chapters immediately!

3.2. Vector Notation

Alice: We'll think have to put in some thought, since not all scalar equations generalize in an obvious way to the vectorial case, but I agree that it would probably be a good guide line.

Okay, to introduce vector notation, let us define:

and

We can then write our second-order differential equation as a single first-order equation in terms of vectors:

or simply:

When written out, this leads to the two equations implicit in:

Bob: In our case of interest, classical mechanics, we will drop the velocity dependence, and study instead the simpler equation:

or, equivalently

Alice: As we will see in the next chapter, we can exploit the simple form of these two equations to get an extra order of accuracy, seemingly for free, using what is called a partitioned Runge-Kutta algorithm. However, in order to put that clever trick in a clear context, in the current section we will be less clever. Rather, we will simply apply the same treatment that we have developed in the previous chapters.

Bob: Fine, but let's not linger too long! And please, let us drop the velocity dependence in the forces. Life is complicated enough as it is.

3.3. One Force Evaluation per Step

Alice: Okay, okay, we'll work with then. Now, everything we have done so far carries over directly to the case of two coupled differential equations. To show this, let us repeat the same derivation, but now in vector form. At the start of a time step, we evaluate the right-hand side of the differential equation at :

This leads to the following dimensionally correct expression:

Combining the last two equations, we have

We can compare this expression with the Taylor series:

where we use the shorthand notations

and

to avoid introducing yet another set of new symbols for the various derivatives of at time 0.

In order to compute

we will need to determine the time derivative of .

3.4. Not So Fast

Bob: presumably that is simply

in analogy to what we did in the first chapter, where we used .

Alice: Not so fast! You have not specified what you mean with that notation. The left-hand side is a vector, while the right-hand side suggests the product of two vectors. What does it mean?

Bob: Hmm. I hadn't thought about that. Good question.

Alice: If it would be an inner product, the left hand side should be a scalar. If, however, it is a tensor product, the left hand side should be a tensor. In neither case does it produce a vector.

Bob: Again, good point. Wel, in case of doubt, write it out! What does it look like in components?

Alice: Let us check. The most intuitive approach would be to start with a small variation in , which can be expressed as

This implies:

Bob: Good! So we do get a vector, after all.

Alice: Yes. And even if you would have allowed me to retain a velocity dependence in the force, we would still have wound up with a vector, but in this case we would have:

where now and . This would give us

Bob: Nice to know, but no thanks, let's stick with position-dependent forces only.

Alice: If you insist. At least the notation above will point the way for further generalizations, whenever we want to go that route.

Bob: Thanks!

3.5. Forward Euler in Vector Form

Alice: Let us introduce the symbol for the right-hand side of Eq. (106), at time 0:

We can then write Eq. (100) as

We demand that Eqs. (99) and Eq. (110) should be equal. The constant term matches trivially, and the first condition arises from the term linear in :

hence

We have no free parameter left, so this leads us to the only possible explicit first-order integration scheme

which is the forward-Euler algorithm, now in vector form. If we write this out in components, we get:

Let us define

where we will use the same symbol for the vector and the last componenent, again to avoid introducing yet more new letters. With this notation, we can write Eq. (111) in a more traditional form as:

We thus recover the forward-Euler algorithm. As before, this scheme is only first-order accurate, since it cannot possibly reproduce the term in Eq. (110).

Bob: Quite a bit of work to regain Forward Euler!

3.6. Two Force Evaluations per Step

Alice: Yes, and we might have guessed the result, but when we move on to using two force evaluations per step, things will undoubtedly get messier.

As before, after a first evaluation of the right-hand side of the differential equation, we can perform a preliminary integration in time, after which we can evaluate the right-hand side again, at a new position:

We can now use a more general expression for the new position, in which we rely on the preliminary information that has been gathered in the two force evoluations. Since each force evaluation has the dimension of a time derivative of the position, we have to multiply each one with a single power of . The coefficient for each term is, as yet, arbitrary, so let us parametrize them as follows:

We can combine these equations, and write them as

We want to determine

where

We have already seen that

and we can also write

The two results that are encoded here, and , can now be plugged back into the definition of

We finally get

where we use again the notation

Bob: A useful function, obviously.

3.7. Putting Everything Together

Alice: To sum up: developing in a Taylor series around gives:

Using this result in Eq. (119), we find for the new position, at the end of our time step:

We can now compare this expression with the Taylor series expansion of the true orbit, as we did in the previous section:

We can write this, as we saw in Eq. (110), in terms of as follows:

Starting with terms to first order in in Eqs. (128) and (130), we have to insist that

which leads to the condition

To second order in , we would like to satisfy:

which can be done through the condition

We conclude that when we allow two force evaluations, we again are left with one free parameter .

Bob: So, now we're done, and we can move on!

3.8. Summary

Alice: Yes, but let us put all our results on the table first.

To summarize, we can write:

Notice that this is exactly the vector generalization of the expressions we found in the case of a first-order differential equation:

Bob: I'd like to see our vector expressions written out in components. That way it will be easier to make contact with the previous chapters.

Alice: My pleasure! Starting with

we find

which leads to the expression for :

We can write this in components as

3.9. Two Examples

Bob: I find this a lot more understandable that the vector notation. And for practical application, let's look at a couple special case. For we find

and for we find

Alice: Ah, that is interesting! This is almost our leapfrog scheme.

Bob: How so?

Alice: The expression for is exactly the same as what we have used when we implemented the leapfrog.

Bob: And so is the expression for . The whole idea of leapfrogging is to advance the velocity with an acceleration that is the exact average of the force calculation at the beginning and at the end of the step.

Alice: Ah, the word exact is important here! While it is true that is the force evaluation at the beginning of the step, the exact force calculation at the end of the step would be . However, in our scheme above, , and the last term is missing.

Bob: Tricky! So we are dealing with an almost-leapfrog scheme, where the last force calculation is based on a predicted value for the new position, instead of the corrected value.

Alice: Yes, that's a good way of putting it. And all this can serve as an invitation to go beyond the straightforward generalizations of the Runge-Kutta schemes for first-order differential equations. It is time to look at more imaginative schemes, that treat position and velocity in different ways!
Previous ToC Up Next