No damping but the solution to simple harmonic oscillator damps?

[member 428835]

I posted yesterday but figured it out; however, a different issue I just detected with the same code arose: namely, why does the solution damp here for an undamped simple harmonic oscillator? I know the exact solution is $\cos (5\sqrt 2 t)$.

global delta alpha beta gamma OMEG
delta      = 0;     % DAMPING
alpha     = 50;   % STIFFNESS
beta       = 0;     % RESTORATION NONLINEARITY
OMEG    = 1.2;  % DRIVING FORCE ANGULAR FREQ
gamma  = 0;     % FORCING

Fs = 100;       % Sampling frequency                  
T  = 1/Fs;      % Sampling period     
L  = 10000;     % Length of signal

[t, x]  = ode45(@duffing,(0:L-1)*T,[0 1]);

plot(t,x(:,2))

and the function ode45 calls is

function xdot = duffing(t,x)
global delta alpha beta gamma OMEG
xdot(1) = -delta*x(1) - alpha*x(2) - beta*x(2)^3 + gamma*cos(OMEG*t);
xdot(2) = x(1);
xdot    = xdot';
end

 

[jedishrfu]

Damping can occur when error is introduced somehow. The error in this case would be subtracting from the system energy causing it to appear as damping.

You could try another ODE solver beside ODE45 to see if that fixes the problem.

When simulating stuff, error can introduce more energy or take away energy depending on the algorithm chosen. As an example, Euler methods are good for systems that are non-periodic such as tracking a particle trajectory and ODE45 is good for periodic systems. In general, though ODE45 always seemed to be the best even though it could be slower in some cases.

You could also try a finer step.

Here's how to choose an ode solver:

https://www.mathworks.com/help/matlab/math/choose-an-ode-solver.html
Here's the MATLAB description of the ode45 solver:

https://www.mathworks.com/help/matlab/ref/ode45.html
which mentions it might not be good for very stiff systems.

As an aside, I found some Matlab notes on using ode45 that might be of interest:

https://www.12000.org/my_notes/matlab_ODE/

 

[jedishrfu]

When we ran simulations of planetary orbits in java code using one of the Euler variants, we would get either a planet crashing into the Sun or zipping off on a journey of adventure and off our screen due to the algorithm introducing small error that manifested as energy loss or energy gain in a closed system.

The trick was to find an algorithm that alternated adding some error and then taking it away and between those two limits produced an effective and accurate simulation over longer periods of time.

 

[Orodruin]

The trick was to find an algorithm that alternated adding some error and then taking it away and between those two limits produced an effective and accurate simulation over longer periods of time.

Another way would be to find an algorithm that respects the important symmetries of the system, ie, one that is explicitly designed such that energy is conserved.

 

[member 428835]

Damping can occur when error is introduced somehow. The error in this case would be subtracting from the system energy causing it to appear as damping.

Can you elaborate on this please?

 

[member 428835]

Another way would be to find an algorithm that respects the important symmetries of the system, ie, one that is explicitly designed such that energy is conserved.

Any ideas for such an algorithm? See, I'm using my developed code to solve the Duffing oscillator, so there are nonlinear components. I just thought I'd benchmark the code with an equation that admits an analytic solution first.

 

[DrClaude]

What are the initial conditions?

For "fun," I implemented the Duffing oscillator in C using RKF45 from GSL, which should be the same integrator as the one used by Matlab, and I get the correct solution for the undamped case.

 

[member 428835]

What are the initial conditions?

For "fun," I implemented the Duffing oscillator in C using RKF45 from GSL, which should be the same integrator as the one used by Matlab, and I get the correct solution for the undamped case.

Hmmmm this is odd. I use $x(0) = 1$ and $\dot x (0) = 0$.