Contact Us!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutSign UpSign In

Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place. Commercial Alternative to JupyterHub.

| Download
Views: 29
Visibility: Unlisted (only visible to those who know the link)
Image: ubuntu2204
Tweet
Share
Share
Kernel: Python 3 (system-wide)

Solving Differential Equations with the Euler-Cromer Method

You will encounter differential equations that cannot be solved analytically. However, there are numerical methods that can be used to solve them with a computer. For many second-order, ordinary differential equations, a fairly simple approach called the Euler-Cromer Method is sufficient. This method was used frequently in Physics 231, but the name may not have been mentioned. (Note that it is referred to as the "last point approximation" or LPA in the refrence below.)

Notation

We will use shorthand notation for time derivatives where

x˙=dxdt\dot{x} = \frac{dx}{dt}

and

x¨=d2xdt2.\ddot{x} = \frac{d^2x}{dt^2}.

In other words, the number of dots indicates the number of time derivatives. We’ll be finding the value of x(t)x(t) at discrete values of tt. Label these times with integer subscripts as

t1=t0+Δtt2=t0+2Δt⋮ti=t0+i  Δtt_1 = t_0 + \Delta t \\ t_2 = t_0 + 2\Delta t \\ \vdots \\ t_i = t_0 + i\;\Delta t

where t0t_0 is the initial time and Δt\Delta t is the time interval between the subsequent solutions. We will also label the values of the function at these time with integers as

x0=x(t0)x1=x(t1)⋮xi=x(ti).x_0 = x(t_0) \\ x_1 = x(t_1) \\ \vdots \\ x_i = x(t_i).

Similarly, the derivatives at step ii are

x˙i=x˙(ti)\dot{x}_i = \dot{x}(t_i)

and

x¨i=x¨(ti).\ddot{x}_i = \ddot{x}(t_i).

The Problem

The form of the second-order differential equation to be solved is

d2xdt2=f(x,dxdt,t),\frac{d^2x}{dt^2} = f\left(x, \frac{dx}{dt}, t\right),

or, in the more compact notation, it is

x¨=f(x,x˙,t).\ddot{x} = f\left(x, \dot{x}, t\right).

With the initial values of the variable (x0x_0) and its first derivative (x˙0\dot{x}_0) are known at the initial time (t0t_0), we want to find xx and x˙\dot{x} at later times.

Euler Method

The simplest numerical method to solve differential equations is the Euler Method. The derivatives at step ii can be approximated as

x˙i=dxdt∣ti≈x(ti+1)−x(ti)Δt=xi+1−xiΔt\dot{x}_i = \frac{dx}{dt}\bigg|_{t_i} \approx \frac{x(t_{i+1}) - x(t_{i})}{\Delta t} = \frac{x_{i+1} - x_i}{\Delta t}

and

x¨i=d2xdt2∣ti≈x˙(ti+1)−x˙(ti)Δt=x˙i+1−x˙iΔt.\ddot{x}_i = \frac{d^2x}{dt^2}\bigg|_{t_i} \approx \frac{\dot{x}(t_{i+1}) - \dot{x}(t_{i})}{\Delta t} = \frac{\dot{x}_{i+1} - \dot{x}_i}{\Delta t} .

Suppose that xx and x˙\dot{x} are known at some time. Rearranging the equation for the first derivative, the approximate the solution after moving a step forward in time is

xi+1≈xi+x˙iΔt.x_{i+1} \approx x_i + \dot{x}_i\Delta t.

To find the solution at another step forward in time, we will need to know the first derivative at that later time. Rearranging the equation for the second derivative, the approximate first derivative after moving a step forward in time is

x˙i+1≈x˙i+x¨iΔt=x˙i+f(xi,x˙i,ti)Δt,\dot{x}_{i+1} \approx \dot{x}_i + \ddot{x}_i\Delta t = \dot{x}_i + f(x_i,\dot{x}_i,t_i)\Delta t,

where we insert the function for the second derivative. These final two equations can be used iteratively to find the approximate solution at later times. Unfortunately, this method is not very accurate after a while. (For classical mechanics problems, this method often doesn't conserve energy.)

Euler-Cromer Method

A simple modification to the method described above greatly improves the accuracy of the numerical solution. First, calculate the value of the second derivative at the current time,

x¨i=f(xi,x˙i,ti),\ddot{x}_i = f(x_i,\dot{x}_i,t_i),

and use it to approximate first derivative after a step forward in time,

x˙i+1≈x˙i+x¨iΔt.\dot{x}_{i+1} \approx \dot{x}_i + \ddot{x}_i\Delta t.

Next, the value of the varible at a later time is approximated by

xi+1≈xi+x˙i+1Δt,x_{i+1} \approx x_i + \dot{x}_{i+1}\Delta t,

where the first derivative at the later time is used. Note the slight difference between this approach and the Euler Method. The current value of x¨\ddot{x} (which may depend on the current values of xx, x˙\dot{x}, and tt) is used to calculate a new value of x˙\dot{x}. However, the new value of x˙\dot{x} is used to calculate the new value of xx.

Implementation in Python

As an example, we can numerically solve the differential equation

d2xdt2=−kx,\frac{d^2x}{dt^2} = -kx,

for k=5k=5 and with initial conditions x=10x=10 and x˙=0\dot{x}=0 at t=0t=0. In the program, the first and second derivatives of xx will be called xdxd and xddxdd, respectively. (Thinks of each "dd" as a dot in the shorthand notation for derivatives.) The first part of the program prepares for the calculation. The pylab library is imported so that the results can be plotted. Constant(s), initial values of xx, xdxd, and xddxdd, time step dtdt, and the final time tftf are set. Three lists (t_listt\_list, x_listx\_list, and xd_listxd\_list) are started with the inital values of tt, xx, and xdxd.

import pylab as pl

# constants
k = 5

# set inital values, time step, & final time
t = 0
x = 10
xd = 0
dt = 0.001
tf = 10

# start lists with the initial values
t_list = [t] 
x_list = [x]
xd_list = [xd]

The main part of the program is the loop that calculates xx and its first derivative (xdxd) at time steps of dtdt using the Euler-Cromer method. The results for each step are appended to the lists. For this example, the calculation continues until the time reaches tftf.

while (t < tf):
    # calculate new values
    xdd = -k*x       # Calculate d2x/dt2 using the current x & dx/dt
    xd = xd + xdd*dt # Use d2x/dt2 to update dx/dt
    x = x + xd*dt   # Use the updated dx/dt to update x
    t = t + dt      # Advance t by a step dt
    # append new values to lists
    t_list.append(t)
    x_list.append(x)
    xd_list.append(xd)

The results were stored in lists so that they can be used. After the calculations are complete, xx and x˙=dx/dt\dot{x}=dx/dt can plotted versus tt as shown below.

# plot the x and dx/dt vs. time (from the lists)
pl.figure()
pl.plot(t_list,x_list, ls='-', c='b')
pl.xlabel('t')
pl.ylabel('x')

pl.figure()
pl.plot(t_list,xd_list, ls='-', c='b')
pl.xlabel('t')
pl.ylabel('dx/dt')

pl.show()
Image in a Jupyter notebookImage in a Jupyter notebook

It is important to check that the time step used was small enough to get accurate results. If dtdt is small enough, changing it slightly won't change the results significantly.

Note that this method can be used for ordinary, second-order differential equations, even if the function depends of something other than time. For example, it could be used to solve

d2zdx2=3z−dzdx\frac{d^2z}{dx^2} = 3z - \frac{dz}{dx}

for z(x)z(x).

Reference

Alan Cromer, "Stable solutions using the Euler approximation," Am. J. Phys. 49, 455-459 (1981). (Note: As mentioned in the article, this method was discovered accidentally by a high school student.)