Hi there,

numerics are usefull especially when nonlinear problems are considered. Here, I want to show how to calculate the trajectory for some simple one-dimensional oscillators.

You can imagine such an oscillator as a little mass attached to a spring. There is a point of rest, where all forces vanish. When the oscillator is displaced, a force appears originating from the spring.

When energy is put in the system, the mass oscillates. Depending on the form of the potential (which basically is the integral of the force with respect to the displacement) we can distinguish harmonic and anharmonic oscillators.

From textbooks (or wikipedia) we know the following items, which we want to check using numerics:

- for the harmonic (linear) oscillator, the trajectory is a pure sine wave of some frequency
- for the anharmonic (nonlinear) oscillator, more frequency components appear:
- there are only odd components for the symmetric nonlinear potential
- there are odd and even components for the asymmetric nonlinear potential

1. The linear (harmonic) oscillator

This is a well known textbook example. The potential energy of the harmonic oscillator is simply a parabula:

Epot = 1/2 k x^2

The resulting force is the negative derivative of the potential:

F = - dE/ds = - k x

This is Hooke's law which is a simple model for springs.

It is also the reason why the harmonic oscillator is called the linear oscillator: the force is directly proportional to the displacement (linear dependency).

Epot = 1/2 k x^2

harmonic potential. The potential energy grows quadratically with the displacement. |

F = - dE/ds = - k x

This is Hooke's law which is a simple model for springs.

Force acting in the harmonic potential |

The equation of motion is can be derived from the force

d^2 x /dt^2 = F/m = -k/m x

where x is the displacement and k the spring constant. We can easily give the solution for the equation:

x(t) = A1 cos( k^2/m^2 t + A2)

where A1 and A2 are two constants that have to be adapted to the initial values of the problem. The numerical solution using ODEs can be done following the steps described in one earlier post. As the linear oscillator is a special case of the more general problem i'll give the code later in this text.

After integration, we get a trajectory looking as the following (k=1, initial position =0, initalspeed = 2)

This looks like a nice, pure sine wave (as expected). To get more insight, we use the fourier fransform

trajectory of the harmonic oscillator |

of the trajectory x(t). You have to be careful using results from octaves ode solvers, because in most cases the output points are not equally spaced in time. Using the interpolation function interp1 can help to linearize the data.

The fourier spectrum looks like the following:

As expected, we only see one single frequency component, which means we indeed have a harmonic oscillation here.

The fourier spectrum looks like the following:

fourier components for the harmonic oscillator: there is only one frequency present |

2. The anharmonic / nonlinear oscillator

Considering the example with the spring, Hooke's law is valid for small displacements. When the displacement gets bigger, additional terms can appear. We want to consider two cases: a symmetric and an antisymmetric nonlinear potential. We give the potential energy:

Epot,sym = 1/2 k x^2 + 1/3 alpha |x|^3

Epot,asym = 1/2 k x^2 + 1/3 alpha x^3

here alpha is a nonlinearity constant. The curves of the potential are shown below:

Again, the forces occuring are the derivatives of the potential with respect to the displacement:

Considering the example with the spring, Hooke's law is valid for small displacements. When the displacement gets bigger, additional terms can appear. We want to consider two cases: a symmetric and an antisymmetric nonlinear potential. We give the potential energy:

Epot,sym = 1/2 k x^2 + 1/3 alpha |x|^3

Epot,asym = 1/2 k x^2 + 1/3 alpha x^3

here alpha is a nonlinearity constant. The curves of the potential are shown below:

Two aharmonic potentials, compared with the harmonic one. They differ in their symmetry |

Fsym = - dEpot,sym/ds = - (k x + alpha |x| x)

Fasym= - dEpot,asym/ds = - (k x + alpha x^2)

We see additional dependencies between the force and the displacement which are of the order of x^2. This is why these oscillators are called nonlinear oscillators. The difference to the linear case becomes clear when we look at the force plot:

Forces derived from the potentials above. Note that the asymmetric force is equal to the symmetric force for positive displacements, but differs for negative displacements. |

For small displacements, the forces are almost identical. They begin to differ considerably, when the displacement grows. For the symmetric potential, the absolute value of the force is equal for -x and x. For the asymmetric potential, the force for negative displacements is smaller than for positive values.

The equations of movements can be written as

(symmetric)

d^2 x /dt^2 = F/m = - (k/m x + alpha/m |x| x)

(asymmetric)

d^2 x /dt^2 = F/m = - (k/m x + alpha/m x^2)

Again, we can easily solve this with the odepkg (see code sample below). I used the same input energy as for the harmonic example resulting in the following trajectories:

The two new trajectories basically also show sine osciallations. As the force for the symmetric potential is stronger for big displacements, the maximum amplitude is smaller. This also results in a higher frequency for the oscillation. For the asymmetric potential, the oscillation is shifted slightly towards negative values in the mean.

We can get the fourier components of this oscillations when analyzing the time traces x(t) of the signals.

We have a closer look on the symmetric harmonic potential:

Besides the fact, that the base frequency is higher than for the harmonic oscillator, we see an additional frequency component! Its value is three times the one of the base frequency (third harmonic). Also, energy is converted into the fifth harmonic.

Trajectories for the nonlinear oscilators (two kinds), compared with the one of the harmonic oscillator |

The two new trajectories basically also show sine osciallations. As the force for the symmetric potential is stronger for big displacements, the maximum amplitude is smaller. This also results in a higher frequency for the oscillation. For the asymmetric potential, the oscillation is shifted slightly towards negative values in the mean.

We can get the fourier components of this oscillations when analyzing the time traces x(t) of the signals.

We have a closer look on the symmetric harmonic potential:

spectrum for the symmetric nonlinear oscillator. Only odd multibles of the base frequency are generated |

A similar behavior can be seen for the asymmetric potential:

Here we have generated additional frequencies at two, three and four times the base frequency. The conversion efficiency is higher than for the symmetric potential.

spectrum for the asymmetric nonlinear oscillator. Here, even and odd multibles of the base frequency are generated |

As a result, the simulations show the expected behavior. We can see that nonlinear oscillators are capable of generating multiples of the base frequency. This is extensively used in nonlinear optics, where light is sent through nonlinear crystals. The crystal structure is responsible for the potential, in most cases they have asymmetric potentials.

When e.g. light from a laser (at a wavelength of 1064 nm) of high intensity passes such a crystal, green light (532 nm) can be generated.

see

for further details.

3. Code to solve the linear and nonlinear oscillators using octave

clear all;

more off;

%

% right hand side of the ODE for the symmetric harmonic potential

%

function dx = rhs_sym(t, x, k, alpha,tend)

t/tend %just prints the progress

xp = x(1); % position

vx = x(2); % velocity

dxp = vx; % position change

dv = -(k * xp + alpha * xp * abs(xp) ); %velocity change

dx = [dxp, dv]; %give the derivative

endfunction

%

% right hand side of the ODE for the asymmetric harmonic potential

%

function dx = rhs_asym(t, x, k, alpha,tend)

t/tend

xp = x(1);

vx = x(2);

dxp = vx;

dv = -(k * xp + alpha * xp^2);

dx = [dxp, dv];

endfunction

% program parameters

k =1;

alpha =0.2;

initialDisplacement = 0;

initialVelocity = 2;

xini = [initialDisplacement, initialVelocity];

tend = 15; %integrate up to that time

% set integration options

options = odeset( 'InitialStep',tend/1e3,

'MaxStep',tend/1e3,

'RelTol',1e-6,

'AbsTol',1e-6,

'NormControl','on')

% linear case, using alpha = 0

[T1, Xlin] = ode45( @rhs_sym, [0, tend], xini, options, k, 0, tend);

%NL case, symetric potential:

[T2, XnlS] = ode45( @rhs_sym, [0, tend], xini, options, k, alpha,\

tend);

%NL case, symetric potential:

[T3, XnlA] = ode45( @rhs_asym, [0, tend], xini, options, k, alpha,\ tend);

% plot the result

plot (T1, Xlin(:,1),'color',[0.5,0.5,0.5],

T2, XnlS(:,1),'.','color',[0,0,1],

T3, XnlA(:,1),'.','color',[1,0,0]);

xlabel 'time'

ylabel 'displacement'

Hint: for the analysis of the frequency components shown above, the integration time has to be set to a higher value and the odeoptions have to be changed to higher accuracy...

%

% right hand side of the ODE for the symmetric harmonic potential

%

function dx = rhs_sym(t, x, k, alpha,tend)

t/tend %just prints the progress

xp = x(1); % position

vx = x(2); % velocity

dxp = vx; % position change

dv = -(k * xp + alpha * xp * abs(xp) ); %velocity change

dx = [dxp, dv]; %give the derivative

endfunction

%

% right hand side of the ODE for the asymmetric harmonic potential

%

function dx = rhs_asym(t, x, k, alpha,tend)

t/tend

xp = x(1);

vx = x(2);

dxp = vx;

dv = -(k * xp + alpha * xp^2);

dx = [dxp, dv];

endfunction

% program parameters

k =1;

alpha =0.2;

initialDisplacement = 0;

initialVelocity = 2;

xini = [initialDisplacement, initialVelocity];

tend = 15; %integrate up to that time

% set integration options

options = odeset( 'InitialStep',tend/1e3,

'MaxStep',tend/1e3,

'RelTol',1e-6,

'AbsTol',1e-6,

'NormControl','on')

% linear case, using alpha = 0

[T1, Xlin] = ode45( @rhs_sym, [0, tend], xini, options, k, 0, tend);

%NL case, symetric potential:

[T2, XnlS] = ode45( @rhs_sym, [0, tend], xini, options, k, alpha,\

tend);

%NL case, symetric potential:

[T3, XnlA] = ode45( @rhs_asym, [0, tend], xini, options, k, alpha,\ tend);

% plot the result

plot (T1, Xlin(:,1),'color',[0.5,0.5,0.5],

T2, XnlS(:,1),'.','color',[0,0,1],

T3, XnlA(:,1),'.','color',[1,0,0]);

xlabel 'time'

ylabel 'displacement'

Hint: for the analysis of the frequency components shown above, the integration time has to be set to a higher value and the odeoptions have to be changed to higher accuracy...

very nice, I was exactly looking for this! Thank you

ReplyDelete