Solving differential equations with Monte Carlo method

[aurimas13]

Hey guys, I have been trying to solve three linear differential equations (dx/dt = -2x +2y + 2; dy/dt = x - 4y + 4z + 1; dz/dt y - 8z + 6) using monte carlo method for almost two weeks now, but still have little idea of where to start and what to do. Thus, I want to ask of you for some guidance and little steps in the code, so I could begin writing the code myself. I believe I understand Monte Carlo method, but have difficulty in implementing it to solve the linear differential equations. The problem is as follows: You are given the three differential equations (dx/dt = -2x +2y + 2; dy/dt = x - 4y + 4z + 1; dz/dt y - 8z + 6) and need to solve them to find X, Y and Z values using Monte Carlo method.

If you need more specifics - let me know and I'll provide them even tho I have given what is written in the problem.

So far my code is this:

#include<iostream>
#include<cmath>
#include<cstdlib>

double rnd();

int main()
{
  long max = 1;
  long n = 0;
  int i,j;

  printf("\nMonte Carlo: Random process - Three pots\n\n");
  printf("Exact result:  X = 3   Y = 2   Z = 1\n\n");
  printf("Iter. \t   X \t    Y \t     Z\n\n");

  for (int j=0; j < 16; j++)
  {
          for (int i = 0; i < max; i++)
          {
            //Monte_Carlo_method;
            n++;
          }
          printf("%5ld %8.4lf %8.4lf %8.4lf \n", max, (double)j / n, (double)j / n, (double)j / n);
          max *= 2;
  }
}

double rnd() {
  return (double)rand() / (double)RAND_MAX;
}

[lastchance]

That can't be the whole of your problem. You won't be able to get a solution of a differential equation unless you have a boundary condition, and you haven't specified one.

Personally, I'd diagonalise the matrix on the RHS and solve three uncoupled first-order equations. If you had to do it numerically I'd use standard Runge-Kutta, not a stochastic method.

The only way I can think of doing it by a stochastic (random-number) method is to basically do the Forward-Euler integration by the Monte Carlo method. Integration by Monte Carlo is fairly straightforward - just find the fraction of your random numbers that lie underneath the curve (or, in your case, surface) which will give you the ratio of your integral to an overall area/volume. But Forward-Euler is a lousy method at the best of times.

EDIT: Are you sure that you have really stated your actual problem? In your code you state the "exact solution", yet any solution of your differential equation is a function of t. The numbers are what you get if you set the derivatives to zero; i.e. the long-term solution.
ARE YOU SURE THAT YOUR ACTUAL PROBLEM IS NOT TO SOLVE 3 SIMULTANEOUS LINEAR EQUATIONS?

Please state what your actual problem says, NOT your paraphrasing of what you think it says.

[aurimas13]

Okey, so the problem is as follows:

Think of a random process and write a Monte Carlo program, which finds the solution for threshold/limit of differential equations system: dx/dt = -2x +2y + 2; dy/dt = x - 4y + 4z + 1; dz/dt y - 8z + 6.

The answers for x,y and z should converge to X=3, Y=2 and Z=1.

Hope that helps. I look forward to responses :)

[lastchance]

Go google "Monte Carlo method linear algebra".

Then pour yourself a stiff drink.

The random walk procedure below is based on
https://math.nist.gov/mcsd/Seminars/2014/2014-11-04-Li-presentation.pdf
Half hour to code - 5 hours to understand the method!
In the general case you would have to precondition the matrix (I've just factored the element sizes).

#include <iostream>
#include <iomanip>
#include <random>
#include <chrono>
#include <vector>
using namespace std;

using vec = vector<double>;
using matrix = vector<vec>;

unsigned seed = chrono::system_clock::now().time_since_epoch().count();
mt19937 gen( seed );
uniform_real_distribution<double> dist( 0.0, 1.0 );


//================================================


double myRand() { return dist( gen ); }


//================================================


vec Markov( int num, const matrix &H, const vec &b )
{
   const double SMALL = 1.0e-4;     // Termination criterion for random walk

   int N = H.size();
   vec result( N );

   // Compute row sums, transition matrix and cumulative probabilities
   vec rowSum( N );
   matrix p( N, vec( N ) );
   matrix sump( N, vec( N, 0 ) );

   for ( int i = 0; i < N; i++ )
   {
      rowSum[i] = 0;
      for ( int j = 0; j < N; j++ ) rowSum[i] += abs( H[i][j] );     
      for ( int j = 0; j < N; j++ )
      {
         p[i][j] = abs( H[i][j] ) / rowSum[i];
         if ( j == 0 ) sump[i][j] = p[i][j];
         else          sump[i][j] = sump[i][j-1] + p[i][j];
      }
   }

   // Now solve by random walk for each coordinate
   for ( int i = 0; i < N; i++ )
   {
      double sumX = 0.0;
      for ( int trial = 1; trial <= num; trial++ )
      {
         int index = i;
         double W = 1, X = b[i];                 // W = multiplied weights; X = displacement on random walk
         while ( abs( W ) > SMALL )
         {
            double r = myRand();                 // Random number on [0,1)
            int j = 0;
            while ( r > sump[index][j] ) j++;    // First column where cumulative probability exceeds r

            W *= rowSum[index];
            if ( H[index][j] < 0 ) W = -W;
            X += W * b[j];
            index = j;
         }
         sumX += X;
      }
      result[i] = sumX / num;
   }

   return result;
}


//================================================


int main()
{
   matrix A = { { -2,  2,  0 },
                {  1, -4,  4 },
                {  0,  1, -8 } };
   vec b = { 2, 1, 6 };


   // Put in Ulam and von Neumann's form: x = Hx + b with size of H pretty small
   int N = A.size();
   double Amax = 0.0;
   for ( int i = 0; i < N; i++ )
   {
      for ( int j = 0; j < N; j++ ) if ( abs( A[i][j] ) > Amax ) Amax = abs( A[i][j] );
   }
   for ( int i = 0; i < N; i++ )
   {
      for ( int j = 0; j < N; j++ ) A[i][j] /= Amax;
      b[i] /= Amax;
   }

   matrix H = A;
   for ( int i = 0; i < N; i++ ) H[i][i]++;


   // Solution by Markov chains / random walks
   cout << "num\tx\ty\tz\n";
   cout.precision( 4 );
   for ( int num = 10; num <= 1000000; num *= 10 )
   {
      vec X = Markov( num, H, b );
      cout << num << '\t';
      for ( double x : X ) cout << x << '\t';
      cout << '\n';
   }
}

num	x	y	z
10	3.149	2.723	0.9211	
100	2.919	2.252	0.994	
1000	3.143	1.928	0.9976	
10000	2.981	1.995	1.003	
100000	2.989	1.998	1.001	
1000000	3.001	2	0.9998

[aurimas13]

Thank you lastchance. I'll google what you proposed and will write my own code.

[aurimas13]

I'm trying to run your proposed code to see if it works on my mac, but get 3 warnings and 2 errors:

 warning: alias declarations are a C++11 extension
      [-Wc++11-extensions]
using vec = vector<double>; 
pratimas_10v2.cpp:9:16: warning: alias declarations are a C++11 extension
      [-Wc++11-extensions]
using matrix = vector<vec>;                ^
pratimas_10v2.cpp:82:11: error: non-aggregate type 'matrix' (aka
      'vector<vector<double> >') cannot be initialized with an initializer list
   matrix A = { { -2,  2,  0 },;
pratimas_10v2.cpp:85:8: error: non-aggregate type 'vec' (aka 'vector<double>')
      cannot be initialized with an initializer list
   vec b = { 2, 1, 6 }; pratimas_10v2.cpp:112:22: 
warning: range-based for loop is a C++11 extension
      [-Wc++11-extensions]
      for ( double x : X ) cout << x << '\t';

I assume the compiler needs updating. Anyway I'll try googling everything, but please let me know if you know how to fix these errors/warnings :)

[lastchance]

Just run it in C++ shell (gear-wheel icon top-right of code).

Then upgrade your compiler.

[aurimas13]

I was right. The compiler needs to be set to C++11 and then it fixes everything.

[aurimas13]

What you mean c++ shell. How I could locate it in mac terminal?

[lastchance]

I don't use a Mac.

On a PC (and, indeed, on an iPad) if you look at a runnable code sample in a post on this forum then, provided the poster has put the sample between code tags, a little gear-wheel icon should come up to the top right of the sample. Click on it.

Sometimes it takes a while to load the shell or run.

[aurimas13]

Understood. Thanks. That's a great feature to have in the forum :)

[aurimas13]

Hi lastchance. I got back to your written code and would like to ask of you (if possible) to explain it step by step as I can't figure everything out. It seems a bit complex :/

Looking foward to your reply.

[lastchance]

Well, I'll do my best. I'm basically relying on somebody's random-walk algorithm at one point.

First your original problem:

Think of a random process and write a Monte Carlo program, which finds the solution for threshold/limit of differential equations system: dx/dt = -2x +2y + 2; dy/dt = x - 4y + 4z + 1; dz/dt y - 8z + 6.

The long-time limit will be when the system goes to steady state. (This system does reach a steady state, but others can grow exponentially.) Then the time derivatives dx/dt, dy/dt and dz/dt tend to zero, so in the long-time limit your equations become the algebraic system
0 = -2x +2y + 2
0 = x - 4y + 4z + 1
0 = y - 8z + 6
In matrix/vector form you could write this as
0 = Ax + b
and you will find A and b hard-coded in my code as

   matrix A = { { -2,  2,  0 },
                {  1, -4,  4 },
                {  0,  1, -8 } };
   vec b = { 2, 1, 6 };

(Note that I would usually write something like Ax=b, so be careful with the signs here.)

The method coded wants that written slightly differently, as x = Hx + b so add x (or Ix, where I is the identity matrix), to either side:
x = (A+I)x + b
or
x = Hx + b
This is done with the lines

   matrix H = A;
   for ( int i = 0; i < N; i++ ) H[i][i]++;

Note that the method won't always converge; you need the "size" of H to be quite small. (Sorry, I'm not going into eigenvalue spectra - just note that simply dividing all elements of A and b by Amax in my lines 90-99 was attempting to keep H small in some sense.)


OK, so you have now got a system of the form
x = Hx + b
Each component of x (i.e. x, y, z) can then be shown to be the "expected value"; i.e. mean of the displacement X at the end of a sequence of random walks, with probabilities at each step related to the elements in H - see slides 10,11,12 in
https://math.nist.gov/mcsd/Seminars/2014/2014-11-04-Li-presentation.pdf
(Sorry, these don't constitute a proof - you will have to do quite a lot of research online to work that out: I don't intend to compress my scribbled notes. I managed to follow the first three pages of
http://www.cs.ubbcluj.ro/~studia-m/2006-1/rosca.pdf
before getting lost, but it adds a bit more detail to the other link.)

Slides 13-27 in the first reference go through (in excruciating detail) a SINGLE random walk. The following code (lines 55-67 of the original) corresponds to a single random walk as set out in, e.g., slide 13.

         int index = i;
         double W = 1, X = b[i];                 // W = multiplied weights; X = displacement on random walk
         while ( abs( W ) > SMALL )
         {
            double r = myRand();                 // Random number on [0,1)
            int j = 0;
            while ( r > sump[index][j] ) j++;    // First column where cumulative probability exceeds r

            W *= rowSum[index];
            if ( H[index][j] < 0 ) W = -W;
            X += W * b[j];
            index = j;
         }

Then I average over a lot of random walks:

         sumX += X;
      }
      result[i] = sumX / num;

Sorry, I'm too shaky on the algorithm itself to write anything more substantive than that.


This seems a horrendously complex problem to do with no other guidance. A Monte Carlo method is definitely not what I would usually think of for solving linear equations - straight Gaussian elimination would do.


Maybe there's some better stochastic method ... perhaps somebody else on the forum can help. My original google search was
"Monte Carlo method linear algebra".


Maybe also I've gone off in completely the wrong direction, so it might be advisable to contact your teacher (nagging gets a surprisingly long way at universities) and ask him/her what was intended. I would hate to send you on a wild goose chase. This seems a real sledgehammer method.

[lastchance]

Are you absolutely sure that your teacher intended a "Monte-Carlo" (i.e. random-number-based) approach, and not an "iterative" approach ... for which there are lots? For example, this is a simple Forward-Euler iteration to steady state. It is definitely NOT a Monte-Carlo approach: there is no randomness anywhere in it.

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

using vec = vector<double>;
using matrix = vector<vec>;

//================================================


int main()
{
   const double EPSILON = 1.0e-6;                // Tolerance for convergence
   const int    ITERMAX = 1000000;               // Need to go home some time

   matrix A = { { -2,  2,  0 },
                {  1, -4,  4 },
                {  0,  1, -8 } };
   vec b = { 2, 1, 6 };

   int N = A.size();
   vec x = b;                // Will hold solution

   // Solve (by timestepping) dX/dt = Ax + b for the long-time limit (i.e. dX/dt -> 0)
   double change = 1;
   double dt = 0.1;
   int iter = 0;
   while ( change > EPSILON && iter < ITERMAX )
   {
      iter++;
      change = 0;
      for ( int i = 0; i < N; i++ )
      {
         // Calculate Ax+b (ith component)
         double RHS = b[i];
         for ( int j = 0; j < N; j++ ) RHS += A[i][j] * x[j];
         x[i] += RHS * dt;
         change += abs( RHS );
      }
   }

   if ( iter == ITERMAX ) cout << "Maximum iterations reached.\n";

   cout << "Solution: ";
   for ( double e : x ) cout << e << " ";
   cout << '\n';
}

Solution: 3 2 1