Stochastic Differential Equations

Multiple Noise Realizations

This notebook is used to approximate solutions to stochastic differential equations.  There are two version of this notebook.    This version approximates the solution of a fixed system with different sets of noise data.  The other version (called "stosolve.nb" and found in a separate directory called "/flow")  fixes the noise data so that it is easy to look at the flow by varying the initial conditions. The program uses a C library to do the numerical approximation and Mathematica as a front end and for plotting and doing statistics with the results.

Acknowledgements

The programming for this Mathematica module was done by David Glickenstein under the supervision of Ruth Williams.  It was developed through a generous contribution to the UCSD Mathematics Department from the GenCorp Foundation of Aerojet, a corporate donor committed to supporting educational improvement.

The underlying C software written for working with the stochastic differential equations was written by Jessica Gaines, while working at the University of Edinburgh, under the direction of Terry Lyons, and using funding provided by EC grant SC1*CT92-0784.  The software is to be used for educational and academic research purposes only.  Anyone interested in more information should email jessica@maths.ed.ac.uk or to write to Jessica Gaines, University of Edinburgh, Dept. of Mathematics and Statistics, JCMB, Kings Buildings, Mayfield Road, Edinburgh EH9 3JZ UK  and then get original copies directly from Edinburgh themselves.

Disclaimer

The author disclaims any responsibility for problems experienced by users as a result of using this software.  The software may not be distributed without written permission from the author.  For more information contact:
David Glickenstein		or		Ruth Williams
glicken @math.ucsd.edu			williams@math.ucsd.edu

at:	UCSD Mathematics Department
	9500 Gilman Drive 0112
	La Jolla, CA 92093-0112

Introduction

This notebook is used to find the Ito solution of \(TraditionalForm
      \`dX\  = \ \(b(X, t)\)\ dt\  + \ \(σ(X, t)\) dW\) where b is the Ito drift, σ is the diffusion, t is time, and W is the noise, and X is a vector (b and σ are matrix valued).  It also can find Stratonovich solutions of \(TraditionalForm
      \`dX\  = \ \(a(X, t)\) dt\  + \ \(σ(X, t)\)∘dW\).  It uses a combination of C and Mathematica.

It is important to be careful when entering input, since C is slightly different from Mathematica.  Make sure to read the directions in each section, as they are a guideline as to how to use the program.  When working with vectors, remember that the jth element of a vector v is v[j] in C and v[[j]] in Mathematica.

The C library was entirely developed by Jessica Gaines and the programs are largely based upon her example programs.

David Glickenstein, UCSD

Initialization

These initialization routines must be run each time a new Mathematica kernel is run, but the code should not be altered at all.

Load Packages

<< Graphics`Graphics3D`

<< Statistics`DescriptiveStatistics`

Turn off unnecessary warnings

Off[General::spell1];

Set colors

\(Blue\  = \ RGBColor[0, 0, 1]; 
Green\  = \ RGBColor[0, 1, 0]; 

    Red\  = \ RGBColor[1, 0, 0]; 
Yellow\  = \ RGBColor[1, 1, 0]; 

    Cyan\  = \ RGBColor[0, 1, 1]; 
Magenta\  = \ RGBColor[1, 0, 1]; 

    White\  = \ RGBColor[1, 1, 1]; 
Black\  = \ RGBColor[0, 0, 0]; \)

Define function PathGraph2D

\(\(PathGraph2D\  = \ 
      Function[
		{l, pathnumber, \ varx, \ vary, \ color}, 
		
        If[\((pathnumber\  <= \ 0)\)\  || \ \((pathnumber\  > \ Length[l])\), 
          
			Print["Error: There are ", \ Length[l], " paths."], 
          \ 
			
          If[\((Min[varx, vary]\  <= \ 0)\)\  || 
              \((Max[varx, \ vary]\  > \ 
                  Length[\(l[\([pathnumber]\)]\)[\([1]\)]])\)\ , 

            						
            Print["Error: There are ", \ 
              Length[\(l[\([pathnumber]\)]\)[\([1]\)]], \ " variables."], 
            
					
            ListPlot[
              Map[Function[q, \ q[\([{varx, \ vary}]\)]], 
							\ 
                l[\([pathnumber]\)]], 
						PlotJoined -> True, \ 

              						PlotStyle\  -> \ color, 
						
              DisplayFunction -> Identity]]]]; \)\)

Define function PathGraph3D

\(\(PathGraph3D\  = \ 
      Function[{l, \ pathnumber, \ varx, \ vary, \ varz, \ color}, \ 		
        If[\((pathnumber\  <= \ 0)\)\  || \ \((pathnumber\  > \ Length[l])\), 
          
			Print["Error: There are ", \ Length[l], " paths."], 
          \ 
			
          If[\((Min[varx, vary, \ varz]\  <= \ 0)\)\  || 
              \((Max[varx, \ vary, \ varz]\  > \ 
                  Length[\(l[\([pathnumber]\)]\)[\([1]\)]])\)\ , 

            						
            Print["Error: There are ", \ 
              Length[\(l[\([pathnumber]\)]\)[\([1]\)]], \ " variables."], 
            
				
            ScatterPlot3D[
              Map[Function[q, \ q[\([{varx, \ vary, \ varz}]\)]], \ 
                l[\([pathnumber]\)]], \ 
					PlotJoined -> True, \ 

              					PlotStyle -> color, \ 
					
              DisplayFunction -> Identity]]]]; \)\)

Open parameters file

We first need to open the file which contains the parameters which determine the system to solve and how we go about solving it.  The file is called "stoparams.c" and is a C file.  It must be in your current directory.  The command I included here opens the file using emacs, but it is not important which text editor is used.  The next step is to edit this file appropriately and MAKE SURE TO SAVE  the file before continuing.

! "emacs stoparams2.c&"

Generate new noise

In order to generate new noise files, you must delete the old ones.  New noise files (called bpaths1, bpaths2, ...) will be generated when you run stosolve.  If you change the seed for the random number generator, you need to run this line or the change will not take effect.

! "/bin/rm bpaths*"

Compile and run stosolve

This command compiles the C code with the parameters you have just entered.  Any time you change the parameters be sure to save the file and then run this command.  It is important that "makefile" be in the current directory.

! make

This command runs the C program.  It reads noise from the (binary) file "bpaths" and stores the solutions in a (text) file called "soln."

! stosolve2

Load solutions

We need to create a list of the points in the solution paths.  This produces list (called p) of each path, reading from the file "soln."  Each path consists of  a list of points.  Each point has first coordinate time and then the coordinates of the path in order.

\(\(
solnfile\  = \ OpenRead["soln"]; 
ml\  = \ {}; 

    r\  = \ empty; 

    While[r\  =!= \ EndOfFile, 
	r\  = \ empty; 
	l\  = \ {}; 
	
      While[\((r\  =!= \ EndOfFile)\)\  && \ \((r\  =!= \ ":")\), 
		
        r\  = \ Read[solnfile, \ Word]; 
		
        If[\((r\  === \ ":")\)\  || \ \((r\  === \ EndOfFile)\), \ 
          ml\  = \ Append[ml, l], l\  = \ Append[l, ToExpression[r]]]
	]]; 
    
Close[solnfile]; 
p\  = \ {}; 

    For[j\  = \ 1, \ j\  <= \ Length[ml], \ \(j++\), 
	q\  = \ {}; 
	
      For[i\  = \ 1, \ i < \ Length[ml[\([j]\)]], \ i\  += 3, \ 
        q\  = \ Append[q, \ \(ml[\([j]\)]\)[\([{i, \ i + 1, \ i + 2}]\)]\ ]]; 
      
		p\  = \ Append[p, q]]; \)\)

Graph paths in 2D and 3D

Introduction to graphing commands

PathGraph2D, PathGraph3D

PathGraph2D and PathGraph3D graph the paths in 2D or 3D and have the following options:
PathGraph2D[l, pathnumber, varx, vary, color], 
PathGraph3D[l, pathnumber, varx, vary, varz, color]
	1. The parameter l is the list of paths which are loaded above in the "Load solutions" section.
	2. The parameter pathnumber is the number of the path which you want to graph (e.g. 1 is the first path).
	3. The parameters varx, vary, varz are the numbers of the coordinates you want to graph, with 1 being time, 2 being the first coordinate of the solution, 3 being the second coordinate in the solution, etc. (e.g. for the phase space of the first two solution coordinates, we use varx = 2 and vary = 3).
	4. The parameter color is the color you want the path.  The options are Red, Blue, Green, Yellow, Cyan, Magenta, White, and Black.  You can also specify RGB values by setting the color to be RGBColor[a,b,c] where a,b,c are the RGB values.

Show

The Show command displays the graphs.  The first parameters are the graphs (generated with PathGraph2D and PathGraph3D) which you want to show on the same axes.  The next parameters are options.  Among the options you must specify DisplayFunction->$DisplayFunction or the graphs will not be displayed.  
When viewing 3D graphs, it is often advisable to include the option ViewPoint->[a,b,c] where [a,b,c] is the point from which you are viewing the graph.  As this point is not at all obvious, the best way to use this option is to go to the Input Menu and select  "3D ViewPoint selector" (or <Shift>-<Mod1>-V).  This window allows you to select your viewpoint by the mouse.  Then press the Paste button to copy the viewpoint to the notebook and use it in the Show command.

Some Tips and Shortcuts

Note that graphics objects can be assigned to variables, so the command 
	graph1 = PathGraph2D[p, 1, 1, 2, Blue];
makes perfect sense to Mathematica and allows you now to run commands like Show[graph1, DisplayFunction->$DisplayFunction].
If you want to show the same thing as you displayed in the previous command (perhaps with a change in ViewPoint or some other option), you can use the command Show[%, DisplayFunction->$DisplayFunction] (and you can change or add options) and the % simply takes same graphs that were shown last time.

Graph

An example of graphing 2 paths in 2D.

\(\(Show[PathGraph2D[p, 2, 2, 3, Blue], \ PathGraph2D[p, 1, 2, 3, Red], 

      	PathGraph2D[p, 3, 2, 3, Green], DisplayFunction -> $DisplayFunction]; 
    \)\)

An example of a graph in 3D

\(\(Show[PathGraph3D[p, 1, 1, 2, 3, Magenta], 
      ViewPoint -> {2.199, \ \(-1.020\), \ 2.361}, 
	
      DisplayFunction -> $DisplayFunction]; \)\)

Statistics

Get final values (lastpcoords)

We make a list called lastpcoords such that listpcoords[[i]] is a list of the final values of the ith coordinate (as described above, i.e. 1 means time, 2 means first coordinate, etc.) of the solutions (one entry for each path).

\(lastp\  = \ Map[Last, \ p]; 
lastpcoords\  = \ {}; 

    For[i = 0, \ i < \ Length[lastp[\([1]\)]], \ \(i++\), 
	
      lastpcoords\  = \ 
        Append[lastpcoords, \ Map[Function[q, q[\([i + 1]\)]], \ lastp]]]; 
    \)

Get statistics

We can have functions Mean, Median, Mode, GeometricMean, HarmonicMean, RootMeanSquare, SampleRange, Variance, StandardDeviation, MeanDeviation, MedianDeviation, InterquartileRange, and others (see Mathematica help, under Add-ons).

Example: find mean of first coordinate of solutions

Mean[lastpcoords[[2]]]

Example: find median of second coordinate of solutions

Median[lastpcoords[[3]]]

Example: find mean solution vector

The command Map takes a function and a list and applies the function to each element of the list, so Map[f, {a,b,c}] is {f(a), f(b), f(c)}.  This allows us to apply statistics to whole vectors.

Map[Mean, lastpcoords]