Numerical solution to differential equations in C++, path to take?

Question

Edit

I am now using odeint. It is fairly simple to use and less memory hungry than my brute force algorithm implementation.

Check my questions here-->http://stackoverflow.com/questions/12060111/using-odeint-function-definition

and here-->http://stackoverflow.com/questions/12150160/odeint-streaming-observer-and-related-questions

---

I am trying to implement a numerical method (Explicit Euler) to solve a set of three coupled differential equations. I have worked with C before, but that was a very long time ago (effectively forgotten everything). I have a pretty good idea on what I want my program to do and also have a rough algorithm.

I am interested in using C++ for this task (picked up Stroustroup's Programming: Principles and Practice using C++). My question is, should I go with arrays or vectors? Vectors seem easier to handle, but I was unable to find how you can make a function return a vector? Is it possible for a function to return more than one vector? At this point, I am familiarizing myself with the C++ syntax.

I basically need my function to return many arrays. I realize that it is not possible in C++, so I can also work with some nested structure such as {{arr1},{arr2},{arr3}..}. Please bear with me as I am a noob and come from programming in Mathematica.

Thanks!

Answer 1

If you want to use C++ for integrating ordinary differential equations and you don't want to reinvent the wheel use odeint. This lib is on its way of becoming the de facto standard for solving ODEs in C++. The code is very flexible and highly optimized and can compete with any handcrafted C-code (and Fortran anyway).

Commenting on you question on returning vectors or arrays: Functions can return vectors and arrays if the are wrapped in a class (like std::array). But this is not recommended, since you make many unnecessary copies (incl. calling the constructors and destructors every time). I assume you want to put your function equation into a c++ function and let it return the resulting vector. For this task it's much better if you pass a reference to a vector to the function and let the function fill this vector. This is also the way how odeint has implemented this.

Answer 2

To make the program do what you wan you could take a look at this code, It may be get you started. I found it very useful, and tested it against mathematica solution, and it is ok. for more information go here

/*
      A simple code for option valuation using the explicit forward Euler method
   for the class Derivative Securities, fall 2010
      http://www.math.nyu.edu/faculty/goodman/teaching/DerivSec10/index.html
      Written for this purpose by Jonathan Goodman, instructor.

      Assignment 8

*/


#include <iostream>
#include <fstream>
#include <math.h>
#define NSPOTS  100   /* The number of spot prices computed  */

/*   A program to compute a simple binomial tree price for a European style put option  */

using namespace std;


//    The pricer, main is at the bottom of the file



void FE(                   // Solve a pricing PDE using the forward Euler method

   double T, double sigma, double r, double K,     // The standard option parameters
   double Smin, double Smax,                       // The min and max prices to return
   int nPrices,                                    // The number of prices to compute between Smin and Smax, 
                                                   // Determines the accuracy and the cost of the computation
   double prices[],                                // An array of option prices to be returned.  
   double intrinsic[],                             // The intrinsic value at the same prices
   double spots[],                                 // The corresponding spot prices, computed here for convenience.
                                                   // Both arrays must be allocated in the calling procedure
   double *SEarly ) {                              // The early exercise boundary

//      Setup for the computation, compute computational parameters and allocate the memory

   double xMin = log(Smin);                                 //  Work in the log variable
   double xMax = log(Smax);
   double dx   = ( xMax - xMin )/ ( (double( nPrices - 1 ) ) );  // The number of gaps is one less than the number of prices
   double CFL  = .8;                                        // The time step ratio
   double dt   = CFL*dx*dx/sigma;                           // The forward Euler time step size, to be adjusted slightly
   int    nTimeSteps = (int) (T/dt);                        // The number of time steps, rounded down to the nearest integer
          nTimeSteps++;                                     // Now rounded up
          dt   = T / ( (double) nTimeSteps );               // Adjust the time step to land precisely at T in n steps.
   int    nx   = nPrices + 2*nTimeSteps;                    // The number of prices at the final time, growing by 2 ...
                                                            // ... each time step
          xMin = xMin - nTimeSteps*dx;                      // The x values now start here
   double *fOld;                                            // The values of the pricing function at the old time
           fOld   = new double [nx];                        // Allocated using old style C++ for simplicity
   double *fNew;                                            // The values of the pricing function at the new time
           fNew   = new double [nx];
   double *V;                                               // The intrinsic value = the final condition
           V      = new double [nx];

//       Get the final conditions and the early exercise values 

   double x;      // The log variable
   double S;      // A stock price = exp(x)
   int    j;
   for ( j = 0; j < nx; j++ ) {
      x = xMin + j*dx;
      S = exp(x);
      if ( S < K ) V[j] = K-S;      // A put struck at K
      else         V[j] = 0;
      fOld[j] = V[j];               // The final condition is the intrinsic value
    }

//      The time stepping loop

   double alpha, beta, gamma;                   // The coefficients in the finite difference formula
   alpha = beta = gamma = .333333333333;        //   XXXXXXXXXXXXXXXXXXXXXXXXXXX
   int jMin = 1;                                // The smallest and largest j ... 
   int jMax = nx - 1;                           // ... for which f is updated.  Skip 1 on each end the first time.
   int jEarly    ;                              // The last index of early exercise
   for ( int k = nTimeSteps; k > 0; k-- ) {     // This is, after all, a backward equation
      jEarly = 0;                               // re-initialize the early exercise pointer
      for ( j = jMin; j < jMax; j++ ) {                                // Compute the new values
         x = xMin + j*dx;                                              // In case the coefficients depend on S
         S = exp(x);
         fNew[j] = alpha*fOld[j-1] + beta*fOld[j] + gamma*fOld[j+1];   // Compute the continuation value
         if ( fNew[j] < V[j] ) {
            fNew[j] = V[j];                                            // Take the max with the intrinsic value
            jEarly  = j;                                               // and record the largest early exercise index
          }
       }
      for ( j = jMin; j < jMax; j++ )                                  // Copy the new values back into the old array
         fOld[j] = fNew[j];
      jMin++;                                                          // Move the boundaries in by one
      jMax--;
    }

//     Copy the computed solution into the desired place

   jMin--;      // The last decrement and increment were mistakes
   jMax++;
   int i = 0;                             // The index into the output array
   for ( j = jMin; j < jMax; j++ ) {      // Now the range of j should match the output array
      x = xMin + j*dx;                                              
      S = exp(x);
      prices[i]    = fOld[j];
      intrinsic[i] = V[j];
      spots[i++]   = S;                     // Increment i after all copy operations
    }
   double xEarly = xMin + jEarly*dx;
         *SEarly = exp(xEarly);             // Pass back the computed early exercise boundary

   delete fNew;       // Be a good citizen and free the memory when you're done.
   delete fOld;
   delete V;
   return;
 }

int main() {

   cout << "Hello " << endl;

   ofstream csvFile;              // The file for output, will be csv format for Excel.
   csvFile.open ("PutPrice.csv");

   double sigma = .3;
   double r     = .003;
   double T     = .5;
   double K     = 100;
   double Smin  = 60;
   double Smax  = 180;
   double prices[NSPOTS];
   double intrinsic[NSPOTS];
   double spots[ NSPOTS];
   double SEarly;

   FE( T, sigma, r,  K, Smin, Smax, NSPOTS, prices, intrinsic, spots, &SEarly ); 

   for ( int j = 0; j < NSPOTS; j++ ) {        //   Write out the spot prices for plotting
      csvFile << spots[j];
      if ( j < (NSPOTS - 1) ) csvFile << ", "; //   Don't put a comma after the last value
    }
   csvFile << endl;

   for ( int j = 0; j < NSPOTS; j++ ) {        //   Write out the intrinsic prices for plotting
      csvFile << intrinsic[j];
      if ( j < (NSPOTS - 1) ) csvFile << ", "; //   Don't put a comma after the last value
    }
   csvFile << endl;


   for ( int j = 0; j < NSPOTS; j++ ) {        //   Write out the computed option prices
      csvFile << prices[j];
      if ( j < (NSPOTS - 1) ) csvFile << ", ";
    }
   csvFile << endl;

   csvFile << "Critical price," << SEarly << endl;

   csvFile << "T ," << T << endl;
   csvFile << "r ," << r << endl;
   csvFile << "sigma ," << sigma << endl;
   csvFile << "strike ," << K << endl;

   return 0 ;

 }

Answer 3

This link might help you, but for ordinary differential equations :

http://www.codeproject.com/KB/recipes/odeint.aspx