APS Excellence in Physics Education Award
November 2019 Science SPORE Prize
November 2011 The Open Source Physics Project is supported by NSF DUE-0442581.

## Different behavior post and replies

Different behavior
Ivo Gienal
9 Posts

Normally I use the Newton 2 version in the ODE Evolution tab: dv/dt = F_Res/m, as I teach my students to use in solving problems. As an example of an anharmonic system I sometimes use an experiment with an air track, where the glider is attached to a spring vertically above the middle of the air track, similar to the quartic. ejs example from the OSP digital library. I introduce a variable F_D_x for the force in the x direction and in the Fixed relations I put: F_D_x = D*x*(1 - h/Math.hypot(x,h)), with D = spring constant and h = the length of the spring in equilibrium position. On the ODE tab I write: dv/dt = - F_D_x/m. With a moderate time step dt = 0.1 the resulting oscillating has an increasing amplitude instead of the expected constant one (anharmonic_oscillator_1.jar). When I put in the ODE Tab dv/dt = - D*x*(1 - h/Math.hypot(x,h)), the amplitude is constant (anharmonic_oscillator_2.jar). Can someone explain the different behavior?

Attached File: anharmonic_oscillator_1.zip

### Replies to Different behavior

Re: Different behavior -
Francisco Esquembre
228 Posts

Ivo,

ODE solvers in EJS (except Euler) evaluate the rate more than once, for different intermediate steps, inside the approximating algorithm to improve the accuracy of the (approximated) solution.

This means that you MUST pass all information about that intermediate state from the rate table. If you delay updating the F_D_x variable until the FIxed relations, the rate cannot be computed correctly for these intermediate states. Hence, you are basically applying Euler method, which results in the increasing amplitude you describe for this problem.

To do it correctly you can either:
- type the complete rate as -D*x*(1-h/Math.hypot(h,x)) in the table
- define a custom method:
double F_D_x (double x) { return D*x*(1-h/Math.hypot(h,x)); }
and type in the rate table: - F_D_x(x);
notice that I pass the x value from within the ODE. EJS knows how to pass the correct value of x for any intermediate computation of the rate.

There is no need to pass v, nor h, nor D, since :
- v is not involved in the computation (if it were, then you would need to pass it, too)
- D and h do not change in the ODE (they are not listed as part of the state and do not depend on x or v)

Paco

OSP Projects: Open Source Physics - EJS Modeling Tracker Physlet Physics Physlet Quantum Physics STP Book