Initial value problems

This topic extends what was done before on integrals to solutions of certain kinds of differential equations, or initial value problems.

1. Euler's method

The problem we solve here is a first order differential equation with initial value,

y' = f (x, y)

with

y (x_{0}) = y_{0}

where the function $f$ and the initial value $y_{0}$ at $x_{0}$ are known, and given by some C++ function and variables.

Our numerical "solution" will not be an alegraic representation of the required function, but will be a table of data giving for certain values $x$ an estimate or approximation of the corresponding value $y$ . This table will only contain such data for certain $x$ -values. Once we have the table we might obtain $y (x)$ for values $x$ outside the table by the method of polynomial interpolation.

The very simplest such method is named after Euler.

Euler's method.

We are given a function for $f$ and an initial value $(x_{0}, y_{0})$ . Take some small quantity $h$ . Suppose $(x_{i}, y_{i})$ has been generated and is on the curve required, or at least not too far from it. Then the gradient at $(x_{i}, y_{i})$ is $f (x_{i}, y_{i})$ and the curve goes through the point $(x_{i}, y_{i})$ so we can approximate the curve by a straight line,

y - y_{i} = f (x_{i}, y_{i}) (x - x_{i})

giving an estimate for the $y$ value at $x_{i + 1} = x_{i} + h$ as approximately

y_{i + 1} = y_{i} + h f (x_{i}, y_{i})

This recurrence generates a table of data, and we expect the table of values $x_{i}, y_{i}$ to be close to the required function.

2. Programming Euler's method

Since it will be convenient to use Euler's method on a range of functions $f$ , we can load the functional package (which requires C++ in at least the 2011 variety). Note that $f$ is a function of two double arguments.

#include <iostream>
#include <functional>

using namespace std;

/** Perform one step of Euler's method.
 *  Given y' = f(x,y) where y(xsub0) = ysub0 
 *  returns the value ysub1 = ysub0 + h*f(xsub0,ysub0)
 */
double eulerstep(const function<double(double,double)>& f, double xsub0, double ysub0, double h) {
  ... 
}

/** Perform n steps of Euler's method.
 *  Given y' = f(x,y) where y(xsub0) = ysub0 
 *  returns the value ysubn given by
 *    ysub(i+1) = ysub(i) + h*f(xsub(i),ysub(i))
 *    xsub(i+1) = xsub(i) + h
 *  as it computes this, the function also tabulates all pairs of values
 *    xsubi , ysubi
 *  printing these to cout.  
 */
double euler(const function<double(double,double)>& f, double xsub0, double ysub0, double h, int n) {
  ... 
}

Exercise.

Program the above functions. Write a suitable main() function that tests them out for the differential equations

$r' = r$ , where $r = r (x)$ and $r (0) = 1$ in the range $x \in [0, 1]$
$s' = - x / s$ , where $s = s (x)$ and $s (0) = 1$ in the range $x \in [0, 1]$

For these you should program suitable functions $f_{r} (x, r)$ and $f_{s} (x, s)$ where $r' = f_{r} (x, r)$ and $s' = f_{s} (x, s)$ .

double f_r(double x, double r) { ... }
double f_s(double x, double s) { ... }

Exercise.

Now comment out the cout lines in the euler function so that it does not produce any output and try to use the Euler method to estimate $r (1)$ and $s (1)$ for the solutions to the differential equations given above. Try numbers of steps equal to 10, 100, 1000, 10000, 100000, 1000000. Use a main similar to the following.

#include <iostream>
#include <functional>
#include <cmath>

using namespace std;

// your functions go here

int main() {
  int nsteps = 0;
  cout << "Number of steps?" << endl;
  cin >> nsteps;
  double h = 1.0 / nsteps;

  double r = euler(f_r,0,1,h,nsteps);
  double r_exact =  // put an expression for the exact answer here
  cout << "r(1)==" << r << "; actual==" << r_exact << " error==" << abs(r-r_exact) << endl;

  double s = euler(f_s,0,1,h,nsteps);
  double s_exact =  // put an expression for the exact answer here
  cout << "s(1)==" << s << "; actual==" << s_exact << " error==" << abs(s-s_exact) << endl;

  return 0;
}

3. Alternatives

There isn't time to go into alternatives to the Euler method here, but suffice it to say here that there are a huge number of them. By and large, the trick is to approximate the curve by a polynomial rather than by a straight line using ideas related to polynomial interpolation. For example, this is what the popular Runge-Kutta RK4 method does. Another clever idea is to vary the step size $h$ : for some non-critical parts of the curve a large step size might be accurate enough but for others a smaller step size may be needed.