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*: Francois Primeau <fprimeau@uci.edu>*Date*: Wed, 9 Dec 2020 07:29:30 -0800

The problem is with the DE itself and not the numerical method perse. The

uniqueness theorem for a differential equation dy/dt = F(t,y) is that both

F(t,y) and \partial F/\partial y be continuous functions in the domain of

interest. For the DE y' = sqrt(2*g*(100-y)) the partial derivative with

respect to y is not continuous at the point y = 100 and indeed the solution

is not unique at the point y = 100. (Thank you JD for clarifying the

distinction between derivatives with respect to t and y which was muddled

in my original response.)

The system X' = F(X) where X = [y;v] and F(X) = [v;-g], does satisfy the

condition for uniqueness. A sufficient condition is that F should have

bounded first partial derivatives with respect to the components of X. In

going from the original 2nd order equation (or equivalently from the system

of 2 first order equations) to the equation y' = sqrt(2*g*(100-y)), some

of the information needed to guarantee uniqueness was lost.

When I learned about existence uniqueness theorems in my ODE classes a

long time ago I remember finding it difficult to care. My (immature)

thinking was that I didn't need to worry about such things because I was

interested in real problems from nature, not artificial problems cooked up

to violate the conditions of the theorem. The present PHYS-L discussion has

given me an opportunity to reconsider my attitude. It also seems to me

that we have an opportunity with this example to construct a homework

problem that illustrates the condition of the theorem with a physical

problem that is quite compelling. I regularly teach a course on

mathematical modeling for Earth System Science students. I'm thinking I

should have a go at creating such a problem starting from the example we

have been discussing. If I'm successful I will call it the Wile E Coyote

problem!

- Francois

On Wed, Dec 9, 2020 at 5:28 AM Carl Mungan via Phys-l <

phys-l@mail.phys-l.org> wrote:

The responses have been helpful and interesting. I have heard of “stiff”

DEs before but never understood what the term meant. Now I see it’s in part

because it can mean many things, but overall it means a numerical

instability in a DE.

But there’s still one part I’d like some more insight about. Is the

problem with the DE itself (and the numerical method used to solve it) or

is the problem with certain points only (what you call fixed points)? It

seems to me this is two quite different questions.

I’ll illustrate by going back to my example. We’ll initially avoid the

fixed point by again starting our ball in motion at y=99. But this time

let’s reverse the velocity direction and throw the ball upward. So let’s

put this code into Mathematica:

s = NDSolve[{y'[t] == Sqrt[2 9.8 (100 - y[t])], y[0] == 99}, y, {t, 0, 4}]

Plot[Evaluate[y[t] /. s], {t, 0, 4}]

Can everyone guess what happens when I execute this?

The answer is it runs fine until it hits the peak. At that point,

Mathematica reports: “"

The plot runs from 99 up to 100 and then abruptly ends. It doesn’t

flatline at the peak. It just ends there.

In contrast, if we run:

s2 = NDSolve[{y''[t] == -9.8, y[0] == 99, y'[0] == Sqrt[2 9.8 (100 -

y[0])]}, y, {t, 0, 4}]

Plot[Evaluate[y[t] /. s2], {t, 0, 4}]

then there is no problem. You get exactly the expected up and down

quadratic.

So indeed, what you all are saying is correct. The problem is that in the

energy approach we have a term like v^2 in the KE. When we differentiate it

w.r.t. y, we get 2 v dv/dy and dv/dy = dv/dt dt/dy = a/v which is well

defined everywhere *except* at v=0. Nevertheless the correct answer is 2 v

a/v = 2 a everywhere including at v=0.

I have to ask once more: Is the problem in how the numerical integration

is performed, or is there really some information missing in the energy

equation. (In the present example, the missing information might be: Use

l’Hopital’s rule at any point that v=0, but I could be oversimplifying the

real issue.)

A possibly related issue is to think about a mass oscillating on a spring.

In the energy equation, I’ve lost information about the sign of the

velocity v. This is particularly problematic at v=0 because zero is

unsigned. But again I might be oversimplifying.

-Carl

On Dec 8, 2020, at 7:20 PM, John Denker via Phys-l <phys-l@mail.phys-l.org> wrote:

number

On 12/8/20 10:22 AM, Francois Primeau via Phys-l wrote:

The problem is that Mathematica knows mathematics but not physics.

Yes.

As stated, your problem is ill posed because you start it off at a

point at which the derivative does not exist.

Yes.

Perhaps there is a way in Mathematica to start off your problem atYes, that's one possibility.

y=100-epsilon for epsilon>0 and then take the limit as epsilon

approaches zero.

Earlier:

an easy way to avoid the fact that the derivative of the square rootThat's even better.

function does not exist at zero is to convert your equation into the

system y’=v and v’=g.

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

FP's answer is correct and concise. Perhaps I can fill in a couple of

details. For starters, let's be clear: Do not confuse the statement:

y'' = -g for all time (which is true)

with the statement:

y'' = -g for all y (which is not true).

Beware: y simply does not exist for y>100.

Similarly: y is not twice differentiable for y>=100.

In more detail:

h = 100 - y # by fiat; new convenient variable

h' = -y' # by differentiating

h' = √(2 g h) for all h>=0 # from conservation

h'' = ½ √(2 g h)¯¹ ⋅ 2 g h' for all h>0 (not including h=0)

= g for all h>0

Mathematica knows mathematics but not physics

Mathematica has no reason to compute h'', and would be unable to compute

it even if it wanted to. So (h,h') = (0,0) is a fixed point.

This is a huge problem, because h'' is important to the physics, and is

presumably deeply built into your intuition. Alas it is not built into

this equation.

This corresponds to Zeno's paradox in reverse. It takes an infinite

of steps to get away from the fixed point.differential

Here's another bit of jargon that may help: This is a /stiff/

equation. Pretty much worst-case stiff.that

https://en.wikipedia.org/wiki/Stiff_equation

Useful practical rule of thumb: Avoid(*) using differential equations

bump up against the boundaries of the domain where the variables aredefined.

very

(*) I was going to say "avoid like sin" or "avoid like the plague"

but we need to find better similes. It turns out people aren't

good at avoiding those things.which

convert your equation into the system y’=v and v’=g.

That would be my recommendation.

That is better math-wise as well as physics-wise.

Among other things, that sets you up to use a symplectic integrator,

will probably outperform anything else rather spectacularly.

https://www.av8n.com/physics/symplectic-integrator.htm

_______________________________________________

Forum for Physics Educators

Phys-l@mail.phys-l.org

https://www.phys-l.org/mailman/listinfo/phys-l

-----

Carl E. Mungan, Professor of Physics 410-293-6680 (O) -3729 (F)

Naval Academy Stop 9c, 572C Holloway Rd, Annapolis MD 21402-1363

mailto:mungan@usna.edu http://usna.edu/Users/physics/mungan/

_______________________________________________

Forum for Physics Educators

Phys-l@mail.phys-l.org

https://www.phys-l.org/mailman/listinfo/phys-l

--

Professor of Earth System Science

Rm 3224 Croul Hall

University of California, Irvine

**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:
**Re: [Phys-L] Mathematica question; stiff differential equation** - Previous by thread:
**Re: [Phys-L] Mathematica question; stiff differential equation** - Next by thread:
**Re: [Phys-L] Mathematica question; stiff differential equation** - Index(es):