Blog Entry © Tuesday, July 29, 2025, Double and Triple Monte Carlo Integration by James Pate Williams, Jr.

#include <math.h>
#include <stdio.h>
#include <stdlib.h>

static double randomRange(double lo, double hi)
{
	return (hi - lo) * (double)rand() / RAND_MAX + lo;
}

static double integrand(double r, double w)
{
	return pow(r, 4.0) * (2.0 - r) * w * w * exp(-r);
}

static double StarkEffectIntegral(double E, int N)
{
	double sum = 0.0;

	for (int i = 0; i <= N; i++)
	{
		double r = randomRange(0.0, 100.0);
		double w = randomRange(-1.0, 1.0);

		sum += integrand(r, w);
	}

	return 100.0 * 2.0 * E * sum / (16.0 * (N - 1));
}

static void firstOrderStarkEffect(double E)
{
	double exact = -3.0 * E;
	int N[9] = {
		1000000, 2000000, 3000000, 4000000,
		5000000, 6000000, 7000000, 8000000,
		9000000 };

	for (int n = 0; n < 9; n++)
	{
		int iN = N[n];
		double integ = StarkEffectIntegral(E, iN);
		double error = 100.0 * fabs(integ - exact) / fabs(exact);

		printf("N = %4ld\tintegral = %13.10lf\t%% error = %13.10lf\n",
			iN, integ, error);
	}

	printf("exact value = %13.10lf\n", exact);
}

static double ee1(int N, double R, double Z)
{
	double pi = 4.0 * atan(1.0);
	double sum = 0.0;

	for (int i = 0; i <= N; i++)
	{
		double r1 = randomRange(1.0e-25, R);
		double r2 = randomRange(0.0, r1);

		sum += R * r1 * r1 * exp(-2.0 * Z * (r1 + r2)) * r2 * r2;
	}

	return 16.0 * pi * pi * sum / (N - 1);
}

static double ee2(int N, double R, double Z)
{
	double pi = 4.0 * atan(1.0);
	double sum = 0.0;

	for (int i = 0; i <= N; i++)
	{
		double r1 = randomRange(1.0e-25, R);
		double r2 = randomRange(r1, R);
		
		sum += R * (R - r2) * r2 * exp(-2.0 * Z * (r1 + r2)) * r1 * r1;
	}

	return 16.0 * pi * pi * sum / (N - 1);
}

static void firstOrderHelium(double Z)
{
	double pi = 4.0 * atan(1.0), R = 25.0;
	double exact = 5.0 * pi * pi / (8.0 * pow(Z, 5.0));

	int N[9] = {
		1000000, 2000000, 3000000, 4000000,
		5000000, 6000000, 7000000, 8000000,
		9000000 };

	for (int n = 0; n < 9; n++)
	{
		int iN = N[n];
		double integ = ee1(iN, R, Z) + ee2(iN, R, Z);
		double error = 100.0 * fabs(integ - exact) / fabs(exact);

		printf("N = %4ld\tintegral = %13.10lf\t%% error = %13.10lf\n",
			iN, integ, error);
	}

	printf("exact value = %13.10lf\n", exact);
}

int main(void)
{
	firstOrderStarkEffect(2.0);
	firstOrderHelium(2.0);
	return 0;
}

#include <math.h>
#include <stdio.h>
#include <stdlib.h>

static double randomRange(double lo, double hi)
{
	return (hi - lo) * (double)rand() / RAND_MAX + lo;
}

static double f(double x, double y, double z)
{
	return pow(sin(x), 2.0) + y * sin(z);
}

static double g(double x, double y, double z)
{
	return x + y * z * z;
}

static double integral(
	double x0, double x1,
	double y0, double y1,
	double z0, double z1,
	double (*f)(double, double, double),
	int N)
{
	double sum = 0.0;

	for (int n = 0; n <= N; n++)
	{
		double x = randomRange(x0, x1);
		double y = randomRange(y0, y1);
		double z = randomRange(z0, z1);

		sum += f(x, y, z);
	}

	return (x1 - x0) * (y1 - y0) * (z1 - z0) *
		sum / (N - 1);
}

int main(void)
{
	double pi = 4.0 * atan(1.0);
	double x0 = 0.0, x1 = pi;
	double y0 = 0.0, y1 = 1.0;
	double z0 = 0.0, z1 = pi;
	double exact = 0.5 * pi * (2.0 + pi);
	int N[9] = {
		1000000, 2000000, 3000000, 4000000,
		5000000, 6000000, 7000000, 8000000,
		9000000 };

	printf("integrand pow(sin(x), 2.0) + y * sin(z)\n");
	printf("x = 0 to pi, y = 0 to 1, z = 0 to pi\n");

	for (int n = 0; n < 9; n++)
	{
		int iN = N[n];
		double integ = integral(
			x0, x1, y0, y1, z0, z1, f, iN);
		double error = 100.0 * fabs(integ - exact) / fabs(exact);

		printf("N = %4ld\tintegral = %13.10lf\t%% error = %13.10lf\n",
			iN, integ, error);
	}

	printf("exact value = %13.10lf\n", exact);

	x0 = -1.0;
	x1 = 5.0;
	y0 = 2.0;
	y1 = 4.0;
	z0 = 0.0;
	z1 = 1.0;
	exact = 36.0;

	printf("integrand x + y * z * z\n");
	printf("x = -1 to 5, y = 2 to 4, z = 0 to 1\n");

	for (int n = 0; n < 9; n++)
	{
		int iN = N[n];
		double integ = integral(
			x0, x1, y0, y1, z0, z1, g, iN);
		double error = 100.0 * fabs(integ - exact) / fabs(exact);

		printf("N = %4ld\tintegral = %13.10lf\t%% error = %13.10lf\n",
			iN, integ, error);
	}

	printf("exact value = %13.10lf\n", exact);
	return 0.0;
}