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?
What ODE solver are you using? I would use an adaptive Runge-Kutta-Fehlberg ODE solver with a tolerance of at least 1.0E-7 for non-linear or chaotic systems. Look at the ODE solution of the Double Pendulum.
Placing the computation of F_D_x in the Fixed relations is wrong. The effect is the same as using the Euler method. You must place this computation in the "Prelim Code" page of the ODE solver page.
This allows the solver to correctly compute the force in the intermediate internal steps that any solver (except Euler) takes.
If you have problems to do this, please attach the example and I can fix it for you.
