No damping but the solution to simple harmonic oscillator damps?

  • #1
joshmccraney
Gold Member
1,957
102
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)##.

Code:
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

Code:
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
 

Answers and Replies

  • #2
12,097
5,778
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/
 
Last edited:
  • Like
Likes joshmccraney
  • #3
12,097
5,778
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.
 
  • #4
Orodruin
Staff Emeritus
Science Advisor
Homework Helper
Insights Author
Gold Member
16,829
6,650
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.
 
  • Like
Likes jedishrfu
  • #5
joshmccraney
Gold Member
1,957
102
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?
 
  • #6
joshmccraney
Gold Member
1,957
102
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.
 
  • #7
DrClaude
Mentor
7,433
3,704
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.
 
  • Like
Likes jedishrfu
  • #8
joshmccraney
Gold Member
1,957
102
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##.
 

Related Threads on No damping but the solution to simple harmonic oscillator damps?

Replies
2
Views
11K
Replies
6
Views
1K
  • Last Post
Replies
1
Views
2K
Replies
3
Views
6K
Replies
6
Views
18K
Top