Unfortunately, I do not own an E-Bow anymore.
Blog Entry © Monday, April 20, 2026, by James Pate Williams, Jr., Vector Analysis Continued and Perhaps Corrected
Blog Entry © Sunday, April 19, 2026, by James Pate Williams, Jr., Scattering from a Spherically Symmetric Potential
Capo 2nd Joseph (Joe) Gay Lead Guitar and James Pate Williams Jr Rythm Guitar a la May 26, 2009 or Earlier
I used a Fender 12-string acoustic electric guitar.
Vector Analysis by James Pate Williams, Jr. Exercises and Supplementary Problems from Introduction to Vector Analysis Fourth Edition© 1979 by Harry F. Davis and Arthur David Snider Selected Exercises from Chapter 1 Pages 48 and 50
// VectorAnalysis.cpp © Tuesday, April 14, 2026
// by James Pate Williams, Jr.
// Reference: Introduction to Vector Analysis Fourth Edition
// © 1979 by Harry F. Davis and Arthur David Snider
#include <iostream>
#include <vector>
static double InnerProduct(
std::vector<double> A,
std::vector<double> B,
int n)
{
double sum = 0.0;
for (int i = 0; i < n; i++)
sum += A[i] * B[i];
return sum;
}
static void VectorProduct(
std::vector<double> A,
std::vector<double> B,
std::vector<double>&C)
{
C.resize(3);
C[0] = A[1] * B[2] - A[2] * B[1];
C[1] = A[2] * B[0] - A[0] * B[2];
C[2] = A[0] * B[1] - A[1] * B[0];
}
static double TripleProduct(
std::vector<double> A,
std::vector<double> B,
std::vector<double> C)
{
double sum0 = 0.0, sum1 = 0.0;
sum0 += A[0] * B[1] * C[2] + A[1] * B[2] * C[0] + A[2] * B[0] * C[1];
sum1 += A[1] * B[0] * C[2] + A[0] * B[2] * C[1] + A[2] * B[1] * C[0];
return sum0 - sum1;
}
static void Exercises_Section_1_13_1_Triple_Products()
{
// triple products (a)
std::vector<double> A1 = { 2, 0, 0 };
std::vector<double> B1 = { 0, 3, 0 };
std::vector<double> C1 = { 0, 0, 5 };
double tp_1 = TripleProduct(A1, B1, C1);
std::cout << "Exercise 1 (a)" << '\t';
std::cout << tp_1 << std::endl;
// triple products (b)
std::vector<double> A2 = { 1, 1, 1 };
std::vector<double> B2 = { 3, 1, 0 };
std::vector<double> C2 = { 0, -1, 5 };
std::cout << "Exercise 1 (b)" << '\t';
double tp_2 = TripleProduct(A2, B2, C2);
std::cout << tp_2 << std::endl;
// triple products (c)
std::vector<double> A3 = { 2, -1, 1 };
std::vector<double> B3 = { 1, 1, 1 };
std::vector<double> C3 = { 2, 0, 3 };
double tp_3 = TripleProduct(A3, B3, C3);
std::cout << "Exercise 1 (c)" << '\t';
std::cout << tp_3 << std::endl;
// triple products (d)
std::vector<double> A4 = { 0, 0, 1 };
std::vector<double> B4 = { 1, 0, 0 };
std::vector<double> C4 = { 0, 1, 0 };
double tp_4 = TripleProduct(A4, B4, C4);
std::cout << "Exercise 1 (d)" << '\t';
std::cout << tp_4 << std::endl;
// volume of a parallelpipped
std::vector<double> A5 = { 3, 4, 1 };
std::vector<double> B5 = { 2, 3, 4 };
std::vector<double> C5 = { 0, 0, 5 };
double tp_5 = TripleProduct(A5, B5, C5);
std::cout << "Exercise 2" << '\t';
std::cout << tp_5 << std::endl;
// volume of a parallelpipped
std::vector<double> A6 = { 3, 2, 1 };
std::vector<double> B6 = { 4, 2, 1 };
std::vector<double> C6 = { 0, 1, 4 };
std::vector<double> D6 = { 0, 0, 7 };
std::vector<double> AB = { -1, 0, 0 };
std::vector<double> AC = { 3, 1, -3 };
std::vector<double> AD = { 3, 2, -6 };
std::cout << "Exercise 3" << '\t';
double tp_6 = TripleProduct(AB, AC, AD);
std::cout << tp_6 << std::endl;
// volume of a tetrahedron
std::vector<double> AB1 = { 1, 1, 0 };
std::vector<double> AC1 = { 1, -1, 0 };
std::vector<double> AD1 = { 0, 0, 2 };
std::cout << "Exercise 4" << '\t';
double tp_7 = TripleProduct(AB1, AC1, AD1) / 6.0;
std::cout << fabs(tp_7) << std::endl;
std::vector<double> P1 = { 0, 0, 0 };
std::vector<double> P2 = { 1, 1, 0 };
std::vector<double> P3 = { 3, 4, 0 };
std::vector<double> P4 = { 4, 5, 0 };
std::vector<double> P5 = { 0, 0, 1 };
std::vector<double> Q1 = { 1, 1, 0 };
std::vector<double> Q2 = { 2, 3, 0 };
std::vector<double> Q3 = { 1, 1, 1 };
std::cout << "Exercise 5" << '\t';
double tp_8 = TripleProduct(Q1, Q2, Q3);
std::cout << fabs(tp_8) << std::endl;
std::vector<double> A10 = { 1, 1, 1 };
std::vector<double> B10 = { 2, 4, -1 };
std::vector<double> C10 = { 1, 1, 3 };
double tp_10 = TripleProduct(A10, B10, C10);
std::vector<double> D10 = { 0 };
VectorProduct(A10, B10, D10);
double magnitude = sqrt(InnerProduct(D10, D10, 3));
std::cout << "Exercise 10" << '\t';
std::cout << tp_10 / magnitude << '\t';
std::cout << 2.0 * sqrt(38.0) / 19.0 << std::endl;
std::vector<double> A11 = { 1, 1, 1 };
std::vector<double> B11 = { 2, 4, -1 };
std::vector<double> C11 = { 1, 1, 3 };
std::vector<double> D11 = { 3, 2, 1 };
std::vector<double> AB11(3), BC11(3), CA11(3), BCCA11(3);
VectorProduct(A11, B11, AB11);
VectorProduct(B11, C11, BC11);
VectorProduct(C11, A11, CA11);
VectorProduct(BC11, CA11, BCCA11);
double Q11 = InnerProduct(AB11, BCCA11, 3);
double A_x = A11[0], A_y = A11[1], A_z = A11[2];
double B_x = B11[0], B_y = B11[1], B_z = B11[2];
double C_x = C11[0], C_y = C11[1], C_z = C11[2];
double term1 = +(A_y * B_z - A_z * B_y) * (B_z * C_x - B_x * C_z) * (C_x * A_y - C_y * A_x);
double term2 = -(A_y * B_z - A_z * B_y) * (B_x * C_y - B_y * C_x) * (C_z * A_x - C_x * A_z);
double term3 = +(A_z * B_x - A_x * B_z) * (B_x * C_y - B_y * C_x) * (C_y * A_z - C_z * A_y);
double term4 = -(A_z * B_x - A_x * B_z) * (B_y * C_z - B_z * C_y) * (C_x * A_y - C_y * A_x);
double term5 = +(A_x * B_y - A_y * B_x) * (B_y * C_z - B_z * C_y) * (C_z * A_x - C_x * A_z);
double term6 = -(A_x * B_y - A_y * B_x) * (B_z * C_x - B_x * C_z) * (C_y * A_z - C_z * A_y);
double P11 = term1 + term2 + term3 + term4 + term5 + term6;
std::cout << "Q = (A x B) . (B x C) x (C x A) = " << Q11 << std::endl;
std::cout << "P = (A x B) . (B x C) x (C x A) = " << P11 << std::endl;
}
static void Exercises_Section_1_14_Vector_Identities()
{
std::vector<double> A11 = { 1, 1, 1 };
std::vector<double> B11 = { 2, 4, -1 };
std::vector<double> C11 = { 1, 1, 3 };
std::vector<double> D11 = { 3, 2, 1 };
std::cout << "Section 1.14 page 50 Exercises Exercise 1" << std::endl;
std::cout << "A = " << A11[0] << '\t' << A11[1] << '\t' << A11[2] << std::endl;
std::cout << "B = " << B11[0] << '\t' << B11[1] << '\t' << B11[2] << std::endl;
std::cout << "C = " << C11[0] << '\t' << C11[1] << '\t' << C11[2] << std::endl;
std::cout << "D = " << D11[0] << '\t' << D11[1] << '\t' << D11[2] << std::endl;
std::cout << "TPI1 = (A x B) x (C x D) = [A, B, D]C - [A, B, C]D = " << std::endl;
double TP1411a = TripleProduct(A11, B11, D11);
double TP1411b = TripleProduct(A11, B11, C11);
std::cout << "[A, B, D] = " << TP1411a << std::endl;
std::cout << "[A, B, C] = " << TP1411b << std::endl;
std::cout << "TPI1_x = " << TP1411a * C11[0] << std::endl;
std::cout << "TPI1_y = " << TP1411a * C11[1] << std::endl;
std::cout << "TPI1_z = " << TP1411a * C11[2] << std::endl;
std::cout << "TPI2_x = " << TP1411b * D11[0] << std::endl;
std::cout << "TPI2_y = " << TP1411b * D11[1] << std::endl;
std::cout << "TPI2_z = " << TP1411b * D11[2] << std::endl;
std::cout << "RHS1 = [A, B, D]C - [A, B, C]D = " << std::endl;
std::vector<double> RHS1(3);
RHS1[0] = TP1411a * C11[0] - TP1411b * D11[0];
RHS1[1] = TP1411a * C11[1] - TP1411b * D11[1];
RHS1[2] = TP1411a * C11[2] - TP1411b * D11[2];
std::cout << "RHS1_x = " << RHS1[0] << std::endl;
std::cout << "RHS1_y = " << RHS1[1] << std::endl;
std::cout << "RHS1_z = " << RHS1[2] << std::endl;
std::vector<double> CD11(3), TD11(3);
std::vector<double> AB11(3), BC11(3), CA11(3), BCCA11(3);
VectorProduct(A11, B11, AB11);
VectorProduct(B11, C11, BC11);
VectorProduct(C11, A11, CA11);
VectorProduct(BC11, CA11, BCCA11);
VectorProduct(A11, B11, AB11);
VectorProduct(C11, D11, CD11);
VectorProduct(AB11, CD11, TD11);
std::cout << "A = " << A11[0] << '\t' << A11[1] << '\t' << A11[2] << std::endl;
std::cout << "B = " << B11[0] << '\t' << B11[1] << '\t' << B11[2] << std::endl;
std::cout << "C = " << C11[0] << '\t' << C11[1] << '\t' << C11[2] << std::endl;
std::cout << "D = " << D11[0] << '\t' << D11[1] << '\t' << D11[2] << std::endl;
std::cout << "A x B = " << AB11[0] << '\t' << AB11[1] << '\t' << AB11[2] << std::endl;
std::cout << "C x D = " << CD11[0] << '\t' << CD11[1] << '\t' << CD11[2] << std::endl;
std::cout << "TD11 = (A x B) x (C x D) = " << std::endl;
std::cout << "TD11_x = " << TD11[0] << std::endl;
std::cout << "TD11_y = " << TD11[1] << std::endl;
std::cout << "TD11_z = " << TD11[2] << std::endl;
VectorProduct(B11, C11, BC11);
VectorProduct(C11, A11, CA11);
VectorProduct(BC11, CA11, D11);
VectorProduct(A11, B11, AB11);
double ip12 = InnerProduct(AB11, D11, 3);
std::cout << "2. Inner Product = " << ip12 << std::endl;
double tp12 = TripleProduct(A11, B11, C11);
std::cout << "2. Triple Product ^ 2 = " << tp12 * tp12 << std::endl;
std::vector<double> ABC11(3), BAC11(3), CAB11(3);
VectorProduct(A11, BC11, ABC11);
VectorProduct(B11, CA11, BAC11);
VectorProduct(C11, AB11, CAB11);
double zx = ABC11[0] + BAC11[0] + CAB11[0];
double zy = ABC11[1] + BAC11[1] + CAB11[1];
double zz = ABC11[2] + BAC11[2] + CAB11[2];
std::cout << "3. Zero Vector = " << zx << ' ' << zy << ' ' << zz;
std::cout << std::endl;
}
int main()
{
Exercises_Section_1_13_1_Triple_Products();
Exercises_Section_1_14_Vector_Identities();
return 0;
}
Blog Entry © Tuesday, April 14, 2026, by James Pate Williams, Jr. Exercises and Supplementary Problems from Introduction to Vector Analysis Fourth Edition © 1979 by Harry F. Davis and Arthur David Snider Selected Exercises from Chapter 1 Page 48
// VectorAnalysis.cpp © Tuesday, April 14, 2026
// by James Pate Williams, Jr.
// Reference: Introduction to Vector Analysis Fourth Edition
// © 1979 by Harry F. Davis and Arthur David Snider
#include <iostream>
#include <vector>
static double InnerProduct(
std::vector<double> A,
std::vector<double> B,
int n)
{
double sum = 0.0;
for (int i = 0; i < n; i++)
sum += A[i] * B[i];
return sum;
}
static void VectorProduct(
std::vector<double> A,
std::vector<double> B,
std::vector<double>&C)
{
C.resize(3);
C[0] = A[1] * B[2] - A[2] * B[1];
C[1] = A[0] * B[2] - A[2] * B[0];
C[2] = A[0] * B[1] - A[1] * B[0];
}
static double TripleProduct(
std::vector<double> A,
std::vector<double> B,
std::vector<double> C)
{
double sum0 = 0.0, sum1 = 0.0;
sum0 += A[0] * B[1] * C[2] + A[1] * B[2] * C[0] + A[2] * B[0] * C[1];
sum1 += A[1] * B[0] * C[2] + A[0] * B[2] * C[1] + A[2] * B[1] * C[0];
return sum0 - sum1;
}
static void Exercises_Section_1_13_1_Triple_Products()
{
// triple products (a)
std::vector<double> A1 = { 2, 0, 0 };
std::vector<double> B1 = { 0, 3, 0 };
std::vector<double> C1 = { 0, 0, 5 };
double tp_1 = TripleProduct(A1, B1, C1);
std::cout << "Exercise 1 (a)" << '\t';
std::cout << tp_1 << std::endl;
// triple products (b)
std::vector<double> A2 = { 1, 1, 1 };
std::vector<double> B2 = { 3, 1, 0 };
std::vector<double> C2 = { 0, -1, 5 };
std::cout << "Exercise 1 (b)" << '\t';
double tp_2 = TripleProduct(A2, B2, C2);
std::cout << tp_2 << std::endl;
// triple products (c)
std::vector<double> A3 = { 2, -1, 1 };
std::vector<double> B3 = { 1, 1, 1 };
std::vector<double> C3 = { 2, 0, 3 };
double tp_3 = TripleProduct(A3, B3, C3);
std::cout << "Exercise 1 (c)" << '\t';
std::cout << tp_3 << std::endl;
// triple products (d)
std::vector<double> A4 = { 0, 0, 1 };
std::vector<double> B4 = { 1, 0, 0 };
std::vector<double> C4 = { 0, 1, 0 };
double tp_4 = TripleProduct(A4, B4, C4);
std::cout << "Exercise 1 (d)" << '\t';
std::cout << tp_4 << std::endl;
// volume of a parallelpipped
std::vector<double> A5 = { 3, 4, 1 };
std::vector<double> B5 = { 2, 3, 4 };
std::vector<double> C5 = { 0, 0, 5 };
double tp_5 = TripleProduct(A5, B5, C5);
std::cout << "Exercise 2" << '\t';
std::cout << tp_5 << std::endl;
// volume of a parallelpipped
std::vector<double> A6 = { 3, 2, 1 };
std::vector<double> B6 = { 4, 2, 1 };
std::vector<double> C6 = { 0, 1, 4 };
std::vector<double> D6 = { 0, 0, 7 };
std::vector<double> AB = { -1, 0, 0 };
std::vector<double> AC = { 3, 1, -3 };
std::vector<double> AD = { 3, 2, -6 };
std::cout << "Exercise 3" << '\t';
double tp_6 = TripleProduct(AB, AC, AD);
std::cout << tp_6 << std::endl;
// volume of a tetrahedron
std::vector<double> AB1 = { 1, 1, 0 };
std::vector<double> AC1 = { 1, -1, 0 };
std::vector<double> AD1 = { 0, 0, 2 };
std::cout << "Exercise 4" << '\t';
double tp_7 = TripleProduct(AB1, AC1, AD1) / 6.0;
std::cout << fabs(tp_7) << std::endl;
std::vector<double> P1 = { 0, 0, 0 };
std::vector<double> P2 = { 1, 1, 0 };
std::vector<double> P3 = { 3, 4, 0 };
std::vector<double> P4 = { 4, 5, 0 };
std::vector<double> P5 = { 0, 0, 1 };
std::vector<double> Q1 = { 1, 1, 0 };
std::vector<double> Q2 = { 2, 3, 0 };
std::vector<double> Q3 = { 1, 1, 1 };
std::cout << "Exercise 5" << '\t';
double tp_8 = TripleProduct(Q1, Q2, Q3);
std::cout << fabs(tp_8) << std::endl;
std::vector<double> A10 = { 1, 1, 1 };
std::vector<double> B10 = { 2, 4, -1 };
std::vector<double> C10 = { 1, 1, 3 };
double tp_10 = TripleProduct(A10, B10, C10);
std::vector<double> D10 = { 0 };
VectorProduct(A10, B10, D10);
double magnitude = sqrt(InnerProduct(D10, D10, 3));
std::cout << "Exercise 10" << '\t';
std::cout << tp_10 / magnitude << '\t';
std::cout << 2.0 * sqrt(38.0) / 19.0 << std::endl;
}
int main()
{
Exercises_Section_1_13_1_Triple_Products();
return 0;
}
Blog Entry © Tuesday, April 7, 2026, by James Pate Williams, Jr., Hydrogen-like Atom Polar and Azimuthal Wavefunctions
Blog Entry © Saturday April 4, 2026, by James Pate Williams, Jr., Hydrogen-like Radial Electron Distribution Functions in CPP
#pragma once
class RadialWaveFunction
{
private:
static double Factorial(int n);
static double Laguerre(double rho, int n, int l);
public:
static double R(double r, int Z, int n, int l);
};
#include <math.h>
#include "RadialWaveFunction.h"
double RadialWaveFunction::Factorial(int n)
{
double factorial = 1.0;
for (int i = 2; i <= n; i++)
factorial *= i;
return factorial;
}
double RadialWaveFunction::Laguerre(double rho, int n, int l)
{
double sum = 0.0;
for (int k = 0; k <= n - l - 1; k++)
{
double factor1 = pow(-1, k + 2 * l + 1);
double factor2 = pow(Factorial(n + l), 2.0);
double factor3 = pow(rho, k);
double numer = factor1 * factor2 * factor3;
double factor4 = Factorial(n - l - 1 - k);
double factor5 = Factorial(2 * l + 1 + k);
double factor6 = Factorial(k);
double denom = factor4 * factor5 * factor6;
sum += numer / denom;
}
return sum;
}
double RadialWaveFunction::R(double r, int Z, int n, int l)
{
double rho = 2.0 * Z * r / n;
double numer1 = pow(2.0 * Z / n, 3.0);
double numer2 = Factorial(n - l - 1);
double denom1 = Factorial(n + n);
double denom2 = pow(Factorial(n + l), 3.0);
double numer3 = -sqrt(numer1 * numer2 /
(denom1 * denom2));
double exp2 = exp(-0.5 * rho);
double rhol = pow(rho, l);
return numer3 * exp2 * rhol * Laguerre(rho, n, l);
}
Blog Entry © Thursday, April 2, 2026, by James Pate Williams, Jr., Matrix Inverses that Exist and Their Characteristic Polynomials
#pragma once
#include <vector>
class Inverse
{
public:
static bool Compute(
int n,
std::vector<std::vector<double>>& M,
std::vector<std::vector<double>>& X);
static void CharacteristicPolynomial(
int n, std::vector<double>& a,
std::vector<std::vector<double>>& C,
std::vector<std::vector<double>>& M,
std::vector<std::vector<double>>& Madj);
};
// Reference: 5.2: The Characteristic Polynomial - Mathematics LibreTexts
#include "Inverse.h"
#include <vector>
bool Inverse::Compute(
int n,
std::vector<std::vector<double>>& M,
std::vector<std::vector<double>>& X) {
double d = 0.0;
int i = 0, j = 0, k = 0, l = 0;
std::vector<double> C(n);
std::vector<std::vector<double>> B;
B.resize(n);
X.resize(n);
for (i = 0; i < n; i++)
{
B[i].resize(n);
X[i].resize(n);
for (j = 0; j < n; j++)
{
B[i][j] = 0.0;
X[i][j] = 0.0;
}
B[i][i] = 1.0;
}
j = -1;
Step2:
j++;
if (j == n)
goto Step6;
for (i = j; i < n; i++) {
if (M[i][j] != 0)
break;
}
if (i == n)
return false;
if (i > j) {
for (l = j; l < n; l++) {
double temp = M[i][l];
M[i][l] = M[j][l];
M[j][l] = temp;
}
for (l = 0; l < n; l++) {
double temp = B[i][l];
B[i][l] = B[j][l];
B[j][l] = temp;
}
}
d = 1.0 / M[j][j];
for (k = j + 1; k < n; k++) {
C[k] = d * M[k][j];
}
for (k = j + 1; k < n; k++) {
for (l = j + 1; l < n; l++) {
M[k][l] -= C[k] * M[j][l];
}
}
for (k = j + 1; k < n; k++) {
for (l = 0; l < n; l++) {
B[k][l] -= C[k] * B[j][l];
}
}
goto Step2;
Step6:
for (i = n - 1; i >= 0; i--) {
std::vector<double> sum(n, 0);
for (j = i + 1; j < n; j++) {
for (k = 0; k < n; k++) {
sum[k] += M[i][j] * X[j][k];
}
}
for (j = 0; j < n; j++)
X[i][j] = (B[i][j] - sum[j]) / M[i][i];
}
return true;
}
double delta(int i, int j) {
return i == j ? 1.0 : 0.0;
}
void Inverse::CharacteristicPolynomial(
int n, std::vector<double>& a,
std::vector<std::vector<double>>& C,
std::vector<std::vector<double>>& M,
std::vector<std::vector<double>>& Madj) {
int i = 0, j = 0, k = 0, l = 0;
std::vector<double> X(n);
std::vector<std::vector<double>> MC;
a.resize(n + 1LL);
C.resize(n);
M.resize(n);
MC.resize(n);
Madj.resize(n);
for (i = 0; i < n; i++) {
C[i].resize(n);
M[i].resize(n);
MC[i].resize(n);
Madj[i].resize(n);
for (j = 0; j < n; j++) {
C[i][j] = delta(i, j);
}
}
a[0] = 1.0;
i = 0;
Step2:
i++;
if (i == n) {
for (j = 0; j < n; j++) {
for (k = 0; k < n; k++) {
double sum = 0.0;
for (l = 0; l < n; l++) {
sum += M[j][l] * C[l][k];
}
MC[j][k] = sum;
}
}
double trace = 0.0;
for (j = 0; j < n; j++) {
trace += MC[j][j];
}
a[n] = -trace / n;
for (j = 0; j < n; j++) {
for (k = 0; k < n; k++) {
Madj[j][k] = pow(-1, n - 1) * C[j][k];
}
}
return;
}
for (j = 0; j < n; j++) {
for (k = 0; k < n; k++) {
double sum = 0.0;
for (l = 0; l < n; l++) {
sum += M[j][l] * C[l][k];
}
MC[j][k] = sum;
}
}
for (j = 0; j < n; j++) {
for (k = 0; k < n; k++) {
C[j][k] = MC[j][k];
}
}
double tr = 0.0;
for (j = 0; j < n; j++) {
tr += C[j][j];
}
a[i] = -tr / i;
for (j = 0; j < n; j++) {
for (k = 0; k < n; k++) {
C[j][k] += a[i] * delta(j, k);
}
}
goto Step2;
}
Blog Entry © April 1, 2026, by James Pate Williams, Jr. More Romberg Integration Results
/*
Author: Pate Willimas c 1995
The following program is a translation of the FORTRAN77
subroutine found in "Elementary Numerical Analysis
An Algorithmic Approach Third Edition" (c) 1980 by S. D.
Conte and Carl de Boor pages 343-344. The program
uses Romberg extrapolation to find the integral of a
function. Also see: https://dlmf.nist.gov/3.5#E10
*/
#include <math.h>
#include <stdio.h>
#include <vector>
#include "Bessel.h"
typedef double real;
long fe;
real t[10][10];
real f(real x)
{
fe++;
return(expl(-x * x));
}
real g(real x)
{
fe++;
std::vector<real> j(4);
Bessel::bessj(x, 0, j);
return exp(-x) * j[0];
}
real romberg(
real a, real b,
int start, int row,
real (*fx)(real))
{
int i, /*j, */k, m;
real h, ratio, sum;
m = start;
h = (b - a) / m;
sum = 0.5 * (fx(a) + fx(b));
if (m > 1)
for (i = 1; i < m; i++)
sum += fx(a + i * h);
t[1][1] = h * sum;
//printf_s("Romberg T-Table\n%9.7lf\n", t[1][1]);
if (row < 2) return(t[1][1]);
for (k = 2; k <= row; k++)
{
h = 0.5 * h;
m *= 2;
sum = 0.0;
for (i = 1; i <= m; i += 2)
sum += fx(a + h * i);
t[k][1] = 0.5 * t[k - 1][1] + sum * h;
for (i = 1; i < k; i++)
{
t[k - 1][i] = t[k][i] - t[k - 1][i];
t[k][i + 1] = t[k][i] - t[k - 1][i] /
(powl(4.0, (real)i) - 1.0);
//for (j = 1; j < k; j++)
//printf_s("%9.7lf ", t[k][j]);
}
//printf_s("\n");
}
if (row < 3) return(t[2][2]);
//printf_s("Table of ratios\n");
for (k = 1; k <= row - 2; k++) {
for (i = 1; i <= k; i++)
{
if (t[k + 1][i] == 0.0)
ratio = 0.0;
else
ratio = t[k][i] / t[k + 1][i];
t[k][i] = ratio;
//printf_s("%5.2lf ", ratio);
}
//printf_s("\n");
}
return(t[row][row - 1]);
}
real CompositeTrapezoidalRule(
real a, real b, int n,
real(*fx)(real)) {
real endPts = 0.5 * (fx(a) + fx(b));
real sum = 0, xk = 0.0;
real h = (b - a) / n;
for (int k = 1; k <= n - 1; k++)
{
xk = a + k * h;
sum += fx(xk);
}
return h * (endPts + sum);
}
double CompositeSimpsonsRule(
double a, double b, int n,
real(*fx)(real))
{
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;
}
void Integrate(int trial, real(*fx)(real))
{
real a = 0.0, b = 0.0;
int row, start;
if (trial == 0)
{
a = 0.0;
b = 1.0;
}
else
{
a = 0.0;
b = 30.0;
}
if (trial == 0)
printf("Romberg integration of f(x) = exp(- x * x)\n");
else
printf("Romberg integration of g(x) = exp(-x) * J0(x)\n");
printf("number of trapezoidal intervals = ");
scanf_s("%d", &start);
printf("number of rows in table (<= 8) = ");
scanf_s("%d", &row);
fe = 0;
printf_s("integral = %14.11lf\tevals = %d\n",
romberg(a, b, start, row, fx), fe);
printf("number of trapezoidal intervals = ");
scanf_s("%d", &row);
fe = 0;
printf_s("integral = %14.11lf\tevals = %d\n",
CompositeTrapezoidalRule(a, b, row, fx), fe);
fe = 0;
printf("number of Simpson's intervals = ");
scanf_s("%d", &row);
printf_s("integral = %14.11lf\tevals = %d\n",
CompositeSimpsonsRule(a, b, row, fx), fe);
if (trial == 0)
printf_s("integral = %14.11lf\n", 0.74682413279);
else
{
real pi = 4.0 * atan(1.0);
real si = sin(pi * 45.0 / 180.0);
printf_s("integral = %14.11lf\n", si);
}
}
int main(void)
{
Integrate(0, f);
Integrate(1, g);
return(0);
}
Numerical Explorations 