James Rohal

For you Mathematica and Matlab junkies here are brushes for SyntaxHighlighter. These has ALL the functions (over 2200 for Mathematica 6.0 and 1300 for Matlab) which tends to make browsers stall when rendering the source code. To make use of this brush, just use language=’mathematica’ or language=’matlab’.

Download for SyntaxHighlighter 1.51: shBrushMathematicaAll.js, shBrushMatlabAll.js

Download for SyntaxHighlighter 2.0: shBrushMathematica.js, shBrushMatlab.js

Matlab Example

[sourcecode language='matlab']
{% Example Monte Carlo Simulation in Matlab
Function: y = x2^2/x1Generate n samples from a normal distribution
r = ( randn(n,1) * sd ) + mu
mu : mean
sd : standard deviationGenerate n samples from a uniform distribution
r = a + rand(n,1) * (b-a)
a : minimum
b : maximum
%}
n = 100000; % The number of function evaluations

% — Generate vectors of random inputs
% x1 ~ Normal distribution N(mean=100,sd=5)
% x2 ~ Uniform distribution U(a=5,b=15)

x1 = ( randn(n,1) * 5 ) + 100;
x2 = 5 + rand(n,1) * ( 15 – 5 );

% — Run the simulation
% Note the use of element-wise multiplication

y = x2.^2 ./ x1;

% — Create a histogram of the results (50 bins)

hist(y,50);

% — Calculate summary statistics

y_mean = mean(y)
y_std = std(y)
y_median = median(y)
[/sourcecode]

Mathematica Example

[sourcecode language='mathematica']
(* STEP 1 *)
(* The amount of time remaining is timeRem *)
timeRem = 10;
path = {p};

(* STEP 2 *)
(* While there is still time left *)
While[timeRem > 0,
(* Assume for this iteration that m is the midpoint of the mirror the photons path intersects with *)
m = Round[p + 2/3 v];
(* Find the smallest positive root. If s != \[Infinity] then s is equal to the time it takes for the photon to reach the intersection point with the mirror, otherwise the photon does not intersect the mirror *)
s = Min[Cases[t /. Solve[(p + t v - m).(p + t v - m) == 1/9, t], _?Positive]];
(* If the time it takes to reach the next intersection is less than the time remaining then that means there is an intersection *)
If[s < timeRem,
(* p is now set to be the point at which it intersects with the mirror and v is determined by our transformation matrix H *)
p += s v;
v = H[p - m].v,
(* Else there was no intersection, so there are two things that may have happened: (1) The photon stops in mid-flight at timeRem = 10. In this case we let s = timeRem and thus p has reached its final destination, or (2) the photon has sufficient time, but the mirror we chose initially as determined by m was not correct. In this case we let s be 2/3 since we will choose the next mirror in the photon’s path *)
s = Min[timeRem, 2/3];
p += s v
];
(* STEP 3 *)
(* Reduce timeRem by s *)
timeRem -= s;
(* Append p to path *)
path = Append[path, p];
];

Print["The final position of our photon = ", Last[path]];
Print["Distance from p to the origin = ", answer = Norm[p]];
Print["Number of digits Mathematica claims to have correct = ", Precision[answer]]
[/sourcecode]