Serious mathematics, written with the reader in mind.
[MsE] Matrix Editions

 Monte Carlo Program This program requires a Pascal compiler. The nine lines beginning function Rand and ending end are a random number generator. The six lines beginning function randomfunction  and ending end define a random function that gives the absolute value of the determinant of a 2 x 2 matric with entries a,b,c,d.  For an n x n matrix, you would enter n squared "seeds".  You can name them what you like; if n=3, you could call them x1, x2, ..., x9 instead of a, b,...i.  In that case you would write x1:=rand(seed); x2:=rand(seed), and so on.  To define the random function you would use the forrmula for the determinant of a 3 x 2 matrix. Here is the program. After the program we give an example showing how to use it. program  montecarlo; const lengthofrun = 100000;  var S,V,x,intguess,varguess,stddev,squarerootlerun,  errfifty,errninety,errninetyfive:longreal;  i,seed,answer:longint; function Rand (var Seed: longint): real; {Generate pseudo random number between 0 and 1} const  Modulus = 65536;       Multiplier = 25173;      Increment = 13849; begin  Seed := ((Multiplier * Seed) + Increment)  mod Modulus;  Rand := Seed / Modulus end; function randomfunction:real; var a,b,c,d:real; begin  a:=rand(seed);b:=rand(seed);c:=rand(seed);d:=rand(seed);  randomfunction:=abs(a*d-b*c); end; begin Seed := 21877; repeat S:=0;V:=0; for i:=1 to lengthofrun do  begin       x:=randomfunction;       S:=S+x;V:=V+sqr(x);  end;         intguess:=S/lengthofrun; varguess:=V/lengthofrun-sqr(intguess); stddev:= sqrt(varguess);  squarerootlerun:=sqrt(lengthofrun);  errfifty:= 0.6745*stddev/squarerootlerun; errninety:= 1.645*stddev/squarerootlerun; errninetyfive:= 1.960*stddev/squarerootlerun; writeln('average for this run = ',intguess); writeln('estimated standard deviation = ',stddev); writeln('with probability 50%, your error is less than ',errfifty);  writeln('with probability 90%, your error is less than ',errninety);  writeln('with probability 95%, your error is less than ',errninetyfive); writeln('another run? 1 with new seed, 2 without'); readln(answer);  if (answer=1) then  begin  writeln('enter a new seed, which should be an integer');  readln(seed);  end; until (answer=0); end. Example of using the Monte Carlo program: Note: In Pascal, x^2 is written sqr(x). To compute the area inside the unit square and above the parabola y=x^2,  you would type function functiomfunction:real; var x,y:real; begin x:=rand(seed);y:=rand(seed); if (y-sqr(x) <0 ) then randomfunction:=0 else randomfunction:=1; end

Last updated March 12, 2003