Category: Numerical Analysis
Blog Entry © Sunday, August 10, 2025, First-Order Perturbation Treatment of the Helium Atom by James Pate Williams, Jr.
Blog Entry © Wednesday, August 6, 2025, C/C++ Orthogonal Polynomials by James Pate Williams, Jr.
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 © 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 © Thursday, July 31, 2025, Numerically Solving Four Second Order Linear Ordinary Differential Equation Boundary Value Problems by James Pate Williams, Jr.
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.
#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;
}
Numerical Explorations 