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] |

*From*: John Denker <jsd@av8n.com>*Date*: Wed, 9 Dec 2020 16:43:51 -0700

Since people seem to be interested in the topic, I'll say a few more words.

1) As various people (including me) have said, the best approach is to

formulate the problem as a function of time: h(t), h'(t), and h''(t).

This is the conventional approach, and it has lots of advantages:

closer to the real physics, symplectic solver, et cetera.

Pedagogical point: It's always better to emphasize the good answer,

and explain why it is good, especially in the introductory course.

There will be plenty of time later for experts to delve into why the

bad answers are bad.

2) Having said that, let's put on our wizard hats and look into the

not-recommended version. It can be salvaged, and the salvage process

sheds considerable light on the root cause of the problem. Here's the

recipe:

The simplified equation of motion is:

h' = √(2 g h) [1]

The RHS is clearly a function of h. Arguably it is indirectly, implicitly

a function of t, since h itself depends on t, but that is useless to us

for present purposes, since we don't yet know the functional form of h(t).

The whole point of the exercise is to *find* h(t).

So there is nothing to lose and much to gain by thinking of it as a function

of h. We could even go so far as to write:

dt/dh = 1/√(2 g h) [2]

and then solve for t(h). In this case there is a convenient solution using

calculus (rather than numerical methods):

∫ dt = ∫ 1/√(2 g h) dh [3]

but let's pretend we didn't notice that.

Obviously integrating equation [2] or [3] numerically is a mess because the

RHS is infinite at h=0. However, we can salvage the situation as follows:

*Tabulate* dt/dh as a function of h, over the range of interest. The graph

is here:

https://www.av8n.com/physics/img48/dt-dh-integrator.png

Our mission, should we decide to accept it, is to find the area under the curve.

(Or in this case, the shaded area to the left of the curve. I plotted h vertically.)

To find the area between h=1 and h=2 we could take the value of dt/dh at h=1

and integrate downward, or we could take the value at dt/dh at h=2 and integrate

upward, or (!) take the value at the midpoint of the interval (h=1.5), which is

a superior estimate of the integrand.

When finding the area between h=0 and h=1 we have only two of those three

options. The integrand is infinite at h=0. However the values at h=1 and

(especially!) h=0.5 are still viable estimates. It's an improper integral

because the integrand diverges, but it is integrable.

End of recipe.

Note that this trick of using an improved estimate of the integrand does not work

if you think of everything as a function of t. That would require knowing h(t)

or at least h'(t), which we don't know. Remember, finding h(t) is the whole point

of the exercise.

This recipe is not something I cooked up special to handle this case. It is

standard good practice to make a graph of what you know. This is mentioned

explicitly in Apostol's _Calculus_ book and probably also in Acton. In a

higher-dimensional space, it turns into a /flow field/ with little vectors

at each point showing which way the derivative goes.

Make a habit of doing this. Encourage students to do this. The harder the

problem, the more useful this is.

It's the same thing with electronics: Draw the circuit diagram. It never

ceases to amaze me when beginners try to figure things out in their head in

situations where world-class experts would draw the diagram.

===

Related point: If you think of it in terms of h(t) and h'(t) and time-stepping

starting from t=0, it is clearly a /stiff/ differential equation, because there

is no possible time-step during which h' does not change significantly during

that interval. It changes by a factor of infinity during that first step.

===

Additional point, affecting both the dh/dt and dt/dh versions: Whenever you

write down a square root, ask yourself why it's not the plus-or-minus square

root. In this case, the energy principle that was used to derive the equation

of motion does *not* distinguish between plus and minus velocity. And at the

critical point, both branches coincide at zero. You know that you are only

interested in positive velocities, but the computer can't read your mind.

Writing the positive square root gives away the answer at all h>0 but *not*

at h=0.

Rule of thumb: Beware the proximity of other solutions. Even if you know you

are not interested in them, the solver doesn't know that (unless you intervene

using very sophisticated methods). This is a fairly common way to stumble

into a stiff differential equation.

In the present case, starting at zero, the solver would be equally happy to

integrate backwards or forwards in time, at least for the first step. This

symmetry is one more way of explaining why it gets stuck.

**Follow-Ups**:**Re: [Phys-L] Mathematica question; stiff differential equation***From:*John Denker <jsd@av8n.com>

**References**:**[Phys-L] simple Mathematica question (one more typo fixed)***From:*Carl Mungan <mungan@usna.edu>

**Re: [Phys-L] simple Mathematica question (one more typo fixed)***From:*Francois Primeau <fprimeau@gmail.com>

**Re: [Phys-L] Mathematica question; stiff differential equation***From:*John Denker <jsd@av8n.com>

**Re: [Phys-L] Mathematica question; stiff differential equation***From:*Carl Mungan <mungan@usna.edu>

- Prev by Date:
**Re: [Phys-L] Mathematica question; stiff differential equation** - Next by Date:
**[Phys-L] A Day in the Life** - Previous by thread:
**Re: [Phys-L] Mathematica question; stiff differential equation** - Next by thread:
**Re: [Phys-L] Mathematica question; stiff differential equation** - Index(es):