Numerical Methods for Solving Initial Value Problems

Overview

I warned you that this blog would be random, and now we will discuss an undergraduate math topic. Recently, I was cataloging some of my old text books and I came across one on differential equations. The title is Differential Equations Computing and Modeling, by C. Henry Edwards and David E. Penney. It occurred to me that while I had done some small computation projects in this class using formulas in Excel and some MATLAB programs, most people do not have access to MATLAB (at least legally), though they should have access to a C or C++ compiler. So, without further ado, here is a small project (single source file) with multiple numerical methods for differential equation solution.

Introduction

First we need to define some terminology.

• ODE - Ordinary differential equation. This is an equation which contains functions of only one independent variable, and derivatives with respect to that variable.
• IVP - Initial value problem. This is an ODE, together with an initial condition. For example, dy/dx = y, y(0) = 1 is an IVP. (Coincidentally, if you solve this for y(1), you will determine Euler's costant e).
• Derivative - What is wrong with you? Go back to Calculus 1. Sorry, that was an inside joke.
  

Next, in order to solve your problem using these numeric methods (as implemented), you need to express your differential equation such that y'=dy/dx is on one side. For example:

y' - y = x

should be expressed as

dy/dx = y + x

Once this is done, you will be implementing dy/dx as a function ( that's a c++ function here) which takes two doubles and returns a double. Specifically, the signature is as below:

// This is your differential equation y'(x,y)... that is
// formulate the function dy/dx as a function of x,y
typedef double (*Y_primed)(double x, double y);

So, to implement the example function that was given, you could do something like the following

// This is the example y', following the signature above
double sample_y_primed(double x, double y) {
return x + y;
}

The Numerical Methods

Numerical methods are used for solving differential equations in an approximate way. There are three numerical methods used in this source. Space will not allow me to cover them thoroughly here, but here are some Wikipedia resources. You can find the algorithms there as well.

• Euler's Method - created by Leonhard Euler, mathematical genius. http://en.wikipedia.org/wiki/Euler_method
• Improved Euler's Method - also called Heun's Method, after Karl M. W. L. Heund. This method is an extension of the Euler Method. http://en.wikipedia.org/wiki/Heun%27s_method
  
• Runge-Kutta Method - A numerical method created by C. Runge and M.W. Kutta which has much lower error for a given number of steps, and thus has better accuracy for the same number of iterations as Euler's Method. Alternatively, a lot of computational efficiency can be gained at the expense of a little bit of accuracy when compared to Euler's Method. http://en.wikipedia.org/wiki/Runge-Kutta_methods

The implemented numerical methods all follow this prototype
double SomeNumericalMethod(Y_primed f, double initial_x, double initial_y, double final_x, unsigned int steps)

... where f, is the y' function covered earlier, initial_x is the x-value at which the initial condition is true, initial_y is the y value at the initial condition, final_x is the x-value at which to calculate the new y value, and steps is the number of iterations. The approximate new y-value, that is the approximate y(final_x), is returned.

Here is a small example:

dy/dx = x + y, y(0) = 1

We shall call ImprovedEulersMethod() function

double y_of_1 = ImprovedEulersMethod(sample_y_primed, 0, 1, 1, 10);

You should find that y_of_1 is near 3.4282, which is approximately equal to the real y-value y(1) = 3.4366. Increase the number of steps (20, 40, 80, ...) to see how the accuracy increases.

Downloading

Here's a link to the C++ source code.

Unexplored Ideas

1. What is the error of each method as a function of the number of iterations used?
2. How many iterations does it take to have n-decimal places accurate?
3. Derive and prove why the special IVPs and final x-values for e, ln(2), and pi arrive at those values.