#pragma once
/*
* MuellersMethod.c (c) Sunday, July 21, 2024 by
* James Pate Williams, Jr. BA, BS, MSwE, PhD
* Translated from the FORTRAN 77 source code
* found in "Numerical Analysis: An Algorithmic
* Approach" by S. D. Conte and Carl de Boor
* Originally coded in FORTRAN IV in 1982 then
* into C# in March 2015 Finished Tuesday,
* July 23, 2024 The complex division method is
* from "A Numerical Library in C for Scientists
* and Engineers" by H. T. Lau Chapter 1
* Translated to better C++ on Saturday, March 21,
* 2026
*/
#include <string>
#include <vector>
class Mueller
{
private:
static double epsilon1, epsilon2;
public:
static double f(std::string name,
double x, int n, int m,
std::vector<double>& zeros);
static double g(
std::string name,
std::string graph,
double x, int n, int m,
std::vector<double>& zeros);
static double OrthogonalPolynomialRoot(
std::string name,
std::string graph, int n,
double epsilon1, double epsilon2,
int maxIt, int m, double x0,
double x1, double x2, std::string& msg,
std::vector<double>& zeros);
static void OrthogonalPolynomialRoots(
std::string name,
std::string graph, int n, int nr,
unsigned int seed, std::string& msg,
std::vector<double>& zeros);
};
#include "Mueller.h"
#include <cmath>
#include <iostream>
#include <string>
#include <vector>
double Mueller::epsilon1 = 1.0e-12;
double Mueller::epsilon2 = 1.0e-15;
double Mueller::f(std::string name,
double x, int n, int m,
std::vector<double>& zeros)
{
double a1n = 1, a2n = 1, a3n = 1, a4n = 1;
double f0 = 1, f1 = x, f2 = 1;
if (name == "Laguerre")
f1 = -x + 1;
for (int i = 2; i <= n; i++)
{
if (name == "Chebyshev")
{
a1n = 1.0;
a2n = 0.0;
a3n = 2.0;
a4n = 1.0;
}
else if (name == "Laguerre")
{
a1n = i - 1 + 1;
a2n = 2 * (i - 1) + 1;
a3n = -1.0;
a4n = i - 1;
}
else if (name == "Legendre")
{
a1n = i + 1 - 1;
a2n = 0.0;
a3n = 2.0 * (i - 1) + 1;
a4n = i - 1;
}
f2 = ((a2n + a3n * x) * f1 - a4n * f0) / a1n;
f0 = f1;
f1 = f2;
// deflate the function
for (int i = 0; i < m; i++)
{
double z = zeros[i];
double denom = x - z;
if (fabs(denom) < epsilon2)
zeros[i] = 1.001 * z;
else
f2 /= denom;
}
}
return f2;
}
double Mueller::g(std::string name,
std::string graph, double x,
int n, int m, std::vector<double>& zeros)
{
double f0 = 1, f1 = x, f2 = 1;
if (graph != "Extrema")
return f(name, x, n, m, zeros);
else
{
if (name == "Laguerre") {
if (x == 1.0)
x -= 0.0000001;
else if (x == -1.0)
x += 0.0000001;
}
double g0x = 1.0, g1x = 1.0, g2x = 1.0;
if (name == "Chebyshev") {
g0x = n;
g1x = -n * x;
g2x = 1.0 - x * x;
}
else if (name == "Laguerre") {
g0x = -n;
g1x = n;
g2x = x;
}
else if (name == "Legendre")
{
g0x = n;
g1x = -n * x;
g2x = 1.0 - x * x;
}
if (fabs(g2x) < epsilon2)
{
std::cout << "g2x is too small!";
return 0.0;
}
f0 = f(name, x, n - 1, 0, zeros);
f1 = f(name, x, n, 0, zeros);
f2 = (g1x * f1 + g0x * f0) / g2x;
// deflate the function
for (int i = 0; i < m; i++)
{
double z = zeros[i];
double denom = x - z;
if (fabs(denom) == 0.0)
zeros[i] = 1.001 * z;
else
f2 /= denom;
}
return f2;
}
}
double Mueller::OrthogonalPolynomialRoot(
std::string name,
std::string graph, int n,
double epsilon1, double epsilon2,
int maxIt, int m, double x0,
double x1, double x2, std::string& msg,
std::vector<double>& zeros)
{
double x3 = 0.0;
int i = 2;
while (true)
{
double h1 = x1 - x0;
if (fabs(h1) < epsilon2)
{
msg = "h1 is too small!";
return 0.0;
}
double h2 = x2 - x1;
if (fabs(h2) < epsilon2)
{
msg = "h2 is too small!";
return 0.0;
}
double f0 = g(name, graph, x0, n, m, zeros);
double f1 = g(name, graph, x1, n, m, zeros);
double f2 = g(name, graph, x2, n, m, zeros);
double g1 = (f1 - f0) / h1;
double g2 = (f2 - f1) / h2;
if (fabs(f2) < epsilon2)
{
x3 = x2;
break;
}
if (fabs(h1 + h2) < epsilon2)
{
msg = "h1 + h2 is too small!";
return 0.0;
}
double k2 = (g2 - g1) / (h1 + h2);
double c1 = g2 + h2 * k2;
double d = c1 * c1 - 4.0 * f2 * k2;
if (d < 0.0)
d = 0.0;
double s1 = sqrt(d);
double m1 = c1 - s1;
double p1 = c1 + s1;
double denom;
if (fabs(m1) > fabs(p1))
denom = m1;
else
denom = p1;
if (fabs(denom) < epsilon2)
{
msg = "denom is too small!";
return 0.0;
}
double h3 = -2.0 * f2 / denom;
x3 = x2 + h3;
double f3 = g(name, graph, x3, n, m, zeros);
if (fabs(f3) < epsilon2)
break;
if (fabs(h3) < epsilon2)
{
msg = "h3 is too small!";
return 0.0;
}
double g3 = (f3 - f2) / h3;
double delta = fabs(x3 - x2);
if (delta < epsilon1)
break;
i++;
if (i > maxIt)
break;
x0 = x1;
x1 = x2;
x2 = x3;
if (fabs(x1 - x0) < epsilon1)
break;
if (fabs(x2 - x1) < epsilon1)
break;
}
return x3;
}
void Mueller::OrthogonalPolynomialRoots(
std::string name,
std::string graph, int n, int nr,
unsigned int seed, std::string& msg,
std::vector<double>& zeros) {
double epsilon1 = 1.0e-10;
epsilon2 = 1.0e-20;
double x0 = -0.5, x1 = 0.0, x2 = 0.5;
double z;
int maxIt = 128;
zeros.resize(nr);
srand(seed);
for (int m = 0; m < nr; m++) {
x0 = (double)rand() / RAND_MAX;
x1 = (double)rand() / RAND_MAX;
x2 = (double)rand() / RAND_MAX;
z = zeros[m] = OrthogonalPolynomialRoot(
name, graph, n, epsilon1, epsilon2,
maxIt, m, x0, x1, x2, msg, zeros);
}
}
// Quadrature1d.cpp (c) Saturday, March 21, 2026
// by James Pate Williams, Jr., BA, BS, MSwE, PhD
// Reference: https://dlmf.nist.gov/3.5#x
#include <iomanip>
#include <iostream>
#include <vector>
#include "Mueller.h"
static void SelectionSort(
std::vector<double>& x,
std::vector<double>& w,
int n) {
for (int i = 0; i < n - 1; i++) {
for (int j = i + 1; j < n; j++) {
if (x[i] > x[j]) {
double t = x[i];
x[i] = x[j];
x[j] = t;
t = w[i];
w[i] = w[j];
w[j] = t;
}
}
}
}
static double fx(double x) {
return x * x * x * exp(-x);
}
static double GaussLegendre(
std::vector<double> x, std::vector<double> w,
double (*fx)(double), int n) {
double integral = 0.0;
for (int i = 0; i < n; i++) {
integral += fx(x[i]) * w[i];
}
return integral;
}
static double Integrate(double (*fx)(double), int n,
unsigned int seed) {
std::string msg;
std::vector<double> x(n);
std::vector<double> weights(n);
Mueller::OrthogonalPolynomialRoots(
"Legendre", "Zeros", n, n, seed, msg, x);
for (int i = 0; i < n; i++) {
double xi = x[i];
double dpn = Mueller::g("Legendre", "Extrema",
xi, n, 0, x);
weights[i] = 2.0 / ((1.0 - xi * xi) * dpn * dpn);
}
SelectionSort(x, weights, n);
std::cout << "Abscissas (x) and Weights (w)" << std::endl;
std::cout << "for Gauss-Legendre Quadrature";
std::cout << std::endl;
std::cout << "x\t\t\tw" << std::endl;
for (int i = 0; i < n; i++) {
std::cout << std::fixed << std::setprecision(15);
std::cout << std::setw(18);
std::cout << x[i] << '\t';
std::cout << std::fixed << std::setprecision(15);
std::cout << std::setw(18);
std::cout << weights[i] << std::endl;
}
return GaussLegendre(x, weights, fx, n);
}
static double SimpsonsRule(
int n, double a, double b,
double (*fx)(double))
{
double h = (b - a) / n;
double h2 = 2.0 * h;
double s = 0.0;
double t = 0.0;
double x = a + h;
for (int i = 1; i < n; i += 2)
{
s += fx(x);
x += h2;
}
x = a + h2;
for (int i = 2; i < n; i += 2)
{
t += fx(x);
x += h2;
}
return h * (fx(a) + 4 * s + 2 * t + fx(b)) / 3.0;
}
int main() {
unsigned int seed = 1;
for (int n = 10; n <= 40; n += 10)
{
double integral = Integrate(fx, n, seed++);
std::cout << std::fixed << std::setprecision(15);
std::cout << std::setw(18);
std::cout << n << '\t' << integral << std::endl;
}
for (int n = 50000; n <= 100000; n += 10000)
{
double integral = SimpsonsRule(n, -1, 1, fx);
std::cout << "Simpson ";
std::cout << n << "\t\t";
std::cout << std::fixed << std::setprecision(15);
std::cout << std::setw(18);
std::cout << integral << std::endl;
}
return 0;
}
Author: jamespatewilliamsjr
My whole legal name is James Pate Williams, Jr. I was born in LaGrange, Georgia approximately 70 years ago. I barely graduated from LaGrange High School with low marks in June 1971. Later in June 1979, I graduated from LaGrange College with a Bachelor of Arts in Chemistry with a little over a 3 out 4 Grade Point Average (GPA). In the Spring Quarter of 1978, I taught myself how to program a Texas Instruments desktop programmable calculator and in the Summer Quarter of 1978 I taught myself Dayton BASIC (Beginner's All-purpose Symbolic Instruction Code) on LaGrange College's Data General Eclipse minicomputer. I took courses in BASIC in the Fall Quarter of 1978 and FORTRAN IV (Formula Translator IV) in the Winter Quarter of 1979. Professor Kenneth Cooper, a genius poly-scientist taught me a course in the Intel 8085 microprocessor architecture and assembly and machine language. We would hand assemble our programs and insert the resulting machine code into our crude wooden box computer which was designed and built by Professor Cooper. From 1990 to 1994 I earned a Bachelor of Science in Computer Science from LaGrange College. I had a 4 out of 4 GPA in the period 1990 to 1994. I took courses in C, COBOL, and Pascal during my BS work. After graduating from LaGrange College a second time in May 1994, I taught myself C++. In December 1995, I started using the Internet and taught myself client-server programming. I created a website in 1997 which had C and C# implementations of algorithms from the "Handbook of Applied Cryptography" by Alfred J. Menezes, et. al., and some other cryptography and number theory textbooks and treatises.
View all posts by jamespatewilliamsjr
Numerical Explorations 