C++ Coding - Euler's Method

Question Consider an initial value ODE of the following form

write a function to return the value of y at x=b given.

For this example we are going to implement the Euler method as given by

where w at the nth step gives an estimate for the value of y at x=b.

As with all programs we start by thinking about what are the parameters and local variables in the problem. It is clear from the specification here that the parameters are a, b, n, as well as the function to be integrated itself although as we are not interested in writing a generic algorithm we can ignore the last one. Start with an empty program with libraries for input/output to screen/file and mathematical functions.

#include <iostream>
#include <fstream>
#include <cmath>
using namespace std;

int main()
{
  return 0;
}

Compile and run the program to check you have everything set up ok.

Next shall declare the external parameters first, followed by the local parameters and then initialise any values at the start of the algorithm.

  // parameters
  int n=10;
  double a=0.,b=1.,alpha=0.;
  // local variables
  double h,x,w;
  // intialise values
  h=(b-a)/(double)(n);
  x=a;
  w=alpha;

Finally we can add in the rest of the algorithm in a couple of lines and output results to the screen at the same time.

  // implement Euler's method
  for(int i=0;i<n;i++)
  {
    x = a + i*h; // update value of x to x_i
    cout << x << " " << w << endl; // output x_i and w_i to screen
    w = w + h*(x*exp(3*x)-2.*w); // update w to w_{i+1}
  }
  cout << b << " " << w << endl; // output x_n and w_n to screen

The output from this program should give y(x=1)=2.7609. Now that the program is running we should check the value of the result for different values of n, and see if the results appear to converge to a constant value as the number of steps n increases.

Now we can move this code into its own function, so that we can be more flexible with obtaining results. If we use the external parameters as arguments to the function and then copy/paste the function in from "main" it will look like:

double eulersMethod(int n,double a,double b,double alpha)
{
  // local variables
  double h,x,w;
  // intialise values
  h=(b-a)/(double)(n);
  x=a;
  w=alpha;
  // implement Euler's method
  for(int i=0;i<n;i++)
  {
    x = a + i*h; // update value of x to x_i
    cout << x << " " << w << endl; // output x_i and w_i to screen
    w = w + h*(x*exp(3*x)-2.*w); // update w to w_{i+1}
  }
  cout << b << " " << w << endl; // output x_n and w_n to screen
  return w;
}

and can be used inside "main" like

  cout << " y(b) ~ " << eulersMethod(1000,0.,1.,0.) << endl;

If you want to output the results of x against w to a file for plotting then we can use the flexibility of the stream variable by including it as an argument to the function. We add this argument as ostream& output (must use a non constant reference) and then where cout appears in the function it can be replaced by output. The final code could look something like:

#include <iostream>
#include <fstream>
#include <cmath>
using namespace std;

double eulersMethod(int n,double a,double b,double alpha,ostream& output)
{
  // local variables
  double h,x,w;
  // intialise values
  h=(b-a)/(double)(n);
  x=a;
  w=alpha;
  // implement Euler's method
  for(int i=0;i<n;i++)
  {
    x = a + i*h; // update value of x to x_i
    output << x << " , " << w << endl; // output x_i and w_i to screen
    w = w + h*(x*exp(3*x)-2.*w); // update w to w_{i+1}
  }
  output << b << " , " << w << endl; // output x_n and w_n to screen
  return w;
}

int main()
{
  // output to screen
  eulersMethod(1000,0.,1.,0.,cout);
  // output to file
  ofstream myFileStream("test.csv");
  eulersMethod(1000,0.,1.,0.,myFileStream);
  myFileStream.close();
  return 0;
}