Category: Numerical Analysis

Blog Entry © Thursday, May 22, 2025, Baseball Ballistics with Simple Drag Computation by James Pate Williams, Jr.

Blog Entry Thursday May 22 2025Download
// https://phys.libretexts.org/Bookshelves/University_Physics/Physics_(Boundless)/3%3A_Two-Dimensional_Kinematics/3.3%3A_Projectile_Motion

#pragma once

class ClassicalBallistics
{
public:
	ClassicalBallistics(void);
	~ClassicalBallistics(void);
	// analytical formula for maximum height
	double H(double V0, double theta0, double g);
	// analytical formula for flight time
	double T(double H, double g);
	// analytical formula for velocity at apex
	double Va(double V0, double theta0);
	// analytical formula for the range
	double L(double Va, double T);
	// analytical formula for time to apex
	double ta(double T);
	// analytical formula for range to apex
	double xa(double L);
	double xa(double H, double L, double theta0);
	double Theta1(double theta0);
	double Theta1(double H, double L, double xa);
	double ClassicalBallistics::V1(double V0);
};


#include "StdAfx.h"
#include "ClassicalBallistics.h"
#include <math.h>

ClassicalBallistics::ClassicalBallistics(void)
{
}
ClassicalBallistics::~ClassicalBallistics(void)
{
}
// analytical formula for maximum height
double ClassicalBallistics::H(double V0, double theta0, double g)
{
	double V02 = V0 * V0, sinTheta0 = sin(theta0);
	double sinTheta02 = sinTheta0 * sinTheta0;

	return V02 * sinTheta02 / (g * 2.0);
}
// analytical formula for flight time
double ClassicalBallistics::T(double H, double g)
{
	return 2.0 * sqrt(2.0 * H / g);
}
// analytical formula for velocity at apex
double ClassicalBallistics::Va(double V0, double theta0)
{
	double cosTheta0 = cos(theta0);
	double V02 = V0 * V0;

	return V0 * cosTheta0;
}
// analytical formula for the range
double ClassicalBallistics::L(double Va, double T)
{
	return Va * T;
}
// analytical formula for time to apex
double ClassicalBallistics::ta(double T)
{
	return 0.5 * T;
}
// analytical formula for range to apex
double ClassicalBallistics::xa(double L)
{
	return 0.5 * L;
}
double ClassicalBallistics::xa(double H, double L, double theta0)
{
	return sqrt(L * H / tan(theta0));
}
double ClassicalBallistics::Theta1(double theta0)
{
	return -theta0;
}
double ClassicalBallistics::Theta1(double H, double L, double xa)
{
	return atan(L * H / pow(L - xa, 2.0));
}
double ClassicalBallistics::V1(double V0)
{ 
	return V0;
}


// http://www.lajpe.org/sep13/04-LAJPE-782_Chudinov.pdf

#pragma once

const double g = 9.81;
// accleration due to gravity in Metric (SI) Units
const double k = 0.000625;
// k is typical for a baseball

class DragBallistics
{

public:

	DragBallistics(void);
	~DragBallistics(void);
    double f(double theta);
    // analytical formula for maximum height
    double H(double V0, double theta0, double g, double k);
    // analytical formula for flight time
    double T(double H, double g);
    // analytical formula for velocity at apex
    double Va(double V0, double theta0, double k);
    // analytical formula for the range
	double L(double Va, double T);
    double ta(double H, double T, double Va, double k);
    double xa(double H, double L, double theta0);
    double Theta1(double H, double L, double xa);
    double V1(double V0, double theta0, double theta1, double k);
};


#include "StdAfx.h"
#include "DragBallistics.h"
#include <math.h>

DragBallistics::DragBallistics(void)
{
}
DragBallistics::~DragBallistics(void)
{
}
double DragBallistics::f(double theta)
{
	double pi = 4.0 * atan(1.0);
	double sinTheta = sin(theta), cosTheta = cos(theta),
		cosTheta2 = cosTheta * cosTheta;
    return (sinTheta / cosTheta2) + log(tan(0.5 * theta + pi / 4.0));
}
// analytical formula for maximum height
double DragBallistics::H(double V0, double theta0, double g, double k)
{
	double V02 = V0 * V0, sinTheta0 = sin(theta0),
		sinTheta02 = sinTheta0 * sinTheta0;
    return V02 * sinTheta02 / (g * (2.0 + k * V02 * sinTheta0));
}
// analytical formula for flight time
double DragBallistics::T(double H, double g)
{
	return 2.0 * sqrt(2.0 * H / g);
}
// analytical formula for velocity at apex
double DragBallistics::Va(double V0, double theta0, double k)
{
	double pi = 4.0 * atan(1.0);
	double cosTheta0 = cos(theta0), cosTheta02 = cosTheta0 * cosTheta0;
    double V02 = V0 * V0, sinTheta0 = sin(theta0);
    return V0 * cosTheta0 / sqrt(1.0 + k * V02 * (sinTheta0 + 
		cosTheta02 * log(tan(0.5 * theta0 + pi / 4.0))));
}
// analytical formula for the range
double DragBallistics::L(double Va, double T)
{
	return Va * T;
}
double DragBallistics::ta(double H, double T, double Va, double k)
{
	return 0.5 * (T - k * H * Va);
}
double DragBallistics::xa(double H, double L, double theta0)
{
	return sqrt(L * H / tan(theta0));
}
double DragBallistics::Theta1(double H, double L, double xa)
{
	return -atan(L * H / pow(L - xa, 2.0));
}
double DragBallistics::V1(double V0, double theta0, double theta1, double k)
{
	double V02 = V0 * V0;
    double sinTheta0 = sin(theta0), cosTheta0 = cos(theta0),
		cosTheta02 = cosTheta0 * cosTheta0;
    return V0 * cosTheta0 / (cos(theta1) * sqrt(1.0 + k * 
		V02 * cosTheta02 * (f(theta0) - f(theta1))));
}


// BaseballBallisticsWin32Console.cpp
// Translated from August 2017 C# application
// May 21, 2025 (c) James Pate Williams, Jr.

#include "stdafx.h"
#include "ClassicalBallistics.h"
#include "DragBallistics.h"
#include <iomanip>
#include <iostream>
#include <math.h>

void PrintResults(char title[],
				  double H, double T, double Va,
				  double L, double Ta, double xa,
				  double Theta1, double V1)
{
	std::cout << title << std::endl;
	std::cout << std::fixed << std::setprecision(2);
	std::cout << "H      = " << H << std::endl;
	std::cout << "T      = " << T << std::endl;
	std::cout << "Va     = " << Va << std::endl;
	std::cout << "L      = " << L << std::endl;
	std::cout << "Ta     = " << Ta << std::endl;
	std::cout << "xa     = " << xa << std::endl;
	std::cout << "theta1 = " << Theta1 << std::endl;
	std::cout << "V1     = " << V1 << std::endl;
}

int _tmain(int argc, _TCHAR* argv[])
{
	while (true)
	{
		char line[128] = { };
		std::cout << "V0 (m / s) or 0 to quit: ";
		std::cin.getline(line, 128);
		double V0 = atof(line);
		if (V0 == 0)
		{
			break;
		}
		std::cout << "Enter angle in degrees: ";
		std::cin.getline(line, 128);
		double theta0 = atof(line);
		double pi = 4.0 * atan(1.0);
		theta0 *= pi / 180.0;
		ClassicalBallistics cBall;
		DragBallistics dBall;
		double cH = cBall.H(V0, theta0, g);
        double dH = dBall.H(V0, theta0, g, k);
        double cT = cBall.T(cH, g);
        double dT = dBall.T(dH, g);
        double cVa = cBall.Va(V0, theta0);
        double dVa = dBall.Va(V0, theta0, k);
        double cL = cBall.L(cVa, cT);
        double dL = dBall.L(dVa, dT);
        double cta = cBall.ta(cT);
        double dta = dBall.ta(dH, dT, dVa, k);
        double cxa = cBall.xa(cL);
        double dxa = dBall.xa(dH, dL, theta0);
        double cTheta1 = 180.0 * cBall.Theta1(theta0) / pi;
        double dTheta1 = 180.0 * dBall.Theta1(dH, dL, dxa) / pi;
        double cV1 = cBall.V1(V0);
        double dV1 = dBall.V1(V0, theta0, pi * dTheta1 / 180.0, k);
		PrintResults("Classical Ballistics",
			cH, cT, cVa, cL, cta, cxa, cTheta1, cV1);
		PrintResults("Drag Ballistics",
			dH, dT, dVa, dL, dta, dxa, dTheta1, dV1);
	}
	return 0;
}

Blog Entry © Wednesday, May 21, 2025, Backpropagation Artificial Neural Network Experiments by James Pate Williams, Jr.

Backpropagation Artificial Neural Network ExperimentsDownload

Learning a Simple Function Using an Artificial Backpropagation Neural Network © Wednesday, May 14, 2025, by James Pate Williams, Jr., BA, BS, Master of Software Engineering, and PhD Computer Science

Learning a Simple Function Using an Artificial Backpropagation Neural NetworkDownload

Four Methods of Numerical Double Integration: Sequential Simpson’s Rule, Multitasking Simpson’s Rule, Sequential Monte Carlo and Multitasking Monte Carlo Methods © Wednesday April 16 – 18, 2025, by James Pate Williams, Jr.

Four Methods of Numerical Double IntegrationDownload

Blog Entry (c) Friday, April 11, 2025, by James Pate Williams, Jr. Multiplication and Division of Finite Power Series

x = 0.25
N = 5
Series Cosine(x)  = 0.968912421711
C++ cos(x)        = 0.968912421711
Series sine(x)    = 0.247403959255
C++ sin(x)        = 0.247403959255
Series Tangent(x) = 0.255341921221
C++ tan(x)        = 0.255341921221
Series sin(2x)    = 0.479425538604
C++ sin(2x)       = 0.479425538604
C++ 2sin(x)cos(x) = 0.479425538604
End app ? y = yes = n
x = 0.5
N = 5
Series Cosine(x)  = 0.877582561890
C++ cos(x)        = 0.877582561890
Series sine(x)    = 0.479425538604
C++ sin(x)        = 0.479425538604
Series Tangent(x) = 0.546302489844
C++ tan(x)        = 0.546302489844
Series sin(2x)    = 0.841470984807
C++ sin(2x)       = 0.841470984808
C++ 2sin(x)cos(x) = 0.841470984808
End app ? y = yes = y
x = 0.75
N = 5
Series Cosine(x)  = 0.731688868808
C++ cos(x)        = 0.731688868874
Series sine(x)    = 0.681638760020
C++ sin(x)        = 0.681638760023
Series Tangent(x) = 0.931596460023
C++ tan(x)        = 0.931596459944
Series sin(2x)    = 0.997494986509
C++ sin(2x)       = 0.997494986604
C++ 2sin(x)cos(x) = 0.997494986604
End app ? y = yes = n
 
x = 1.00
N = 5
Series Cosine(x)  = 0.540302303792
C++ cos(x)        = 0.540302305868
Series sine(x)    = 0.841470984648
C++ sin(x)        = 0.841470984808
Series Tangent(x) = 1.557407730344
C++ tan(x)        = 1.557407724655
Series sin(2x)    = 0.909297423159
C++ sin(2x)       = 0.909297426826
C++ 2sin(x)cos(x) = 0.909297426826
End app ? y = yes = n
x = 1.25
N = 5
Series Cosine(x)  = 0.315322332275
C++ cos(x)        = 0.315322362395
Series sine(x)    = 0.948984616456
C++ sin(x)        = 0.948984619356
Series Tangent(x) = 3.009569952151
C++ tan(x)        = 3.009569673863
Series sin(2x)    = 0.598472085108
C++ sin(2x)       = 0.598472144104
C++ 2sin(x)cos(x) = 0.598472144104
End app ? y = yes = n
x = 1.50
N = 5
Series Cosine(x)  = 0.070736934117
C++ cos(x)        = 0.070737201668
Series sine(x)    = 0.997494955682
C++ sin(x)        = 0.997494986604
Series Tangent(x) = 14.101472846329
C++ tan(x)        = 14.101419947172
Series sin(2x)    = 0.141119469924
C++ sin(2x)       = 0.141120008060
C++ 2sin(x)cos(x) = 0.141120008060
End app ? y = yes = y


// DivMulPowerSeries.cpp (c) Tuesday, April 8, 2025
// by James Pate Williams, Jr.
// https://math.libretexts.org/Bookshelves/Calculus/Calculus_(OpenStax)/10%3A_Power_Series/10.02%3A_Properties_of_Power_Series
// https://en.wikipedia.org/wiki/Formal_power_series

#include <iostream>
#include <iomanip>
#include <string>
#include <vector>

std::vector<double> Multiplication(
	std::vector<double> c,
	std::vector<double> d,
	int N) {
	std::vector<double> e(N + 1);

	for (int n = 0; n <= N; n++) {
		double sum = 0.0;

		for (int k = 0; k <= n; k++) {
			sum += c[k] * d[n - k];
		}

		e[n] = sum;
	}

	return e;
}

std::vector<double> Division(
	std::vector<double> a,
	std::vector<double> b,
	int N) {
	double a0 = a[0];
	std::vector<double> c(N + 1);

	for (int n = 0; n <= N; n++) {
		double sum = 0.0;

		for (int k = 1; k <= n; k++) {
			sum += a[k] * c[n - k];
		}

		c[n] = (b[n] - sum) / a0;
	}

	return c;
}

double Factorial(int n) {
	double nf = 1.0;

	for (int i = 2; i <= n; i++) {
		nf *= i;
	}

	return nf;
}

std::vector<double> Cosine(double x, int N, double& fx) {
	std::vector<double> series(N + 1);

	fx = 0.0;

	for (int n = 0; n <= N; n++) {
		int argument = 2 * n;
		double coeff = pow(-1, n) / Factorial(argument);
		series[n] = coeff;
		fx += coeff * pow(x, argument);
	}

	return series;
}

std::vector<double> Sine(double x, int N, double& fx) {
	std::vector<double> series(N + 1);
	
	fx = 0.0;

	for (int n = 0; n <= N; n++) {
		int argument = 2 * n + 1;
		double coeff = pow(-1, n) / Factorial(argument);
		series[n] = coeff;
		fx += coeff * pow(x, argument);
	}

	return series;
}

std::vector<double> Tangent(double x, int N, double& fx) {
	double fc = 0.0, fs = 0.0;
	std::vector<double> seriesC = Cosine(x, N, fc);
	std::vector<double> seriesS = Sine(x, N, fs);
	std::vector<double> seriesT = Division(seriesS, seriesC, N);
	fx = fs / fc;

	return seriesT;
}

std::vector<double> Sine2x(double x, int N, double& fx) {
	double fc = 0.0, fs = 0.0;
	std::vector<double> seriesC = Cosine(x, N, fc);
	std::vector<double> seriesS = Sine(x, N, fs);
	std::vector<double> series2 = Multiplication(seriesS, seriesC, N);
	fx = 2.0 * fs * fc;

	return series2;
}

int main()
{
	while (true) {
		char line[128] = { };
		std::cout << "x = ";
		std::cin.getline(line, 127);
		std::string str1(line);
		double x = std::stod(str1);
		std::cout << "N = ";
		std::cin.getline(line, 127);
		std::string str2(line);
		int N = std::stoi(str2);
		double cx = 0.0, sx = 0.0, tx = 0.0, xx = 0.0;
		std::vector<double> cSeries = Cosine(x, N, cx);
		std::vector<double> sSeries = Sine(x, N, sx);
		std::vector<double> tSeries = Tangent(x, N, tx);
		std::vector<double> xSeries = Sine2x(x, N, xx);
		std::cout << std::fixed << std::setprecision(12);
		std::cout << "Series Cosine(x)  = " << cx << std::endl;
		std::cout << "C++ cos(x)        = " << cos(x) << std::endl;
		std::cout << "Series sine(x)    = " << sx << std::endl;
		std::cout << "C++ sin(x)        = " << sin(x) << std::endl;
		std::cout << "Series Tangent(x) = " << tx << std::endl;
		std::cout << "C++ tan(x)        = " << tan(x) << std::endl;
		std::cout << "Series sin(2x)    = " << xx << std::endl;
		std::cout << "C++ sin(2x)       = " << sin(x + x) << std::endl;
		std::cout << "C++ 2sin(x)cos(x) = " << 2.0 * sin(x) * cos(x);
		std::cout << std::endl;
		std::cout << "End app ? y = yes = ";
		std::cin.getline(line, 127);

		if (line[0] == 'Y' || line[0] == 'y') {
			break;
		}
	}

	return 0;
}

Blog Entry, Tuesday, April 8, 2025, (c) James Pate Williams, Jr. PROBLEMS from “Mathematical Methods in the Physical Sciences Second Edition” (c) 1983 by Mary L. Boas CHAPTER 1 SECTION 13

Boas Problems Section 13Download

Approximation of the Ground-State Total Energy of a Beryllium Atom © Sunday, March 30 to Tuesday April 1, 2025, by James Pate Williams, Jr., BA, BS, Master of Software Engineering, PhD Computer Science

Approximation of the Ground-State Total Energy of a Beryllium AtomDownload

Blog Entry © Sunday, March 29, 2025, by James Pate Williams, Jr., BA, BS, Master of Software Engineering, PhD Slater Determinant Coefficients for Z = 2 to 4

Enter the atomic number Z (2 to 6 or 0 to quit): 2
2       1       1       +       a(1)b(2)
1       0       0       -       a(2)b(1)
# Even Permutations = 1
Enter the atomic number Z (2 to 6 or 0 to quit): 3
6       3       1       +       a(1)b(2)c(3)
5       2       0       -       a(1)b(3)c(2)
4       2       0       -       a(2)b(1)c(3)
3       1       1       +       a(2)b(3)c(1)
2       1       1       +       a(3)b(1)c(2)
1       0       0       -       a(3)b(2)c(1)
# Even Permutations = 3
Enter the atomic number Z (2 to 6 or 0 to quit): 4
24      12      0       +       a(1)b(2)c(3)d(4)
23      11      1       -       a(1)b(2)c(4)d(3)
22      11      1       -       a(1)b(3)c(2)d(4)
21      10      0       +       a(1)b(3)c(4)d(2)
20      10      0       +       a(1)b(4)c(2)d(3)
19      9       1       -       a(1)b(4)c(3)d(2)
18      9       1       -       a(2)b(1)c(3)d(4)
17      8       0       +       a(2)b(1)c(4)d(3)
16      8       0       +       a(2)b(3)c(1)d(4)
15      7       1       -       a(2)b(3)c(4)d(1)
14      7       1       -       a(2)b(4)c(1)d(3)
13      6       0       +       a(2)b(4)c(3)d(1)
12      6       0       +       a(3)b(1)c(2)d(4)
11      5       1       -       a(3)b(1)c(4)d(2)
10      5       1       -       a(3)b(2)c(1)d(4)
9       4       0       +       a(3)b(2)c(4)d(1)
8       4       0       +       a(3)b(4)c(1)d(2)
7       3       1       -       a(3)b(4)c(2)d(1)
6       3       1       -       a(4)b(1)c(2)d(3)
5       2       0       +       a(4)b(1)c(3)d(2)
4       2       0       +       a(4)b(2)c(1)d(3)
3       1       1       -       a(4)b(2)c(3)d(1)
2       1       1       -       a(4)b(3)c(1)d(2)
1       0       0       +       a(4)b(3)c(2)d(1)
# Even Permutations = 12
Enter the atomic number Z (2 to 6 or 0 to quit):

// AOPermutations.cpp : This file contains the 'main' function.
// Program execution begins and ends there.
// Copyright (c) Saturday, March 29, 2025
// by James Pate Williams, Jr., BA, BS, MSwE, PhD
// Signs of the atomic orbitals in a Slater Determinant

#include <algorithm>
#include <iostream>
#include <string>
#include <vector>

int main()
{
    char alpha[] = { 'a', 'b', 'c', 'd', 'e', 'f' }, line[128] = {};
    int factorial[7] = { 1, 1, 2, 6, 24, 120, 720 };

    while (true)
    {
        int col = 0, counter = 0, row = 0, sign = 1, t = 0, Z = 0, zfact = 0;
        int numberEven = 0;
        std::cout << "Enter the atomic number Z (2 to 6 or 0 to quit): ";
        std::cin.getline(line, 127);
        std::string str(line);
        Z = std::stoi(str);

        if (Z == 0)
        {
            break;
        }

        if (Z < 2 || Z > 6)
        {
            std::cout << "Illegal Z, please try again" << std::endl;
            continue;
        }

        zfact = factorial[Z];

        std::vector<char> orb(Z);
        std::vector<int> tmp(Z), vec(Z);

        for (int i = 0; i < Z; i++)
        {
            orb[i] = alpha[i];
            vec[i] = i + 1;
        }

        do
        {
            for (int i = 0; i < (int)vec.size(); i++)
            {
                tmp[i] = vec[i];
            }

            t = 0;

            do
            {
                t++;
            } while (std::next_permutation(tmp.begin(), tmp.end()));

            std::cout << t << '\t' << t / 2 << '\t';
            std::cout << (t / 2 & 1) << '\t';

            if (Z == 2 || Z == 3)
            {
                if ((t / 2 & 1) == 0)
                {
                    std::cout << "-\t";
                }

                else
                {
                    std::cout << "+\t";
                    numberEven++;
                }
            }

            else
            {
                if ((t / 2 & 1) == 1)
                {
                    std::cout << "-\t";
                }

                else
                {
                    std::cout << "+\t";
                    numberEven++;
                }
            }

            for (int i = 0; i < Z; i++)
            {
                std::cout << orb[i] << '(' << vec[i] << ')';
            }

            row++;
            std::cout << std::endl;

            if (zfact != 2 && row == zfact)
            {
                std::cout << std::endl;
                break;
            }

            row %= Z;
        } while (std::next_permutation(vec.begin(), vec.end()));

        std::cout << "# Even Permutations = ";
        std::cout << numberEven << std::endl;
    }

    return 0;
}