Second order differential equations

Introduction

This week, we will extend our discussion of solving differential equations to include second order differential equations -- those involving a second derivative of one of the variables.

Netwon's second law, F=ma, is naturally a second order differential equation as the acceleration is a second derivative of the position, x. Last week, we solved Newton's Second Law for the special case where the force was either a constant or a function of the velocity. In that case, we could solve Newton's second law as a first order differential equation by writing everything as a function of the velocity, v.

This week, we will move to the more general situation where the force can be a function of the position and velocity as well as allowing for a free function of time.

Now that we are working with more complicated equations, our basic algorithm does not work as well. This is a common problem in programming and highlights the fact that you should not implicitly trust any answer produced by a computer, you should always think about whether the answer is sensible and where errors may have occurred due to the way the program works.

Syntax Summary

Nothing this week.

Worksheet

The Harmonic Oscillator

The example that we will solve is the simple harmonic oscillator (for example, a mass on a spring). The basic equation is $$F = -kx$$

Using Newton's second law, this can be written as

$$ ma = -kx$$

so

$$d^2x/dt^2 = dv/dt = - (k/m) x $$ In order to make use of the Euler method that we learned last week, we can re-write this as two coupled first order differential equations: $$dv/dt = -(k/m) x\\ dx/dt = v$$

Euler's Method

We can then solve each of these using the Euler method discussed last week by writing

$$v(t + dt) = v(t) - (k/m) x(t)\;dt\\ x(t + dt) = x(t) + v(t)\;dt $$

Of course, we still need to specify the initial position and velocity. In this example, we take x(0) = 5m and v(0) = 0 m/s.

Here is the code to do it:

# set constants
k = 10
m = 1
t_max = 10.0
#Set the number of iterations to be used in the for loop
no_of_iterations=1000

# set time step so that the loop will always iterate until t=t_max seconds 
dt = t_max/no_of_iterations


# make arrays to store data
t = zeros(no_of_iterations)
x = zeros(no_of_iterations)
v = zeros(no_of_iterations)

# set initial conditions
t[0] = 0
x[0] = 5
v[0] = 0

# evolve
for i in arange(1,no_of_iterations):
    t[i] = dt * i
    v[i] = v[i-1] - dt *k/m*x[i-1]
    x[i] = x[i-1] + dt * v[i-1]

Don't forget to add comments and breaks in your code to make it easy to read and edit. This program will get quite long so it is important to keep it neatly structured.

While this is the most straightforward method of solving a second order differential equation, it turns out that it is not very accurate. The first exercise involves solving this equation and examining the accuracy. It is, of course, possible to get reasonable results but it requires a very small time step. In general it is preferable to use a more sophisticated numerical method instead.

Energy Conservation in Euler's Method

In a perfect harmonic oscillator with no friction, the total energy is conserved. The total energy is $$\frac{1}{2} mv^2+ \frac{1}{2} kx^2$$ As the system evolves the energy is converted from kinetic to potential and back again.

Ideally, it should also be (approximately) conserved in a numerical scheme. However, for Euler's method, this is not the case. The numerical evolution of the system with Euler's method actually leads to an energy that increases exponentially! This means that the system is not stable and, even with a small time step, the solution will blow up if we run for a long time.

To see that the energy in Euler's method increases, we can use the evolution equations above to see what happens:

• Suppose in Euler's method we have a position x[i-1], velocity v[i-1] and total energy E[i-1] at time step i-1.
• We want to calculate the energy at time step i, E[i] = 0.5 mv[i]² + 0.5 kx[i]².
• We then substitute x[i] = x[i-1] + dt v[i-1] and v[i] = v[i-1] - dt *k/m*x[i-1] into the equation for E[i].
• After a bit of algebra, we find that E[i] = E[i-1] (1 + dt² (k/m) ).
• So, the energy at each time step changes by an amount proportional to dt². What's worse is that the energy will always increase -- the error term is strictly positive.

Over one period, $$T=2π (m/k)^{1/2}$$ the number of time steps will be equal to T/dt and at each step the energy increase is proportional to dt². The total energy increase due to the finite time step used in Euler's method, over one cycle, will be proportional to dt. This is not physical, it is due to the numerical method we have used to evolve the equations.

The Euler-Cromer Method

There is a straightforward modification to the Euler method which greatly improves the behaviour of the numerical solutions. In this, the Euler-Cromer method, we solve the equations

$$v(t + dt) = v(t) - (k/m) x(t)\;dt\\ x(t + dt) = x(t) + v(t + dt)\;dt $$

Note that we now use the velocity at time t + dt when evaluating the position. Inside the loop in the code, we evaluate x[i] as

x[i] = x[i-1] + dt * v[i]

It's a very simple change in the code, but it makes a significant difference to the results. To understand why this improves results, we can look again at energy conservation.

Energy Conservation in Euler-Cromer Method

As before, we can use the evolution equations to see how the energy changes with time:

• Suppose in the Euler-Cromer method we have a position x[i-1], velocity v[i-1] and total energy E[i-1] at time step i-1.
• We want to calculate the energy at time step i, E[i] = 0.5 mv[i]² + 0.5 kx[i]².
• We then substitute x[i] = x[i-1] + dt * v[i] and v[i] = v[i-1] - dt *k/m*x[i-1] to replace x[i] and v[i] in the expression for total energy.
• But, the x[i] expression contains a v[i], so we have to substitute again to get rid of it in favour of x[i-1] and v[i-1].
• After a few lines of algebra, we get E[i] = E[i-1] + dt² k/m ( mv[i-1]² - kx[i-1]²) + dt³( ...)
• As before the energy changes by a term proportional to dt² at each time step. But the term isn't necessarily positive. In fact, for a harmonic oscillator, the dt² term will oscillate over a cycle and give a zero total contribution.

It turns out that over a cycle, the energy is conserved up to dt³. This makes Euler-Cromer significantly more accurate than Euler's method.

Damping

In real life, no harmonic oscillator will keep going for ever. There will be a damping force which slows down the system and causes it to lose energy. Typically, this is modelled as an additional force which is proportional to the velocity (and opposes the motion):

$$F = - kx - cv$$

In this case, we can re-write the equations of motion as $$ v(t + dt) = v(t) - dt (k/m) x(t) - dt (c/m) v(t)\\ x(t + dt) = x(t) + dt v(t + dt) $$

The size of the damping term will have an effect on the motion of the oscillator. In particular:

• If c = 2 (k·m)^1/2, the system is said to be "critically damped". In this case, the oscillator will return to equilibrium as quickly as possible without oscillating.
• If c > 2 (k·m)^1/2, the system is said to be "overdamped". It returns to equilibrium without oscillating, but more slowly than a critically damped system.
• If c < 2 (k·m)^1/2, the system is said to be "underdamped". It performs oscillations with decreasing amplitude which finally reaches zero.

Driving

Finally, we want to add a driving force to the harmonic oscillator. Now the force becomes

$$F = - kx - cv + F_0(t)$$

We can now write the equations of motion as

$$ v(t + dt) = v(t) - dt\; (k/m)\; x(t) - dt\; (c/m)\; v(t) + dt\; F_0(t)/m\\ x(t + dt) = x(t) + dt\; v(t + dt) $$

The most common form of driving force is a periodic force:

$$F_0(t) = A\; sin(ω_0t)$$

Depending upon the frequency ω of the driving force, you can get quite different responses of the oscillator. In particular, driving at a frequency close to the natural frequency of the oscillator can lead to a large increase in the amplitude of the oscillations, a phenomenon known as resonance.

Excercises