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] simple Mathematica question (one more typo fixed)

In a mathematical sense there is no problem at all, because y(t)=100 is a perfectly valid solution to the energy-conservation differential equation.

This "unphysical" solution comes from taking Newton's 2nd Law for a single conservative force:
[$m\ddot{y} = -U'(y)$]

and multiplying both sides by the velocity in order to get a "first integral":
[$m\ddot{y}\dot{y} = -U'(y)\dot{y} \Longrightarrow \frac{d}{dt}\left[\frac{1}{2}m\dot{y}^2 + U(y)\right] = 0$]

In doing so, you introduce another possible (spurious) solution to the equation, in which [$\dot{y} = 0$] identically. Multiplying both sides of the N2 equation by zero effectively erases the dynamics.

- Craig

On 12/8/20 12:17 PM, Carl Mungan via Phys-l wrote:


Since people seem to be in a talkative mood today, maybe somebody can help me a with a simple (?) matter. It involves dropping a rock. I’m using Mathematica but you don’t have to be an expert in that to see the issue.

If we use conservation of mechanical energy, then we find the downward speed of the rock (choosing upward to be +y with the ground at y=0) is v_y = - Sqrt[2*g*(y0-y)].

Suppose we choose y0 = 100 m and g = 9.8 m/s/s and ask Mathematica to do a numerical solution from t = 0 to 4 s:

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

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

Can you guess what happens here? It’s a failure. The rock just sits at rest in the air at 100 m.

If I instead start the evaluation at 99 m height, so the rock has already fallen 1 m, then it works just fine:

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

But I want to start the rock at 100 m. How do I nudge the rock so it starts falling?

I must be missing something. My actual problem is like this, but more complicated. The above is just a test case for debugging purposes.

Who can tell me what I need to do? -Carl

ps: Converting to acceleration does work but I don’t see why I should have to do that:

s2 = NDSolve[{y''[t] == -9.8, y[0] == 100, y'[0] == 0}, y, {t, 0, 4}]

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

Carl E. Mungan, Professor of Physics 410-293-6680 (O) -3729 (F)
Naval Academy Stop 9c, 572C Holloway Rd, Annapolis MD 21402-1363 <><> <><>
Forum for Physics Educators<>

Craig C. Wiegert, Associate Dept. Head
Associate Professor of Physics Office: 215 Physics Bldg.
Department of Physics & Astronomy Phone: 706-542-4023
University of Georgia Fax: 706-542-2492
Athens, GA 30602-2451 Email:<>