Category: Numerical Analysis

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

Historical Attendance Data for the Four Services Offered by the First United Methodist Church of LaGrange, Georgia, Compiled and Analyzed by James Pate Williams, Jr.

Church Attendance BPNN 2006 to 2008Download

Blog Entry © Thursday, June 5, 2025, by James Pate Williams, Jr., Analytic, Numeric, and Siacci’s Method for Solving Ballistic Trajectory Problems (Point-Mass Projectile Motion)

Blog Entry Thursday June 5 2025Download

Blog Entry (c) Monday, June 2, 2025, Exercises from “Exterior ballistics, 1935” by James Pate Williams, Jr.

Below are some Exercises from the textbook Exterior ballistics, 1935 by former Lieutenant Commander Ernest Edward Herrman of the United States Navy Naval Academy in Annapolis, Maryland, see Chapter Four:

Density         Published Density

1.183473        1.183472
1.306584        1.306582
1.202410        1.202408
0.694163        0.694162

Exercises from Exterior ballistics, 1935
By Lieutenant Commander Ernest Edward Herrmann
Exercise 1 Page 41

Book            Mine

0.157100        0.157563
0.492680        0.492810
0.517340        0.517829
0.675240        0.675580
0.844370        0.844597
1.018110        1.018197
1.014010        1.014027
1.094820        1.094910

Percentage Differences Exercise 1 Page 41

Percentage Difference [1] = 0.294283%
Percentage Difference [2] = 0.026383%
Percentage Difference [3] = 0.094477%
Percentage Difference [4] = 0.050340%
Percentage Difference [5] = 0.026880%
Percentage Difference [6] = 0.008545%
Percentage Difference [7] = 0.001676%
Percentage Difference [8] = 0.008220%

Three Altitude Related Density Calculations
Rho1 is from Exterior ballistics, 1935
Rho2 is from NASA
Rho3 is from Wikipedia

h ft    Rho1 SI         Rho2 SI         Rho3

0       1.000000        1.226614        1.208993
1000    0.968902        1.215727        1.174010
2000    0.938772        1.204914        1.139806
3000    0.909578        1.194175        1.106369
4000    0.881292        1.183510        1.073689
5000    0.853886        1.172917        1.041752
6000    0.827332        1.162397        1.010547
7000    0.801604        1.151949        0.980062
8000    0.776676        1.141574        0.950286

Three Altitude Related Ballistic Density Calculations
Rho1 is from Exterior ballistics, 1935
Rho2 is from NASA
Rho3 is from Wikipedia
STP = Temperature = 59 F Pressure = 29.53 Humidity = 78%

Temp    Press   h ft    Rho1 SI         Rho2 SI         Rho3

65      29.60   1000    0.968902        1.203134        1.161848
85      29.75   18000   0.566291        0.992154        0.656245
57      30.25   8000    0.776676        1.174418        0.977626
69      29.80   13000   0.663193        1.077573        0.801819
32      30.15   15000   0.622587        1.157366        0.822138

Exercise 3 Problems 1 - 4 m below means mine
Ra Ram are the Mayevski Retardations
Rf Rfm are the forces of Air Resistence
V Vm are the Applicable Velocities

Ra      Ram     Rf      Rfm     V       Vm

567.20  567.22   229.28  229.19 2626    2626
323.85  323.86   503.50  503.29 3114    3114
 91.58   91.58  2477.50 2476.50 2862    2862
 62.02   62.02  4049.40 4047.83 2584    2584

Exercise 4 - Problems 1 to 5
Ra      Ram     BC              Density

293.82  306.84   3.110358       0.989712
160.92  185.47   4.737834       1.009200
119.06  123.90   6.991926       0.952504
 73.58   84.81  10.328256       0.987977
 60.52   69.59  12.442559       1.080786

Exercise 5 - Problems 1 to 4
Ra      Gv      i-book          i-mine

511.76   803.10 1.001500        0.959553
334.04  1071.14 0.602720        0.578983
 90.62   922.06 0.601870        0.577421
 57.30   799.71 0.616870        0.593688

Test Case from Exterior ballistics, 1935
Temperature 84 F        Pressure 29.90 In Hg
Density SI = 0.959428   Book Density = 0.960
Altitude 18,000 Feet
Density SI = 0.999439   Book Density = 0.991

The differences above are probably due to the fact that Lieutenant Commander Herrmann used logarithms and logarithm tables and perhaps a slide-rule. I use double precision C/C++ real numbers.