Creative Commons License Copyright © Michael Richmond. This work is licensed under a Creative Commons License.

Using Euler's Method to solve Ordinary Differential Equations

See Sections 1.1, 25.1-25.2 of your textbook for more information on these methods.

A differential equation is one which expresses the change in one quantity in terms of others. Thus,

          y = x  
is not a diff eqn, whereas

         -- = 2x
is a diff eqn.

Diff eqns occur very frequently in all branches of physics, and so we must devise ways to deal with them. Here are several:

We'll discuss analytic solutions and Euler's method this week.

Analytic solutions

Sometimes, when you're lucky, the physical situation has an analytic solution. For example, consider the relationship between velocity and time, for the case of constant acceleration:

     --  =  a         =  -g

      v  =  at + v0   =  -gt + v0
Pretty simple! In this case, we can even find an exact relationship between velocity, position and time:

      --  =  v  =  -gt

               1     2
      y   = - --- g t    +  y0

If you can find an exact analytic solution, use it. That's always the best way to go, since you can figure out the quantities of interest (position, velocity) for any time, with absolutely no error.

Unfortunately, in many situations, the physical relationships between quantities is more complicated, such as

     dv                    2
     --  =  a  =  -g  + K v
In most of these cases, there is no exact solution -- or, if there is one, it might be very hard to figure out. Under these circumstances, we can apply numerical methods to understand the behavior of the system -- but we always fall back on analytic solutions to check the numerical results.

Euler's method

Let's consider the following real-life problem: Joe Geronimo jumps out of an airplane at an altitude of H = 100 meters. What are his altitude y(t) and velocity v(t) as functions of time? How long until he hits the ground?

We'll start by making a simplification: assume that there is no air resistance. We'll include air resistance and other complications in the future.

In this simplified situation,

     d y        dv
    ------   = ----  =  a  =  -g  =  -9.8 m/s^2
       2        dt

This can be solved analytically!

      v(t)  =  v0  -  g*t
where v0 is a constant of integration. We use the initial conditions of the problem to set its value. In this case, Joe's speed at time t = 0, when he jumps out of the plane, is v = 0 m/s. So that means v0 = 0 m/s, and so

      v(t)  =   0  -  g*t

Likewise, we can solve for Joe's position as a function of time:

      y(t)  =  y0  +  v0*t  -  0.5*g*t^2
Once again, there are constants of integration. We know that at the start, Joe's height is H = 100 m, and his velocity is v = 0 m/s, so we have

      y(t)  =   H  +   0    -  0.5*g*t^2

Okay, that's the analytic solution. We could use it to solve for Joe's speed and altitude at any time, and we could also figure out the time at which Joe's altitude is zero (ouch!).

Now, let's see how we can solve the problem numerically. Consider Joe's velocity first. We know

      --  =  -g
and we want to find v(t).

Hmmm. Suppose that we know Joe's velocity at some particular instant t0. As we can for any arbitrary function, we can expand his velocity as a function of time in a Taylor series:

                                        1             2
     v(t)  =  v(t0)  +  v'*(t - t0)  +  - v''*(t - t0)   +   ....

In this particular case, we happen to know the first derivative of velocity, v'. If we ignore all the terms of second, third, and higher order, we can calculate Joe's speed an instant later:

     v(t)  =  v(t0)  +  v'*(t - t0)

Note that this solution only gives the correct answer if the second, third and higher-order derivatives are all zero. However, if we pick a small enough time interval, then it will give us a result which is close to the true answer, even without the terms of higher order.

When dealing with numerical solutions, we often write quantities with little subscripts to indicate "current" and "next" values:

         t    =   current time

         t    =   next time

         v    =   current velocity

         v    =   velocity at next time

         v'   =   current derivative of velocity
In which case, we can write:

       v      =  v    +  v' * (t    - t )
        i+1       i       i     i+1    i

In many cases, it helps to define a "timestep:" the amount of time which passes between each of the moments at which we calculate position, velocity, or other physical quantities.

         dt  =  t     - t

Code to solve differential equations almost always contains a loop somewhere in which quantities at the current time are used to calculate quantities at the next instant:


       know   t    v    v'
               i    i    i 

       v     =  v   +  v' * dt
        i+1      i      i     

       v'    =  (something)

       t     =  t   +  dt
        i+1      i

       % now prepare to go through loop for the next moment

       v     =  v 
        i        i+1

       v'    =  v'
        i        i+1

       t     =  t 
        i        i+1


This approach is called Euler's method: it uses a first-order approximation to the Taylor series to calculate new values for physical quantities.

Now, in this particular case, as Joe falls without air resistance, we know

      v' =   --  =  -g 

            d v
     v'' =  ---- =   0
and all higher derivatives of velocity are also zero. So the Taylor series turns into

     v(t)  =  v(t0)  +  (-g)*(t - t0)  +  0  +  0 ...
     v(t)  =  v(t0)  +  (-g)*(t - t0)  
which means that Euler's method

     v     =   v     +  (-g) * (t    - t )
      i+1       i                i+1    i
will give the exact value for velocity. Lucky us! Most of the time, the second (and higher) order term of the Taylor series is not zero, and so Euler's method yields only approximately correct values.

For example, consider the height of Joe above the ground as a function of time: y(t). We can write a Taylor series for it:

                                        1             2
     y(t)  =  y(t0)  +  y'*(t - t0)  +  - y''*(t - t0)   +   ....
Now, in this case,

     y'    =   v

     y''   =   a  =  -g

     y'''  =   0           (and all higher derivs, too)

The second derivative of position with respect to time is not zero (in fact, it's a constant, -9.8 meters per second squared). So, if we apply Euler's method to calculate position as a function of time:

     y     =  y   +  v * (t    - t )  
      i+1      i           i+1    i

Let's try that and see what happens:

  t          a           v             y  
 (s)       (m/s^2)     (m/s)          (m)
  0        -9.8          0            100

  1        -9.8        -9.8        

  2        -9.8       -19.6

  3        -9.8       -29.4

  4        -9.8       -39.2


      Q:   If we use this scheme, 
           how far above the ground is Joe at t = 4 s?

      Q:   Estimate when Joe hits the ground.

But this is not the right answer: the analytic solution indicates that Joe reaches the ground at t = 4.52 s. Uh-oh! What can we do!? Well, we have several options:

  1. We could try to decrease the step size in time. The error we make in using Euler's method is
                      1                  2
            error  =  - * g * (t    - t )
                      2         i+1    i
    If we make the timestep small, then the error -- which goes as the square of the timestep -- will become very small; small enough, perhaps, that it will make no significant difference in our result.

  2. We could try to include the second-order term in the Taylor expansion explicitly in our calculations:
                                              1                   2
         y     =  y   +  v * (t    - t )  -  --- * g * (t    - t )
          i+1      i           i+1    i       2          i+1    i
    In this case, because the second derivative is constant and all higher derivatives are zero, this ought to give us the exact value for the next position. In general, however, the second- and higher-order derivatives are not so simple. Often, the higher derivatives change with time.

  3. We could use a numerical technique more sophisticated than Euler's method. The Runge-Kutta methods, for example, try to estimate the values of higher-order terms in the Taylor series, without evaluating them explicitly. They, too, don't give results which are exactly correct -- but they are certainly more accurate than Euler's method. On the other hand, they are a bit more complicated to put into code, and so it's more likely that the coder (or someone reading the code) might make a mistake...

This week, let's stick with Euler's method, for both velocity and position, to solve the problem. The initial conditions are:

and you know that, without air resistance,

     ----  =  a  =  -g  =  -9.8 m/s^2

Creative Commons License Copyright © Michael Richmond. This work is licensed under a Creative Commons License.