Category: Numerical Analysis

Blog Entry © Monday August 4, 2025, by James Pate Williams, Jr. Bessel Functions of the First and Second Kind, Their Derivatives and Zeros

Blog Entry Monday August 04 2025Download

Blog Entry © Friday, August 1, 2025, by James Pate Williams, Jr. Numerically Solving a Two-Dimensional Elliptic Partial Differential Equation (PDE) Boundary Value Problem (BVP) and a New MP3

Blog Entry Friday August 1 2025Download

Blog Entry © Thursday, July 31, 2025, Numerically Solving Four Second Order Linear Ordinary Differential Equation Boundary Value Problems by James Pate Williams, Jr.

Blog Entry Thursday July 31 2025Download

Blog Entry (c) Wednesday, July 30, 2025, by James Pate Williams, Jr. Four Root Finding Algorithms and C Source

function exp(-x) - sin(0.5 * pi * x)
The interval found by bisection is:
a = +0.4435735341  b = +0.4435735341
bisection # iterations = 50
The root by regula falsi:
x = +0.4435735341  f(x) = -0.0000000000
regula falsi # iterations = 8
The root by the secant method:
x = +0.4435735341  f(x) = +0.0000000000
secant method # iterations = 9
The root by Newton's Method:
x = +0.4435735341  f(x) = -0.0000000000
Newton's Method # iterations = 6

D:\roots\x64\Release\roots.exe (process 28152) exited with code 0 (0x0).
Press any key to close this window . . .

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

typedef double real;

static real f(real x)
{
	double pi2 = 2.0 * atan(1.0);
	return(exp(-x) - sin(pi2 * x));
}

static real g(real x)
{
	double pi2 = 2.0 * atan(1.0);
	return(-exp(-x) - pi2 * cos(pi2 * x));
}

static int bisection(
	real(*f)(real), real* a, real* b, real xtol, int* flag)
{
	real error, fa, fm, xm;
	int n = 0;

	fa = (*f)(*a);
	if (fa * (*f)(*b) > 0.0)
	{
		*flag = -1;
		return(n);
	}
	error = fabsl(*b - *a);
	while (error > xtol)
	{
		n++;
		error *= 0.5;
		if (error <= xtol)
		{
			*flag = 0;
			return(n);
		}
		xm = 0.5 * (*a + *b);
		if (xm + error == xm)
		{
			*flag = 1;
			return(n);
		}
		fm = (*f)(xm);
		if (fa * fm > 0.0)
		{
			*a = xm;
			fa = fm;
		}
		else
			*b = xm;
	}
	*flag = 2;
	return(n);
}

static int regula(
	real (*f)(real), real a, real b, real xtol,
	real ftol, int ntol, real* w, int* flag)
{
	int n = 0;
	real fa, fb, fw, signfa, prvsfw;

	fa = f(a);
	if (fa >= 0.0) signfa = +1.0; else signfa = -1.0;
	fb = f(b);
	if (signfa * fb >= 0.0)
	{
		*flag = -1;
		return n;
	}
	*w = a;
	fw = fa;
	for (n = 0; n <= ntol; n++)
	{
		if (fabs(a - b) <= xtol)
		{
			*flag = 0;
			return n;
		}
		if (fabs(fw) <= ftol)
		{
			*flag = 1;
			return n;
		}
		*w = (fa * b - fb * a) / (fa - fb);
		if (fw >= 0.0) prvsfw = +1.0; else prvsfw = -1.0;
		fw = f(*w);
		if (signfa * fw > 0.0)
		{
			a = *w;
			fa = fw;
			if (fw * prvsfw > 0.0) fb = 0.5 * fb;
		}
		else
		{
			b = *w;
			fb = fw;
			if (fw * prvsfw > 0.0) fa = 0.5 * fa;
		}
	}
	*flag = 2;
	return n;
}

static int secantMethod(
	real(*f)(real), real xtol, real ftol,
	int ntol, real xm1, real x0, real* w)
{
	real fm1 = f(xm1), f0 = f(x0);
	real df = fabs(fm1 - f0), f1 = 0.0;
	real dx = fabs(xm1 - x0), x1 = 0.0;
	int n = 0;

	while (n < ntol && df > ftol && dx > xtol)
	{
		x1 = (f0 * xm1 - fm1 * x0) / (f0 - fm1);
		f1 = f(x1);
		df = fabs(f(x1) - f0);
		dx = fabs(x1 - x0);
		fm1 = f0;
		f0 = f1;
		xm1 = x0;
		x0 = x1;
		n++;
	}

	*w = x1;
	return n;
}

static int NewtonsMethod(
	real(*f)(real), real(*g)(real),
	real xtol, real ftol, int ntol,
	real* w)
{
	// f is the function
	// g is the function's derivative
	// xtol is root's tolerance
	// ftol is the function's tolerance
	// ntol is the maximum # of iterations
	real f1 = 0.0, g1 = 0.0, x0 = *w, x1 = 0.0;
	real f0 = f(x0), g0 = g(x0);
	real deltaX = DBL_MAX, deltaF = DBL_MAX;
	int n = 0;

	while (n < ntol && deltaX > xtol && deltaF > ftol)
	{
		x1 = x0 - f0 / g0;
		f1 = f(x1);
		g1 = g(x1);
		deltaX = fabs(x1 - x0);
		deltaF = fabs(f1 - f0);
		f0 = f1;
		g0 = g1;
		x0 = x1;
		n++;
	}

	*w = x1;
	return n;
}

int main(void)
{
	int flag = 0, ntol = 0;
	real a0 = 0, b0 = 0, ftol = 0, w1 = 0, w2 = 1.0;
	real a1 = 0, b1 = 0, w3 = 0, xtol = 0;

	a0 = 0.0;
	b0 = 1.0;
	ntol = 128;
	ftol = 1.0e-15;
	xtol = 1.0e-15;
	int its1 = bisection(f, &a0, &b0, xtol, &flag);
	a1 = 0.0;
	b1 = 1.0;
	int its2 = regula(f, a1, b1, xtol, ftol, ntol, &w1, &flag);
	int its3 = secantMethod(f, ftol, xtol, ntol, 0.0, 1.0, &w2);
	int its4 = NewtonsMethod(f, g, xtol, ftol, ntol, &w3);
	printf("function exp(-x) - sin(0.5 * pi * x)\n");
	printf("The interval found by bisection is:\n");
	printf("a = %+13.10lf  b = %+13.10lf\n", a0, b0);
	printf("bisection # iterations = %ld\n", its1);
	printf("The root by regula falsi:\n");
	printf("x = %+13.10lf  f(x) = %+13.10lf\n", w1, f(w1));
	printf("regula falsi # iterations = %ld\n", its2);
	printf("The root by the secant method:\n");
	printf("x = %+13.10lf  f(x) = %+13.10lf\n", w2, f(w2));
	printf("secant method # iterations = %ld\n", its3);
	printf("The root by Newton's Method:\n");
	printf("x = %+13.10lf  f(x) = %+13.10lf\n", w2, f(w3));
	printf("Newton's Method # iterations = %ld\n", its4);
	return(0);
}

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

Blog Entry Tuesday July 29 2025Download
#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;
}

Blog Entry © Sunday, July 20, 2025, by James Pate Williams, Jr. Example and Two Exercises from a Numerical Analysis Textbook

Happy 59th Anniversary of the historic moon landing by Neil Armstrong, Buzz Aldrin, and Michael Collins in Apolla 11.

Blog Entry July 20 2025Download

Newtonian Gravity and Cosmology © Friday, July 18, 2025, by James Pate Williams, Jr. in Conjunction with the Microsoft Copilot AI

Newtonian Gravity and CosmologyDownload

Historical Attendance Data for the First United Methodist Church (FUMC) of LaGrange, Georgia for the Year 2006 and Modeling the Data Using Various Polynomials and Least Squares Fitting by James Pate Williams, Jr., BA, BS, Master of Software Engineering and Doctor of Philosophy in Computer Science

Church Data Polynomial Fitting 2006Download