Numerical approximations to first order ODEs

This page is intended to be a supplement to the lecture on numerical approximations to first order ODEs for the UNC course MATH383: First Course in Differential Equations. The general first order ODE

    \[ \frac{dy}{dt} = f(t,y) \]

with an initial condition

    \[ y(t_0) = y_0 \]

defines the general first order initial value problem. Here we will visualize numerical approximations to these initial value problems using Euler’s method (a.k.a. forward Euler). We will visualize these approximations with a Python script. Python is a highly user-friendly programming language, and you won’t actually need to know how to program in order to use the script. The requirements for running the script are a working installation of Python 2.7, along with the SciPy stack of modules. For instructions on installing these requirements, look at Getting Started with Python and SciPy.

Running the script

Download the script and rename it to “forward_euler.py”. When you run the script the first time, you should get a plot that looks like this:

first_run

By default, the script visualizes numerical approximations to the first order ODE

\frac{dy}{dt} = -y,

with the initial condition

y(t_0) = 1.

This could describe the mass of an isotope undergoing radioactive decay.

On the plot:

  • The black arrows show the slope field for the given right-hand side.
  • The blue dots show the approximate solution to the initial value problem using forward Euler. The blue lines are only there to join the blue dots as a visual aid.
  • The red line shows an approximate solution to the initial value problem using a highly precise numerical approximation. This is precise enough that it will be visually indistinguishable from the exact solution.

The buttons above the plot allow you to move the plot, zoom in and out, and save the figure.

Playing with the script

The code between lines 9 and 24 can be modified to change the right-hand side, the initial condition, the step size for forward Euler, and the final time.

Changing the right-hand side and initial condition

First, let’s change our ODE to

\frac{dy}{dt} = \frac{1- y^2 + t}{2 - t}.

In Python, any line that starts with “#” is a comment line and is ignored by the Python interpreter. We can change our ODE by first commenting out line 12, by putting a “#” before “return -y”. Second, we uncomment line 13, by removing the “#” before “#return (1.0 – y**2 + t) / (2.0 – t)”. Lines 12 and 13 should now read

 #return -y
 return (1.0 - y**2 + t) / (2.0 - t)

Now let’s change our initial condition to

y(t_0) = 0.

We can do this by commenting out line 16, and uncommenting line 17 so that the code reads

#y_0 = 1
y_0 = 0

If we run the script now, we should get something that looks like this:

second_run

Changing the final time

Now let’s change the final time to

t_n = 1.8.

This problem now resembles the equation for the velocity in the rocket flight example from class, albeit with very simple coefficients. The mass of the rocket is given by

m(t) = 2 - t,

and the equation will be valid until the rocket runs out of fuel at time t=1.8 s. In the script we can set this by changing line 23 to

t_n = 1.8

If we run the script now, we should get something that looks like this:

thrid_run

Changing the step size

Now let’s change the step size for forward Euler. If we make the step size larger, by changing line 20 to

h = 0.3

we get a less precise approximate solution:

fourth_run If we make the step size smaller, by changing line 20 to

h = 0.03

we get a more precise approximate solution:

fifth_run

Changing whatever you feel like

The code is meant to be played around with. Change the initial value problem, the step size, the end time, or any of the guts of the code to your heart’s content.