How can i cover the stochastic integral on the negative side

z = 1; (*z is constant set to 1*)
f[j_] := j/(z^2 + j^2)^(3/2);
Integrate[f[j], j]
xlast = 5.0; xinit = -2.5;
pltf = Plot[f[j], {j, xinit, xlast}, Filling -> Axis];
Print[pltf]
nintegratef = NIntegrate[f[j], {j, xinit, xlast}];
Print[nintegratef]

findmax = FindMaximum[f[j], {j, xlast}]
fmax = findmax[[1]];
L = xlast – xinit;
A = fmax*L; (* build the rectangular box with dimension fmax x L *)

totalcountmax = 10000;
For[i = 1, i <= totalcountmax, i++, inbox[i] = {0, 0}; missed[i] = {0, 0}; ]; inboxcount = 0; countmissed = 0; For[i = 1, i <= totalcountmax, i++, xrand = RandomReal[{xinit, xlast}]; yrand = RandomReal[{0, fmax}]; If[0 <= yrand <= f[xrand], inbox[i] = {xrand, yrand}; inboxcount = inboxcount + 1; ,(*else, record the missed points*) missed[i] = {xrand, yrand}; countmissed = countmissed + 1; ];(*end if*) Countinbox[i] = inboxcount; Countmissed[i] = countmissed; ];(*end for i*) tinbox[q_] := Table[inbox[i], {i, 1, q}]; tmissed[q_] := Table[missed[i], {i, 1, q}]; gp[q_] := Graphics[{Green, PointSize[0.01], Point[tinbox[q]]}]; gpm[q_] := Graphics[{Yellow, PointSize[0.01], Point[tmissed[q]]} ]; pltfj[q_] := Plot[f[j], {j, xinit, xlast}, Frame -> True,
PlotLabel -> {q, Style[Countmissed[q], Yellow],
Style[Countinbox[q], Green]
}];

Manipulate[Show[pltfj[q], gpm[q], gp[q]], {q, 1, totalcountmax, 1}]
integratestochastic = A*inboxcount/totalcountmax;
Print[“{nintegratef,integratestochastic,A,totalcountmax,inboxcount}”, \
{nintegratef, integratestochastic, A, totalcountmax, inboxcount}];

=================

=================

1 Answer
1

=================

You have to change your For loop to account for positive and negative regions, and make the box twice as big. First change the definition of A,

A = 2 fmax*L;

Then modify your For loop to be

posinboxcount = 0;
neginboxcount = 0;
countmissed = 0;
For[i = 1, i <= totalcountmax, i++, xrand = RandomReal[{xinit, xlast}]; yrand = RandomReal[{-fmax, fmax}]; If[0 <= Sign[f[xrand]] yrand <= Sign[f[xrand]] f[xrand], inbox[i] = {xrand, yrand}; If[Sign[yrand] == 1, posinboxcount = posinboxcount + 1;, neginboxcount = neginboxcount + 1;] ,(*else, record the missed points*) missed[i] = {xrand, yrand}; countmissed = countmissed + 1; ];(*end if*) Countinbox[i] = posinboxcount + neginboxcount; Countmissed[i] = countmissed; ];(*end for i*) and finally your definition for the integral, integratestochastic = A*(posinboxcount - neginboxcount)/totalcountmax; Now you will get If you want to write this more compactly, taking advantage of the built-in functions in Mathematica, then you could calculate the area like region = Rectangle[{xinit, -fmax}, {xlast, fmax}]; sample = RandomPoint[region, 100000]; poscounts = Select[sample, 0 < #[[2]] < f[#[[1]]] &]; negcounts = Select[sample, 0 > #[[2]] > f[#[[1]]] &];
Area[region] (Length@poscounts – Length@negcounts)/Length[sample]

(* 0.173436 *)

And you can make a graphics like

xlast = 5.0;
xinit = -2.5;
z = 1;
f[j_] := j/(z^2 + j^2)^(3/2);
nintegratef = NIntegrate[f[j], {j, xinit, xlast}];
findmax = FindMaximum[f[j], {j, xlast}]
fmax = findmax[[1]];
region = Rectangle[{xinit, -fmax}, {xlast, fmax}];
sample = RandomPoint[region, 10000];
Manipulate[
Module[{poscounts, negcounts, missed},
poscounts = Select[sample[[;; q]], 0 < #[[2]] < f[#[[1]]] &]; negcounts = Select[sample[[;; q]], 0 > #[[2]] > f[#[[1]]] &];
missed = Complement[sample[[;; q]], poscounts, negcounts];
Column[{
Show[
Plot[f[j], {j, xinit, xlast}, Frame -> True, ImageSize -> 400],
Graphics[{Red, PointSize[0.01], Point@negcounts,
Blue, PointSize[0.01], Point@poscounts,
Yellow, PointSize[0.01], Point@missed}]],
Grid[{{“Npoints = “, q},
{“Positive counts = “, Length[poscounts]},
{“Negative counts = “, Length[negcounts]},
{“Calculated area = “,
Area[region] (Length@poscounts – Length@negcounts)/q},
{“Actual area = “, nintegratef}}]
}]
]
, {{q, Length@sample}, 1, Length@sample, 1}]