Integration

Simpson One Dimensional - Implementation 6.1

Simpson’s method - takes a function f, and integrates it from a to b in Ns steps - there is no check that Ns is even, which it must be to correctly handle the boundaries.

Simpson[f_, a_, b_, Ns_] := Module[{Ia, j, dx},
dx = (b-a)/Ns;
Ia = 0.0;
For[j =0, j ≤ Ns-2, j = j+2,
Ia = Ia + (1/3) dx (f[a + j dx] + 4 f[a + (j+1) dx] + f[a + (j+2) dx]);
];
Return[Ia];
]

A simple test to ensure basically correct behavior:

f[x_] := x Sin[x]

test = Simpson[f,0.0,1.0,100]

exact = N[-Cos[1] + Sin[1]]

test - exact

The Simpson method must be used with an even number of grid points, otherwise the upper limit, Ns-2 in the sum (the For loop above), is not achieved - as an example, the same test as above, but with 99 grid points:

testbad = Simpson[f,0.0,1.0,99]

Higher Dimension

Two dimensions - Implementation 6.2, 6.3

Here we send in f, the limits in x, ax, bx, the number of steps for the x-integration, Nx, and all the same information for the y-integration.

Simpson2D[f_,ax_,bx_,Nx_,ay_,by_,Ny_] := Module[{Iy},
Iy[x_] := Module[{Igrandy},
Igrandy[y_] := f[x,y];
Return[Simpson[Igrandy,ay,by,Ny]];
];
Return[Simpson[Iy,ax,bx,Nx]];
]

f[x_,y_] := Sin[x] Cos[y](x^2 + y^2)

test = Simpson2D[f, 0,10,100,0,10,100]

Compare with the built-in numerical integrator:

ntest = NIntegrate[Sin[x] Cos[y](x^2 + y^2),{x,0,10},{y,0,10}]

Abs[test - ntest]

They’re relatively close, and we only used one hundred points in the Simpson’s method, so there’s plenty of room for improvement.

The next function integrates the user-specified function f(x,y) over a circle of radius R - the circular boundaries are handled in the limits for the y-integration.  There is a final division by π (R^2) to return the average value of f over the circle.

CircleAvg[f_,R_, Ns_] := Module[{Iy},
Iy[x_] := Module[{Igrandy},
Igrandy[y_] := f[x,y];
Return[Simpson[Igrandy,-Sqrt[R^2 -x^2],Sqrt[R^2 - x^2],Ns]];
];
Return[Simpson[Iy,-R,R,Ns]/(Pi R^2)];
]

First we’ll check with the obvious constant function:

f[x_,y_] := 1

CircleAvg[f,1.0,100]

f[x_,y_] := (x y)^2

CircleAvg[f,1.0,200]

In circular coordinates, this integral is just (radius one):  (1/π)∫(r Cos[T] r Sin[T])^2 rdrdT, which ends up being 1/24:

N[1/24]

So the error is:

Abs[CircleAvg[f,1.0,200] - N[1/24]]

Three dimensions - Implementation 6.4

More of the same here - we have a function f, now of three variables, then specify the (Cartesian) limits and number of steps in each of the three directions.

Simpson3D[f_,ax_,bx_,Nx_,ay_,by_,Ny_, az_,bz_,Nz_] := Module[{Ixy},
Ixy[z_] := Module[{Igrandxy},
Igrandxy[x_,y_] := f[x,y,z];
Return[Simpson2D[Igrandxy,ax,bx,Nx,ay,by,Ny]];
];
Return[Simpson[Ixy,az,bz,Nz]];
]

Here’s a test function:

f[x_,y_,z_] := x^2 - x y + z^3

test = Simpson3D[f, 0,1,32,-1/2,1/2,32,-1/4,3/4,32]

The correct number here (computing the integral directly and inserting the limits) is: 79/192, so we can compute the residual here:

test - N[79/192]

Electric field outside of a cube of (constant) charge - Figure 6.4

We’ll make a coarse approximation, with the number of steps set at 8 in each direction.

Ns = 8;

Next, the functions we’ll be using to calculate the x, y, and z components of the electric field:

Ex[x_,y_,z_,pref_] := Module[{igrand,ret},
igrand[xp_,yp_,zp_] := (pref) ((x-xp)/( (x-xp)^2 + (y-yp)^2 + (z-zp)^2)^(3/2));
ret = Simpson3D[igrand,-1/2,1/2,Ns,-1/2,1/2,Ns,-1/2,1/2,Ns];
If[ret === Indeterminate, ret = 0;];
Return[ret];
]

Ey[x_,y_,z_,pref_] := Module[{igrand,ret},
igrand[xp_,yp_,zp_] := (pref) ((y-yp)/( (x-xp)^2 + (y-yp)^2 + (z-zp)^2)^(3/2));
ret = Simpson3D[igrand,-1/2,1/2,Ns,-1/2,1/2,Ns,-1/2,1/2,Ns];
If[ret === Indeterminate, ret = 0;];
Return[ret];
]

Ez[x_,y_,z_,pref_] := Module[{igrand,ret},
igrand[xp_,yp_,zp_] := (pref) ((z-zp)/( (x-xp)^2 + (y-yp)^2 + (z-zp)^2)^(3/2));
ret = Simpson3D[igrand,-1/2,1/2,Ns,-1/2,1/2,Ns,-1/2,1/2,Ns];
If[ret === Indeterminate, ret = 0;];
Return[ret];
]

Now we’ll make a grid of values (step size ds) and compute the electric field at each point.

ds = N[1/4];

Inside the cube, one can get zero separations, that’s okay provided we ignore them:

Off[Power::infy];
Off[Infinity::indet];

EEtab = Table[{ {x,y,z},{Ex[x,y,z,1],Ey[x,y,z,1],Ez[x,y,z,1]}},{x,-2,2,ds},{y,-2,2,ds},{z,-2,2,ds}];

On[Power::infy];
On[Infinity::indet];

G1 = ListVectorPlot3D[EEtab,VectorStyle->GrayLevel[0],VectorPoints->Coarse]

G0 = Graphics3D[{Glow[Gray],Black,Cuboid[{-1/2,-1/2,-1/2},{1/2,1/2,1/2}]}];

Show[G0,G1]

Monte Carlo

We’ll seed the random number generator, so that the output always matches.

SeedRandom[11111]

Circle areas - Implementation 6.5

Below, we throw darts at a circular dartboard set inside a square with side length 2 R.

MCArea[R_, Ns_] := Module[{hit, index,x,y},
hit = 0;
For[index= 1, index ≤ Ns, index=index+1,
x = RandomReal[{-R,R}];
y = RandomReal[{-R,R}];
If[Sqrt[x^2 + y^2] ≤ R, 
hit = hit + 1;
];
];
Return[(hit/Ns) (2 R)^2];
]

MCArea[2.0,100000]

N[Pi (2)^2]

Integrals

Two different views of Monte Carlo integration - the first is the standard story for Monte Carlo, and the second is (6.80)

MCInt[f_, a_,b_, fmax_, Ns_] := Module[{hit,index,pt,val},
hit = 0;
For[index = 1, index ≤ Ns, index=index+1,
pt = RandomReal[{a,b}];
val = RandomReal[{0,fmax}];
If[val ≤ f[pt],
hit = hit+1;
];
];
Return[(hit/Ns) (b-a) fmax];
];

MCInt2[f_,a_,b_, Ns_] := Module[{retval,pt, dx,index},
dx = (b-a)/Ns;
retval = 0.0;
For[index = 1, index ≤ Ns, index=index+1,
pt = RandomReal[{a,b}];
retval = retval + dx f[pt];
];
Return[retval];
]

f[x_] := Sin[x]/x

N[MCInt[f,0,1,1,100000]]

MCInt2[f,0,1,100000]

The two Monte Carlo forms agree, and both approximate the “actual” value:

N[SinIntegral[1]]

Example - Figure 6.9

Here, we’ll make a plot of the mean, standard deviation (on either side of the mean), and the actual approximations generated by 1000 runs of Monte Carlo to approximate the integral of Sin[x]/x as above - this gives us a way to attach error bars to the Monte Carlo approximation:

data = Table[N[MCInt[f,0,1,1,1000]],{j,1,1000}];

Show[{ListPlot[data,PlotStyle->GrayLevel[0.5]], Plot[Mean[data],{j,1,1000},PlotStyle->GrayLevel[0.1]],Plot[Mean[data] + Sqrt[Variance[data]],{j,1,1000},PlotStyle->GrayLevel[0.1]],Plot[Mean[data] - Sqrt[Variance[data]],{j,1,1000},PlotStyle->GrayLevel[0.1]], Plot[SinIntegral[1],{j,1,1000},PlotStyle->GrayLevel[0]]}]

Mean[data]

Transformation - Figure 6.10

SeedRandom[11111]

We start with a flat distribution on [0,1], and end with a quadratic distribution via (6.93).

f[y_] := -6 y (y-1)

y[x_] :=(1/4) (2+((1-I (Sqrt[3]))/(((-1+2 x+2 (Sqrt[(-1+x) x])))^(1/3)))+(1+I (Sqrt[3])) (((-1+2 x+2 (Sqrt[(-1+x) x])))^(1/3)))

Plot[y[x],{x,0,1}]

ydata = Table[Chop[y[RandomReal[{0,1}]]],{j,1,50000}];

Histogram[ydata,ChartBaseStyle->GrayLevel[0.5]]