Chronology Current Month Current Thread Current Date
[Year List] [Month List (current year)] [Date Index] [Thread Index] [Thread Prev] [Thread Next] [Date Prev] [Date Next]

Re: [Phys-l] accurate numerical solution of equations of motion



Because I am a huge fan of using spreadsheets for quick (and I mean QUICK) and (often not even so) dirty numerical solutions of DE's, this really caught my eye. So I put together a side-by-side Euler/ Denker solution of the Kepler problem and was genuinely surprised and disappointed to find that, while both seemed to perform equally well at short dt's, the Euler method did much better at longer dt's, i.e. far less precession of the orbital ellipse.

At first I thought I must have implemented the Denker method incorrectly, but except for the "getting started" step in which one has to revert to the standard Euler method to generate v(2), it was precisely as prescribed.

On the other hand, I had not used Denker's prescription of the Euler method, but rather the so-called Euler-Cromer method [see "Stable solutions using the Euler approximation," Alan Cromer, 49 (5), 455-9 (1981)] in which Denker's Eq 2 is replaced with

x(3) = x(2) + v(3) Δt

This method is so vastly superior to the standard Euler method and so simple to implement that I tend not even to remember that I'm doing it any more. Indeed, when I went back and put my Euler method solution into standard form it was clear that the Denker method was the hands-down winner.

So it seems to me from at least my Kepler solution, that the Denker method is almost as good as the Euler-Cromer method, but not as easy to implement. I haven't done the analysis to know whether this should be expected to be the case in general or not.

John Mallinckrodt
Cal Poly Pomona

On Oct 21, 2009, at 11:55 PM, John Denker wrote:

Hi --

Many times over the years I've needed to write a program
to find numerical solutions for an equation of motion such
as F=ma.

We can rewrite the equation of motion as

(d/dt) (d/dt) x = (1/m) F(x)

where x is a vector. The RHS is considered known and we
have to integrate the equation to find x(t).

Nowadays, it is popular to use the simplest possible method,
i.e. first-order Euler, as described below. To achieve any
reasonable accuracy with this method, it is necessary to use
a very small stepsize Δt along with a very large number of
steps. An accuracy of 1 ppm could easily require a million
iterations.

If you are writing in c++ and running on a modern desktop, a
million iterations can be done in less time than it takes to
talk about it, hence the popularity of the brute-force approach.

However ... I would like to talk about the old-school approach.
I reckon a few people on this list might be interested. Even
if you don't need it right now, it is worth remembering the
general idea, so that you can look up the details if/when you
ever have the need.

For example, suppose you are writing for some sort of embedded
system that has a thousand or a million times less CPU power
than your desktop. Or suppose you want to implement the
solver using a spreadsheet (instead of c++) and you don't want
the spreadsheet to be -- literally -- a million lines long. Or
suppose you are getting killed by accumulated roundoff errors,
in which case using a smaller stepsize and a larger number of
iterations will make things worse not better. My point is that
it is possible to use methods where each step is accurate to
second order in the stepsize, so you can achieve part-per-million
accuracy using only a thousand steps or so.

Equation [1] is a second order differential equation. Beware
that the term "second order" in this paragraph means something
completely different from the same term as used in the previous
paragraph. In this paragraph, it means that equation [1] involves
the second derivative. The simplest way to proceed is:
1) First we integrate the acceleration (a) to find the
velocity (v). This is only a first-order differential
equation.
2) Then, given (a) and (v), we integrate to find the
position (x).

Step (2) is easy. Whereas the first-order Euler method would set

x(3) = x(2) + v(2) Δt [2]

we instead use second-order Euler, that is,

x(3) = x(2) + v(2) Δt + 0.5 a(2) (Δt)^2 [3]

This is correct to second order in Δt, i.e. the errors will be
third order in Δt. So with a thousand steps we can achieve ppm
accuracy.

We also need an accurate way of doing step (1), i.e. integrating
the acceleration to find the velocity. We cannot use the same
trick as in step (2), because we don't know the second derivative
of velocity. In analogy to equation [2] we could write

v(3) = v(2) + a(2) Δt [4]

but there is no usable analogy to equation [3]. And equation [4]
itself is not accurate enough for present purposes.

Instead we use the scheme illustrated in the following figure:
http://www.av8n.com/physics/img48/stride-2.png

The key equation is

v(3) = v(1) + 2 a(2) Δt [5]

That is, rather than computing v(3) in terms of the most
recent known v -- i.e. v(2) -- we compute it in terms of the
slightly older v -- i.e. v(1).

For some non-experts, it goes against intuition to use older
data when newer data is available, but you can see from the
figure (lower right) that v(2) + a(2) Δt is significantly
different from v(3). In contrast, it is easy to show that
equation [5] (illustrated in the lower left corner of the
figure) is correct to second order in Δt ... even though
no explicit (Δt)^2 term appears in the equation. The errors
are of order (Δt)^3.

Homework: Do the math to verify these assertions.

Forsooth, for strictly circular motion, equation [5] is
good to all orders. That is to say, it is exact. For motion
that is approximately circular, such as for an elliptical
orbit with "small" eccentricity, the errors in equation [5]
are something "small" times (Δt)^3 ... which is very nice.

Beware: If you are using stride=2 to calculate v from a, it
is a Bad Idea to use stride=2 to calculate x from v. Do you
see why? If you don't believe me, try it sometime and see
what happens.

If you want to try out these techniques, I suggest setting up
a spreadsheet to solve the Kepler problem, i.e. F(x) ∝ 1/x•x.
This has the nice property that you can check the numerical
solution against the well-known exact analytical solution.
-- check that energy is conserved to ppm accuracy
-- check that angular momentum is conserved to ppm accuracy
-- check that the orbit closes on itself, as opposed to
spiraling out or spiraling in.
-- check that the numerical solution is in fact an ellipse
to ppm accuracy.

Another nice property is that it is easy to use the spreadsheet
to make a plot of the solution.

Solving the Kepler problem numerically won't tell you anything
you didn't already know ... except that it gives you confidence
in the method. You can then use the numerical method for other
problems that don't have a nice analytic solution.

Example:
http://www.av8n.com/fly/turns-on-pylons.htm

================

I know it is somewhat eccentric to use elegant old-school numerical
techniques on a spreadsheet. Many old-school numerical analysts
despise spreadsheets. But bear with me; I'm not as crazy as
you might think. There is a segment of the population that is
intimidated by c++ and is not going to learn it anytime soon ...
but is not intimidated by spreadsheets. The folks in this segment
are never going to read big fancy books on numerical methods ...
but you can teach them a couple of good ideas, like using second-
order Euler for part of the problem, and using stride=2 for the
other part. These ideas allow them to get really good results
from their spreadsheets.

Pedagogical remark: Almost all real-world physics jobs involve
a ton of computing. IMHO there should be more of an effort to
integrate computing into the physics curriculum at all levels.
I'm not sure _how_ to do this, given how overcrowded the schedule
already is ... but still it needs to be done somehow. Otherwise
the product will suffer.

==============================

There's a lot more I could say about this if anybody is interested,
but I'll stop here for now.
_______________________________________________
Forum for Physics Educators
Phys-l@carnot.physics.buffalo.edu
https://carnot.physics.buffalo.edu/mailman/listinfo/phys-l