Tag: coding
Blog Entry © Saturday, May 16, 2026, by James Pate Williams, Jr. Some More Linear Algebra Examples
Blog Entry © Thursday, May 14, 2026, by James Pate Williams, Jr. More Numerical Integration Results
// NumericalIntegrals.cpp (c) Thursday, May 14, 2026
// by James Pate Williams, Jr., BA, BS, MSwE, PhD
#include <iomanip>
#include <iostream>
#include <vector>
#include <stdlib.h>
static double f(double x) {
return sin(x);
}
static double MonteCarlo(double a, double b,
double (*f)(double), int n){
double sum = 0;
for (int i = 0; i < n; i++) {
double x = (b - a) * (double)rand() / RAND_MAX;
sum += f(x);
}
return (b - a) * sum / n;
}
static double Factorial(int n) {
double factorial = 1.0;
for (int i = 2; i <= n; i++)
factorial *= i;
return factorial;
}
static double Series(double a, double b, int n)
{
double sumA = 0.0, sumB = 0.0;
int sign = 1;
for (int i = 0; i <= n; i++) {
sumA += sign * pow(a, 2 * i + 2) /
Factorial(2 * i + 2);
sign *= -1;
}
sign = 1;
for (int i = 0; i <= n; i++) {
sumB += sign * pow(b, 2 * i + 2) /
Factorial(2 * i + 2);
sign *= -1;
}
return sumB - sumA;
}
static double CompositeTrapezoidalRule(
double a, double b, int n) {
double pi = 4.0 * atan(1.0);
double endPts = 0.5 * (f(a) + f(b));
double sum = 0, xk = 0.0;
double h = (b - a) / n;
for (int k = 1; k <= n - 1; k++) {
xk = a + k * h;
sum += f(xk);
}
return h * (0.5 * endPts + sum);
}
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;
}
static void Romberg(double a, double b,
double (*f)(double), int mStart, int nRow,
std::vector<std::vector<double>>& T) {
int m = mStart;
double h = (b - a) / m;
double sum = 0.5 * (f(a) + f(b));
if (m > 1) {
for (int i = 1; i <= m - 1; i++) {
sum += f(a + i * h);
}
}
T[0][0] = sum * h;
std::cout << "romberg t-table" << std::endl;
std::cout << std::fixed;
std::cout << std::setprecision(5) << T[0][0];
std::cout << std::endl;
if (nRow < 2)
return;
for (int k = 2; k <= nRow; k++) {
h /= 2.0;
m *= 2;
sum = 0.0;
for (int i = 1; i <= m; i += 2) {
sum += f(a + i * h);
}
T[k][1] = 0.5 * T[k - 1LL][1] + sum * h;
for (int j = 1; j <= k - 1; j++) {
T[k - 1LL][j] = T[k][j] - T[k - 1LL][j];
T[k][j + 1LL] = T[k][j] + T[k - 1LL][j] /
(pow(4.0, j) - 1.0);
}
for (int j = 1; j <= k; j++) {
std::cout << std::fixed;
std::cout << std::setprecision(5);
std::cout << T[k][j] << '\t';
}
std::cout << std::endl;
}
if (nRow < 3) {
return;
}
std::cout << "table of ratios" << std::endl;
double ratio = 0.0;
for (int k = 1; k <= nRow - 2; k++) {
for (int j = 1; j <= k; j++) {
if (T[k + 1LL][j] == 0.0) {
ratio = 0.0;
}
else {
ratio = T[k][j] / T[k + 1LL][j];
}
T[k][j] = ratio;
}
for (int j = 1; j <= k; j++) {
std::cout << std::fixed;
std::cout << std::setprecision(5);
std::cout << T[k][j] << '\t';
}
std::cout << std::endl;
}
}
double MonteCarloVolume(double R, int n)
{
double pi = 4.0 * atan(1.0), pi2 = 2.0 * pi;
double R2 = R * R, sum = 0;
for (int i = 0; i < n; i++)
{
double r = R2 * (double)rand() / RAND_MAX;
double t = pi * (double)rand() / RAND_MAX;
double p = pi2 * (double)rand() / RAND_MAX;
sum += r * r * sin(t);
}
return R * pi * pi2 * sum / n;
}
int main()
{
srand(1);
std::vector<std::vector<double>> T;
T.resize(35);
for (int i = 0; i < 35; i++) {
T[i].resize(35);
}
Romberg(0.0, 2.0, f, 2, 7, T);
std::cout << std::setprecision(11);
std::cout << "analytic integral of sine = " << -cos(2.0) + cos(0.0);
std::cout << std::endl;
std::cout << "simpson's rule integral = " << SimpsonsRule(500, 0, 2.0, f);
std::cout << std::endl;
std::cout << "monte carlo integral = " << MonteCarlo(0.0, 2.0, f, 2130);
std::cout << std::endl;
std::cout << "infinite series integral = " << Series(0.0, 2.0, 16);
std::cout << std::endl;
double integral = CompositeTrapezoidalRule(0.0, 2.0, 175000000);
std::cout << "romberg integral = " << integral << std::endl;
std::cout << "actual spherical volume = " << 4.0 * 4.0 * atan(1.0) / 3.0;
std::cout << std::endl;
double volume = MonteCarloVolume(1.0, 1000000);
std::cout << "approx spherical volume = " << volume;
std::cout << std::endl;
}
Blog Entry © Wednesday, May 13, 2026, by James Pate Williams, Jr. Adaptive n-Quadrature Versus Monte Carlo Integration
Blog Entry © Monday, May 11, 2026, by James Pate Williams, Jr., Laplace Equation in a Solid Cylinder
Blog Entry © Tuesday, May 5, 2026, by James Pate Williams, Jr., Solution of the Laplace Equation on a Rectangle
Blog Entry © Sunday, April 19, 2026, by James Pate Williams, Jr., Scattering from a Spherically Symmetric Potential
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;
}
Numerical Explorations 