Elementary Physics – Page 3 – Numerical Explorations

Exercises Szabo Ostlund Download

#include <complex>
#include <vector>

class CmpLinearAlgebra
{
public:
    static void CmpPrintMatrix(
        int m, int n,
        std::vector<std::vector<std::complex<double>>>& Ac);
    static void CmpAddition(
        size_t m, size_t n,
        std::vector<std::vector<std::complex<double>>>& A,
        std::vector<std::vector<std::complex<double>>>& B,
        std::vector<std::vector<std::complex<double>>>& C);
    static void CmpSubtraction(
        size_t m, size_t n,
        std::vector<std::vector<std::complex<double>>>& A,
        std::vector<std::vector<std::complex<double>>>& B,
        std::vector<std::vector<std::complex<double>>>& C);
    static void CmpMultiply(
        size_t m, size_t n, size_t p,
        std::vector<std::vector<std::complex<double>>>& A,
        std::vector<std::vector<std::complex<double>>>& B,
        std::vector<std::vector<std::complex<double>>>& C);
    static void CmpAnticommutator(
        size_t n,
        std::vector<std::vector<std::complex<double>>>& A,
        std::vector<std::vector<std::complex<double>>>& B,
        std::vector<std::vector<std::complex<double>>>& C);
    static void CmpCommutator(
        size_t n,
        std::vector<std::vector<std::complex<double>>>& A,
        std::vector<std::vector<std::complex<double>>>& B,
        std::vector<std::vector<std::complex<double>>>& C);
    static void CmpAdjoint(
        size_t m, size_t n,
        std::vector<std::vector<std::complex<double>>>& Ac,
        std::vector<std::vector<std::complex<double>>>& Ad);
    static std::complex<double> CmpDeterminant(
        bool& failure, int n,
        std::vector<std::vector<std::complex<double>>>& A);
    static bool CmpGaussianElimination(
        int m, int n,
        std::vector<std::vector<std::complex<double>>>& A,
        std::vector<std::complex<double>>& b,
        std::vector<std::complex<double>>& x,
        std::vector<size_t>& pivot);
    static bool CmpGaussianFactor(
        int n, std::vector<std::vector<std::complex<double>>>& M,
        std::vector<size_t>& pivot);
    static bool CmpGaussianSolution(
        int n, std::vector<std::vector<std::complex<double>>>& M,
        std::vector<std::complex<double>>& b,
        std::vector<std::complex<double>>& x,
        std::vector<size_t>& pivot);
    static bool CmpSubstitution(
        int m, int n, std::vector<std::vector<std::complex<double>>>& M,
        std::vector<std::complex<double>>& b,
        std::vector<std::complex<double>>& x,
        std::vector<size_t>& pivot);
    static bool CmpInverse(
        int n, std::vector<std::vector<std::complex<double>>>& M,
        std::vector<std::vector<std::complex<double>>>& Mi);
    static void CmpCharPolyAndAdjoint(
        int n,
        std::vector<std::vector<std::complex<double>>>& C,
        std::vector<std::vector<std::complex<double>>>& I,
        std::vector<std::vector<std::complex<double>>>& M,
        std::vector<std::vector<std::complex<double>>>& adjoint,
        std::vector<std::complex<double>>& a);
    static void CmpMatrixKernel(
        int m, int n,
        std::vector<std::vector<std::complex<double>>>& M,
        std::vector<std::vector<std::complex<double>>>& X,
        size_t& r);
    static void CmpMatrixImage(
        int m, int n,
        std::vector<std::vector<std::complex<double>>>& M,
        std::vector<std::vector<std::complex<double>>>& N,
        std::vector<std::vector<std::complex<double>>>& X,
        int rank);
};

#include <vector>

class DblLinearAlgebra
{
public:
    static void DblPrintMatrix(
        int m, int n, std::vector<std::vector<double>>& A);
    static void DblAddition(
        size_t m, size_t n,
        std::vector<std::vector<double>>& A,
        std::vector<std::vector<double>>& B,
        std::vector<std::vector<double>>& C);
    static void DblSubtraction(
        size_t m, size_t n,
        std::vector<std::vector<double>>& A,
        std::vector<std::vector<double>>& B,
        std::vector<std::vector<double>>& C);
    static void DblMultiply(
        size_t m, size_t n, size_t p,
        std::vector<std::vector<double>>& A,
        std::vector<std::vector<double>>& B,
        std::vector<std::vector<double>>& C);
    static void DblAnticommutator(
        size_t n,
        std::vector<std::vector<double>>& A,
        std::vector<std::vector<double>>& B,
        std::vector<std::vector<double>>& C);
    static void DblCommutator(
        size_t n,
        std::vector<std::vector<double>>& A,
        std::vector<std::vector<double>>& B,
        std::vector<std::vector<double>>& C);
    static double DblDeterminant(
        bool& failure, int n,
        std::vector<std::vector<double>>& A);
    static bool DblGaussianElimination(
        int m, int n, std::vector<std::vector<double>>& A,
        std::vector<double>& b, std::vector<double>& x,
        std::vector<size_t>& pivot);
    static bool DblGaussianFactor(
        int n, std::vector<std::vector<double>>& M,
        std::vector<size_t>& pivot);
    static bool DblGaussianSolution(
        int n, std::vector<std::vector<double>>& M,
        std::vector<double>& b, std::vector<double>& x,
        std::vector<size_t>& pivot);
    static bool DblSubstitution(
        int n, std::vector<std::vector<double>>& M,
        std::vector<double>& b, std::vector<double>& x,
        std::vector<size_t>& pivot);
    static bool DblInverse(
        int n, std::vector<std::vector<double>>& M,
        std::vector<std::vector<double>>& A);
    static void DblCharPolyAndAdjoint(
        int n,
        std::vector<std::vector<double>>& C,
        std::vector<std::vector<double>>& I,
        std::vector<std::vector<double>>& M,
        std::vector<std::vector<double>>& adjoint,
        std::vector<double>& a);
    static void DblMatrixKernel(
        int m, int n,
        std::vector<std::vector<double>>& M,
        std::vector<std::vector<double>>& X,
        size_t& r);
    static void DblMatrixImage(
        int m, int n,
        std::vector<std::vector<double>>& M,
        std::vector<std::vector<double>>& N,
        std::vector<std::vector<double>>& X,
        int rank);
};

// Exercises from "Modern Quantum Chemistry An Introduction to Advanced
// Electronic Structure Theory" by Attila Szabo and Neil S. Ostlund
// https://chemistlibrary.wordpress.com/wp-content/uploads/2015/02/modern-quantum-chemistry.pdf
// Program (c) Tuesday, August 12, 2025 by James Pate Williams, Jr.

#include <complex>
#include <iomanip>
#include <iostream>
#include <vector>
#include "DblLinearAlgebra.h"
#include "CmpLinearAlgebra.h"

int main()
{
	double AArcb[3][3] = { { 2, 3, -1 }, { 4, 4, -3 }, { -2, 3, -1 } };
	double AArso[3][3] = { { 1, 1, 0 }, { 1, 2, 2 }, { 0, 2, -1 } };
	double BBrso[3][3] = { { 1, -1, 1 }, { -1, 0, 0 }, { 1, 0, 1} };
	double BBr[3][3] = { { 1, -1, 1 }, { -1 , 0, 0 }, { 1, 0, 1 } };
	double AAcr[3][3] = { { 1, 2, 3 }, { 4, 5, 6 }, { 7, 8, 9 } };
	double AAci[3][3] = { { 1, 1, 2 }, { 3, 0, 1 }, { 0, 2, 4 } };
	double BBcr[3][3] = { { 1, 0, 1 }, { 1 , 1, 0 }, { 0, 1, 1 } };
	double BBci[3][3] = { { 1, 2, 3 }, { 4, 5, 6 }, { 7, 8, 9 } };
	int m = 3, n = 3, p = 3;

	std::vector<double> br(3);
	std::vector<size_t> pivot(3);
	std::vector<std::vector<double>> Arcb(3, std::vector<double>(3));
	std::vector<std::vector<double>> Arso(3, std::vector<double>(3));
	std::vector<std::vector<double>> Brso(3, std::vector<double>(3));
	std::vector<std::vector<double>> Br(3, std::vector<double>(3));
	std::vector<std::vector<double>> Cr(3, std::vector<double>(3));
	std::vector<std::vector<double>> Ai(3, std::vector<double>(3));
	std::vector<std::vector<double>> Ari(3, std::vector<double>(3));
	std::vector<std::vector<std::complex<double>>> Ac(3,
		std::vector<std::complex<double>>(3));
	std::vector<std::vector<std::complex<double>>> Bc(3,
		std::vector<std::complex<double>>(3));
	std::vector<std::vector<std::complex<double>>> Cc(3,
		std::vector<std::complex<double>>(3));
	std::vector<std::vector<std::complex<double>>> Dc(3,
		std::vector<std::complex<double>>(3));
	std::vector<std::vector<std::complex<double>>> Ec(3,
		std::vector<std::complex<double>>(3));
	std::vector<std::vector<std::complex<double>>> Fc(3,
		std::vector<std::complex<double>>(3));
	std::vector<std::vector<std::complex<double>>> Gc(3,
		std::vector<std::complex<double>>(3));
	
	for (int i = 0; i < m; i++)
	{
		for (int j = 0; j < p; j++)
		{
			Arcb[i][j] = AArcb[i][j];
			Arso[i][j] = AArso[i][j];
			Brso[i][j] = BBrso[i][j];
			Ac[i][j]._Val[0] = AAcr[i][j];
			Ac[i][j]._Val[1] = AAci[i][j];
		}
	}

	for (int i = 0; i < p; i++)
	{
		for (int j = 0; j < n; j++)
		{
			Br[i][j] = BBr[i][j];
			Bc[i][j]._Val[0] = BBcr[i][j];
			Bc[i][j]._Val[1] = BBci[i][j];
		}
	}
	
	DblLinearAlgebra::DblMultiply(3, 3, 3, Arcb, Br, Cr);
	std::cout << "Ar * Br = Cr Conte & de Boor" << std::endl;
	DblLinearAlgebra::DblPrintMatrix(3, 3, Cr);
	std::cout << std::endl;
	CmpLinearAlgebra::CmpMultiply(3, 3, 3, Ac, Bc, Cc);
	std::cout << "Ac * Bc = Cc" << std::endl;
	CmpLinearAlgebra::CmpPrintMatrix(3, 3, Cc);
	std::cout << std::endl;
	// Exercise 1.2
	std::cout << "Exercise 1.2 page 5 Commutator" << std::endl;
	DblLinearAlgebra::DblCommutator(3, Arso, Brso, Cr);
	DblLinearAlgebra::DblPrintMatrix(3, 3, Cr);
	std::cout << std::endl;
	std::cout << "Exercise 1.2 page 5 Anticommutator" << std::endl;
	DblLinearAlgebra::DblAnticommutator(3, Arso, Brso, Cr);
	DblLinearAlgebra::DblPrintMatrix(3, 3, Cr);
	std::cout << std::endl;
	CmpLinearAlgebra::CmpAdjoint(3, 3, Cc, Dc);
	std::cout << "Exercise 1.3 page 6 Cc adjoint" << std::endl;
	CmpLinearAlgebra::CmpPrintMatrix(3, 3, Dc);
	std::cout << std::endl;
	CmpLinearAlgebra::CmpAdjoint(3, 3, Ac, Ec);
	CmpLinearAlgebra::CmpAdjoint(3, 3, Bc, Fc);
	CmpLinearAlgebra::CmpMultiply(3, 3, 3, Fc, Ec, Gc);
	std::cout << "Exercise 1.3 page 6 Bc adjoint * Ac adjoint" << std::endl;
	CmpLinearAlgebra::CmpPrintMatrix(3, 3, Gc);
	std::cout << std::endl;
	std::cout << "Ar matrix" << std::endl;
	DblLinearAlgebra::DblPrintMatrix(3, 3, Arcb);
	bool inv = DblLinearAlgebra::DblInverse(n, Arcb, Ai);
	std::cout << std::endl;
	std::cout << "Ar Conte & de Boor inverse = " << inv << std::endl;
	DblLinearAlgebra::DblPrintMatrix(3, 3, Ai);
	std::cout << std::endl;
	std::cout << "Ar * Ar inverse" << std::endl;
	DblLinearAlgebra::DblMultiply(3, 3, 3, Arcb, Ai, Ari);
	DblLinearAlgebra::DblPrintMatrix(3, 3, Ari);
	std::cout << std::endl;
	std::cout << "Ac" << std::endl;
	CmpLinearAlgebra::CmpPrintMatrix(3, 3, Ac);
	std::cout << std::endl;
	inv = CmpLinearAlgebra::CmpInverse(3, Ac, Bc);
	std::cout << "Ac inverse = " << inv << std::endl;
	CmpLinearAlgebra::CmpPrintMatrix(3, 3, Bc);
	CmpLinearAlgebra::CmpMultiply(3, 3, 3, Ac, Bc, Cc);
	std::cout << std::endl;
	std::cout << "Ac * AC inverse" << std::endl;
	CmpLinearAlgebra::CmpPrintMatrix(3, 3, Cc);
}

#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; }

Category: Elementary Physics

Blog Entry © Saturday, August 16, 2025, by James Pate Williams, Jr. Some More Elementary Quantum Chemistry

Blog Entry © Friday, August 15, 2025, by James Pate Williams, Jr. Some Elementary Quantum Chemistry

Blog Entry (c) Wednesday, August 13, 2025, by James Pate Williams, Jr. Exercises from an Online Textbook

Blog Entry © Monday, August 11, 2025, by James Pate Williams, Jr. Van der Waals Interaction Between Two Hydrogen Atoms

Blog Entry © Sunday, August 10, 2025, First-Order Perturbation Treatment of the Helium Atom by James Pate Williams, Jr.

Blog Entry © Thursday, August 7, 2025, z-Component of Cylindrical Velocity Based on an Equation in Classical Electrodynamics Second Edition © 1975 by John David Jackson

Blog Entry © Saturday, August 2, 2025, by James Pate Williams, Jr. Orthonormal Transverse Magnetic and Transverse Electric Fields Win32 Desktop C/C++ Application in the Release Configuration

Blog Entry © Tuesday, July 29, 2025, Double and Triple Monte Carlo Integration by James Pate Williams, Jr.

Blog Entry Thursday, July 24, 2025, © James Pate Williams, Jr. Comparison of Data Between Two Sources

Newtonian Gravity and Cosmology © Friday, July 18, 2025, by James Pate Williams, Jr. in Conjunction with the Microsoft Copilot AI