“Babe Ruth is generally considered the owner of the record for the longest home run in MLB history with a 575-foot bomb launched at Navin Field in Detroit in 1921.” – https://www.msn.com/en-us/sports/mlb/what-is-the-longest-home-run-in-mlb-history/ar-AA1dGwlZ
As you can see, I estimated the pitch velocity at 90 miles per hour and Babe Ruth’s (Sultan of Swing) at 90 miles per hour also. My analytic calculations yield a range of the baseball’s trajectory as about 576 feet.
// "Numerical Computation 2: Methods, Software,
// and Analysis" by Christoph W. Ueberhuber
// Chapter 17 Random Numbers
// "The Art of Computer Programming Volume 2"
// "Seminumerical Algorithms Second Edition"
// "Chapter 3 RANDOM NUMBERS" by Donald E. Knuth
// https://en.wikipedia.org/wiki/Mersenne_Twister
using System.Collections.Generic;
namespace SimplePRNGs
{
class PRNGs
{
private long AMxn1, AMxn;
private long AMyn1, AMyn;
private long AMk;
private long AMm = 34359738368;
private long LCG0z0, LCG0z1;
private long LCG1z0, LCG1z1;
private long LCG2z0, LCG2z1, LCG2z2;
private long MFG0z0, MFG0z1;
private readonly long LCG0m = 4294967296;
private readonly long LCG2m = 2147483647;
private readonly List<long> AMV = new();
private long MTindex;
private long[] MT;
private LaggedFibRng fibRng;
public void SetSeedLCG0(long z0)
{
LCG0z0 = z0;
}
public void SetSeedLCG1(long z0)
{
LCG1z0 = z0;
}
public void SetSeedLCG2(long z0, long z1)
{
LCG2z0 = z0;
LCG2z1 = z1;
}
public long LCG0()
{
LCG0z1 = (69069 * LCG0z0) % LCG0m;
LCG0z0 = LCG0z1;
return LCG0z1;
}
public long LCG1()
{
LCG1z1 = (69069 * LCG1z0 + 1) % LCG0m;
LCG1z0 = LCG1z1;
return LCG1z1;
}
public long LCG2()
{
LCG2z2 = (1999 * LCG2z1 + 4444 * LCG2z0) % LCG2m;
LCG2z0 = LCG2z1;
LCG2z1 = LCG2z2;
return LCG2z2;
}
public void SetSeedMFG0(long z0, long z1)
{
MFG0z0 = z0;
MFG0z1 = z1;
}
public long MFG0()
{
long MFG0z2 = (MFG0z1 + MFG0z0) % LCG0m;
MFG0z0 = MFG0z1;
MFG0z1 = MFG0z2;
return MFG0z2;
}
public void ComputeNextXY()
{
AMxn1 = (3141592653 * AMxn + 2718281829) % AMm;
if (AMxn1 < 0)
AMxn1 += AMm;
AMyn1 = (2718281829 * AMyn + 3141592653) % AMm;
if (AMyn1 < 0)
AMyn1 += AMm;
}
public void AMSeed(long k, long x0, long y0)
{
long AMTxn1, AMTxn = x0;
AMxn = x0;
AMyn = y0;
AMk = k;
for (int i = 0; i < k; i++)
{
AMTxn1 = (3141592653 * AMTxn + 2718281829) % AMm;
if (AMTxn1 < 0)
AMTxn1 += AMm;
AMTxn = AMTxn1;
AMV.Add(AMTxn1);
}
}
public long AlgorithmM()
{
ComputeNextXY();
AMxn = AMxn1;
AMyn = AMyn1;
long j = (AMk * AMyn1) / AMm;
long r = AMV[(int)j];
AMV[(int)j] = AMxn1;
if (r < 0)
r += AMm;
return r;
}
public void MTInitialization(long seed)
{
long f = 6364136223846793005;
long n = 312, w = 64;
MTindex = n;
MT = new long[n];
MT[0] = seed;
for (int i = 1; i < n; i++)
MT[i] = f * (MT[i - 1] ^ (MT[i - 1] >> (int)(w - 2))) + i;
}
public long MTExtractNumber()
{
unchecked
{
long n = 312;
long c = (long)0xFFF7EEE000000000;
long b = 0x71D67FFFEDA60000;
long d = 0x5555555555555555;
long u = 29, s = 17, t = 27, l = 43;
if (MTindex == n)
MTTwist();
long y = MT[MTindex];
y ^= ((y >> (int)u) & d);
y ^= ((y << (int)s) & b);
y ^= ((y << (int)t) & c);
y ^= (y >> (int)l);
MTindex++;
return y;
}
}
public void MTTwist()
{
unchecked
{
long n = 312, m = 156, r = 31;
long a = (long)0xB5026F5AA96619E9;
MTindex = n + 1;
long lower_mask = (1 << (int)r) - 1;
long upper_mask = ~lower_mask;
for (int i = 0; i < n; i++)
{
long x = (MT[i] & upper_mask) |
(MT[(i + 1) % n] & lower_mask);
long xA = x >> 1;
if (x % 2 != 0)
xA ^= a;
MT[i] = MT[(i + m) % n] ^ xA;
}
}
MTindex = 0;
}
public void LaggedFibRngSeed(int seed)
{
fibRng = new LaggedFibRng(seed);
}
public long LaggedFibonacci(long modulus)
{
long lo = fibRng.Next();
long hi = fibRng.Next();
long rs = ((hi << 31) | lo) % modulus;
return rs;
}
}
}
//https://learn.microsoft.com/en-us/archive/msdn-magazine/2016/august/test-run-lightweight-random-number-generation
// modified by current author James Pate Williams, Jr. on August 30, 2023
using System.Collections.Generic;
namespace SimplePRNGs
{
public class LaggedFibRng
{
private const int k = 606; // Largest magnitude"-index"
private const int j = 273; // Other "-index"
private const long m = 4294967296; // 2^32
private readonly List<long> vals = null;
private long seed;
public LaggedFibRng(int seed)
{
vals = new List<long>();
for (int i = 0; i < k + 1; ++i)
vals.Add(i);
if (seed % 2 == 0) vals[0] = 11;
// Burn some values away
for (int ct = 0; ct < 1000; ++ct)
{
long dummy = Next();
}
} // ctor
public long Next()
{
// (a + b) mod n = [(a mod n) + (b mod n)] mod n
long left = vals[0] % m; // [x-big]
long right = vals[k - j] % m; // [x-other]
long sum = (left + right) % m; // prevent overflow
if (sum < 0)
seed = sum + m;
else
seed = sum;
vals.Insert(k + 1, seed); // Add new val at end
vals.RemoveAt(0); // Delete now irrelevant [0] val
return seed;
}
}
}
// "Numerical Computation 2: Methods, Software,
// and Analysis" by Christoph W. Ueberhuber
// Chapter 17 Random Numbers
// "The Art of Computer Programming Volume 2"
// "Seminumerical Algorithms Second Edition"
// "Chapter 3 RANDOM NUMBERS" by Donald E. Knuth
using System.Collections.Generic;
namespace SimplePRNGs
{
class PRNGs
{
private long AMxn1, AMxn;
private long AMyn1, AMyn;
private long AMk;
private long AMm = 34359738368;
private long LCG0z0, LCG0z1;
private long LCG1z0, LCG1z1;
private long LCG2z0, LCG2z1, LCG2z2;
private long MFG0z0, MFG0z1;
private readonly long LCG0m = 4294967296;
private readonly long LCG2m = 2147483647;
private readonly List<long> AMV = new();
public void SetSeedLCG0(long z0)
{
LCG0z0 = z0;
}
public void SetSeedLCG1(long z0)
{
LCG1z0 = z0;
}
public void SetSeedLCG2(long z0, long z1)
{
LCG2z0 = z0;
LCG2z1 = z1;
}
public long LCG0()
{
LCG0z1 = (69069 * LCG0z0) % LCG0m;
LCG0z0 = LCG0z1;
return LCG0z1;
}
public long LCG1()
{
LCG1z1 = (69069 * LCG1z0 + 1) % LCG0m;
LCG1z0 = LCG1z1;
return LCG1z1;
}
public long LCG2()
{
LCG2z2 = (1999 * LCG2z1 + 4444 * LCG2z0) % LCG2m;
LCG2z0 = LCG2z1;
LCG2z1 = LCG2z2;
return LCG2z2;
}
public void SetSeedMFG0(long z0, long z1)
{
MFG0z0 = z0;
MFG0z1 = z1;
}
public long MFG0()
{
long MFG0z2 = (MFG0z1 + MFG0z0) % LCG0m;
MFG0z0 = MFG0z1;
MFG0z1 = MFG0z2;
return MFG0z2;
}
public void ComputeNextXY()
{
AMxn1 = (3141592653 * AMxn + 2718281829) % AMm;
if (AMxn1 < 0)
AMxn1 += AMm;
AMyn1 = (2718281829 * AMyn + 3141592653) % AMm;
if (AMyn1 < 0)
AMyn1 += AMm;
}
public void AMSeed(long k, long x0, long y0)
{
long AMTxn1, AMTxn = x0;
AMxn = x0;
AMyn = y0;
AMk = k;
for (int i = 0; i < k; i++)
{
AMTxn1 = (3141592653 * AMTxn + 2718281829) % AMm;
if (AMTxn1 < 0)
AMTxn1 += AMm;
AMTxn = AMTxn1;
AMV.Add(AMTxn1);
}
}
public long AlgorithmM()
{
ComputeNextXY();
AMxn = AMxn1;
AMyn = AMyn1;
long j = (AMk * AMyn1) / AMm;
long r = AMV[(int)j];
AMV[(int)j] = AMxn1;
if (r < 0)
r += AMm;
return r;
}
}
}
namespace PointMassBallistics
{
public class TableEntry : IComparable<TableEntry>
{
public double range, elevationDegrees, elevationMinutes,
angleFallDegrees, angleFallMinutes, timeOfFlight,
strikingVelocity, maximumOrdinate;
public int CompareTo(TableEntry other)
{
if (elevationDegrees < other.elevationDegrees &&
elevationMinutes < other.elevationMinutes)
return -1;
if (elevationDegrees > other.elevationDegrees &&
elevationMinutes > other.elevationMinutes)
return +1;
return 0;
}
}
}
// Solves the following system of first order
// ordinary differential equations. Formulas
// are from "Elementary Numerical Analysis:
// An Algorithmic Approach 3rd Edition" by S.
// D. Conte & Carl de Boor (c) 1980 8.12 page 398.
// Extended from two to four equations.
// See https://apps.dtic.mil/dtic/tr/fulltext/u2/a439796.pdf
// Also view https://eugeneleeslover.com/USN-GUNS-AND-RANGE-TABLES/OP-770-1.html
using System;
using System.Collections.Generic;
namespace PointMassBallistics
{
class RungeKutta4
{
private double BC;
private readonly double g = 32.17405;
private readonly double[] GarveN = {
2, 3, 5, 3, 2, 1.7, 1.55 };
private readonly double[] log10K = {
5.66989 - 10, 2.77344 - 10, 6.80187 - 20,
2.98090 - 10, 6.11926 - 10, 7.09620 - 10, 7.60905 - 10 };
private readonly double[] K = new double[7];
private int zone;
static private double Density(double y)
{
return Math.Pow(10, -0.00001372 * y);
}
static private int ComputeIndex(double v)
{
int index;
if (v > 3600)
index = 6;
else if (v > 2600 && v <= 3600)
index = 5;
else if (v > 1800 && v <= 2600)
index = 4;
else if (v > 1370 && v <= 1800)
index = 3;
else if (v > 1230 && v <= 1370)
index = 2;
else if (v > 790 && v <= 1230)
index = 1;
else
index = 0;
return index;
}
public double MayevskiRetardation(double v, int zone)
{
// See Exterior Ballistics 1935 by Ernest Edward Herrmann
// Garve function
return K[zone] * Math.Pow(v, GarveN[zone]);
}
private double Vx(double syn, double vxn, double vyn)
{
double v = Math.Sqrt(vxn * vxn + vyn * vyn);
zone = ComputeIndex(v);
double E = Density(syn) * MayevskiRetardation(v, zone) / BC;
return -E * vxn / v;
}
private double Vy(double syn, double vxn, double vyn)
{
double v = Math.Sqrt(vxn * vxn + vyn * vyn);
zone = ComputeIndex(v);
double E = Density(syn) * MayevskiRetardation(v, zone) / BC;
return -E * vyn / v - g;
}
static private double Sx(double vxn)
{
return vxn;
}
static private double Sy(double vyn)
{
return vyn;
}
public void Solve(
double t0, double t1,
double vx0, double vy0,
double sx0, double sy0,
double BC, int nSteps, ref List<double> lt,
ref List<double> lvx, ref List<double> lvy,
ref List<double> lsx, ref List<double> lsy)
{
double k1, k2, k3, k4;
double l1, l2, l3, l4;
double m1, m2, m3, m4;
double n1, n2, n3, n4;
double h = (t1 - t0) / nSteps, tn = t0;
double vxn = vx0, vyn = vy0, sxn = sx0, syn = sy0;
int n = 1;
for (int i = 0; i < log10K.Length; i++)
K[i] = Math.Pow(10, log10K[i]);
this.BC = BC;
lt.Add(tn);
lvx.Add(vxn);
lvy.Add(vyn);
lsx.Add(sxn / 3);
lsy.Add(syn);
while (true)
{
tn = t0 + n * h;
k1 = h * Vx(syn, vxn, vyn);
l1 = h * Vy(syn, vxn, vyn);
m1 = h * Sx(vxn);
n1 = h * Sy(vyn);
k2 = h * Vx(syn + 0.5 * n1, vxn + 0.5 * k1, vyn + 0.5 * l1);
l2 = h * Vy(syn + 0.5 * n1, vxn + 0.5 * k1, vyn + 0.5 * l1);
m2 = h * Sx(vxn + 0.5 * m1);
n2 = h * Sy(vyn + 0.5 * n1);
k3 = h * Vx(syn + 0.5 * n2, vxn + 0.5 * k2, vyn + 0.5 * l2);
l3 = h * Vy(syn + 0.5 * n2, vxn + 0.5 * k2, vyn + 0.5 * l2);
m3 = h * Sx(vxn + 0.5 * m2);
n3 = h * Sy(vyn + 0.5 * n2);
k4 = h * Vx(syn + n3, vxn + k3, vyn + l3);
l4 = h * Vy(syn + n3, vxn + k3, vyn + l3);
m4 = h * Sx(vxn + m3);
n4 = h * Sy(vyn + n3);
vxn = vx0 + (k1 + 2 * k2 + 2 * k3 + k4) / 6.0;
vyn = vy0 + (l1 + 2 * l2 + 2 * l3 + l4) / 6.0;
sxn = sx0 + (m1 + 2 * m2 + 2 * m3 + m4) / 6.0;
syn = sy0 + (n1 + 2 * n2 + 2 * n3 + n4) / 6.0;
vx0 = vxn;
vy0 = vyn;
sx0 = sxn;
sy0 = syn;
n++;
lt.Add(tn);
lvx.Add(vxn);
lvy.Add(vyn);
lsx.Add(sxn / 3);
lsy.Add(syn);
if (syn <= 1.0e-2)
break;
}
}
}
}
// Solves the following system of first order
// ordinary differential equations. Formulas
// are from "Elementary Numerical Analysis:
// An Algorithmic Approach 3rd Edition" by S.
// D. Conte & Carl de Boor (c) 1980 8.12 page 398.
// Extended from two to four equations.
// See "Ordinary Differential Equations
// from Calculus to Dynamical Systems"
// by Virginia W. Noonburg for the exact
// solution, view pages 167 - 170.
using System.Collections.Generic;
namespace DoubleSpring
{
class RungeKutta4
{
private double omega2;
static private double Dx1(double x3)
{
return x3;
}
static private double Dx2(double x4)
{
return x4;
}
private double Dx3(double x1, double x2)
{
return -2 * omega2 * x1 + omega2 * x2;
}
private double Dx4(double x1, double x2)
{
return omega2 * x1 - omega2 * x2;
}
public void Solve(
double t0, double t1,
double x10, double x20,
double x30, double x40,
double omega, int nSteps, ref List<double> lt,
ref List<double> lx1, ref List<double> lx2,
ref List<double> lx3, ref List<double> lx4)
{
double k1, k2, k3, k4;
double l1, l2, l3, l4;
double m1, m2, m3, m4;
double n1, n2, n3, n4;
double h = (t1 - t0) / nSteps, tn = t0;
double x1n = x10, x2n = x20, x3n = x30, x4n = x40;
int n = 1;
omega2 = omega * omega;
lt.Add(tn);
lx1.Add(x1n);
lx2.Add(x2n);
lx3.Add(x3n);
lx4.Add(x4n);
while (tn <= t1)
{
tn = t0 + n * h;
k1 = h * Dx1(x3n);
l1 = h * Dx2(x4n);
m1 = h * Dx3(x1n, x2n);
n1 = h * Dx4(x1n, x2n);
k2 = h * Dx1(x3n + 0.5 * k1);
l2 = h * Dx2(x4n + 0.5 * l1);
m2 = h * Dx3(x1n + 0.5 * m1, x2n + 0.5 * n1);
n2 = h * Dx4(x1n + 0.5 * m1, x2n + 0.5 * n1);
k3 = h * Dx1(x3n + 0.5 * k2);
l3 = h * Dx2(x4n + 0.5 * l2);
m3 = h * Dx3(x1n + 0.5 * m2, x2n + 0.5 * n2);
n3 = h * Dx4(x1n + 0.5 * m2, x2n + 0.5 * n2);
k4 = h * Dx1(x3n + k3);
l4 = h * Dx2(x4n + l3);
m4 = h * Dx3(x1n + m3, x2n + n3);
n4 = h * Dx4(x1n + m3, x2n + n3);
x1n = x10 + (k1 + 2 * k2 + 2 * k3 + k4) / 6.0;
x2n = x20 + (l1 + 2 * l2 + 2 * l3 + l4) / 6.0;
x3n = x30 + (m1 + 2 * m2 + 2 * m3 + m4) / 6.0;
x4n = x40 + (n1 + 2 * n2 + 2 * n3 + n4) / 6.0;
x10 = x1n;
x20 = x2n;
x30 = x3n;
x40 = x4n;
n++;
lt.Add(tn);
lx1.Add(x1n);
lx2.Add(x2n);
lx3.Add(x3n);
lx4.Add(x4n);
}
}
}
}
// Two-tank Mixing Problem
// x'(t) = f(t, x, y)
// y'(t) = g(t, x, y)
// x(0) = 5, y(0) = 1
// See “Ordinary Differential Equations
// from Calculus to Dynamical Systems”
// by Virginia W. Noonburg for the exact
// solution, view pages 165 – 167. Also
// see “Elementary Numerical Analysis:
// An Algorithmic Approach Third
// Edition” by S. D. Conte and Carl de
// Boor pages 398 – 491 for the Runge-
// Kutta-4 equations.
#include "RK4System.h"
void RK4System::Solve(
real t0, real t1,
real x0, real z0,
real(*f)(real, real, real),
real(*g)(real, real, real),
int nSteps, vector<real>& tv,
vector<real>& xv, vector<real>& yv)
{
real k1, k2, k3, k4, l1, l2, l3, l4;
real h, tn, xn, yn, xn1, yn1;
h = (t1 - t0) / nSteps;
tn = t0;
xn = x0;
yn = z0;
tv.push_back(tn);
xv.push_back(xn);
yv.push_back(yn);
for (unsigned int n = 1; n <= nSteps; n++)
{
tn = t0 + n * h;
k1 = h * f(tn, xn, yn);
l1 = h * g(tn, xn, yn);
k2 = h * f(tn + 0.5 * h, xn + 0.5 * k1, yn + 0.5 * l1);
l2 = h * g(tn + 0.5 * h, xn + 0.5 * k1, yn + 0.5 * l1);
k3 = h * f(tn + 0.5 * h, xn + 0.5 * k2, yn + 0.5 * l2);
l3 = h * g(tn + 0.5 * h, xn + 0.5 * k2, yn + 0.5 * l2);
k4 = h * f(tn + h, xn + k3, yn + l3);
l4 = h * g(tn + h, xn + k3, yn + l3);
xn1 = xn + (k1 + 2 * k2 + 2 * k3 + k4) / 6.0;
yn1 = yn + (l1 + 2 * l2 + 2 * l3 + l4) / 6.0;
xn = xn1;
yn = yn1;
tv.push_back(tn);
xv.push_back(xn);
yv.push_back(yn);
}
}
#pragma once
// Solves the following system of first order
// ordinary differential equations. Formulas
// are from "Elementary Numerical Analysis:
// An Algorithmic Approach" by S. D. Conte &
// Carl de Boor (c) 1980 8.12 page 398.
// x' = f(t, x, y)
// y' = g(t, x, y)
#include <vector>
using namespace std;
typedef long double real;
class RK4System
{
public:
static void Solve(
real t0, real t1,
real x0, real y0,
real(*f)(real, real, real),
real(*g)(real, real, real),
int nSteps, vector<real>& tv,
vector<real>& xv, vector<real>& yv);
};
#include "RungeKutta4.h"
#include <iomanip>
#include <iostream>
using namespace std;
int main()
{
int nSteps;
vector<Point> solution1, solution2;
vector<Point> solution3, solution4;
while (true)
{
cout << "n-steps (zero to exit app) = ";
cin >> nSteps;
if (nSteps == 0)
break;
cout << "n\t\txn\tyn2p(xn)\tyn4p(xn)\tEuler\t\tTaylor\t\tDFDX2\t\tDFDX4\r\n";
RungeKutta4::ComputeEuler(nSteps, 1, 2, -1, solution1);
RungeKutta4::ComputeRK2(nSteps, 1, 2, -1, solution2);
RungeKutta4::ComputeRK4(nSteps, 1, 2, -1, solution3);
RungeKutta4::ComputeTaylor(nSteps, 1, 2, -1, solution4);
for (int i = 0; i <= nSteps; i++)
{
POINT pt1 = solution1[i];
POINT pt2 = solution2[i];
POINT pt3 = solution3[i];
POINT pt4 = solution4[i];
cout << setw(2) << pt1.n << "\t";
cout << setprecision(8) << fixed << pt1.x << "\t";
cout << setprecision(8) << fixed << pt3.y << "\t";
cout << setprecision(8) << fixed << pt4.y << "\t";
cout << setprecision(8) << fixed << pt1.y << "\t";
cout << setprecision(8) << fixed << pt2.y << "\t";
cout << setprecision(8) << fixed << pt2.derivative << "\t";
cout << setprecision(8) << fixed << pt3.derivative << "\r\n";
}
}
}
#pragma once
#include <vector>
using namespace std;
typedef struct Point
{
int n;
long double derivative, x, y;
} POINT, *PPOINT;
class RungeKutta4
{
private:
static long double FX(
long double x,
long double y);
static long double DFDX(
long double x,
long double y,
long double yp);
public:
static void ComputeEuler(
int nSteps,
long double x0,
long double x1,
long double y0,
vector<Point>& pts);
static void ComputeRK2(
int nSteps,
long double x0,
long double x1,
long double y0,
vector<Point>& pts);
static void ComputeRK4(
int nSteps,
long double x0,
long double x1,
long double y0,
vector<Point>& pts);
static void ComputeTaylor(
int nSteps,
long double x0,
long double x1,
long double y0,
vector<Point>& pts);
};
// Translated from FORTRAN 77 source code found
// in "Elementary Numerical Analysis:
// An Algorithmic Approach" by S. D. Conte and Carl de
// Boor. Translator: James Pate Williams, Jr. (c)
// August 20, 2023
#include "RungeKutta4.h"
long double RungeKutta4::FX(
long double x, long double y)
{
long double x2 = x * x;
long double term1 = 1.0 / x2;
long double term2 = -y / x;
long double term3 = -y * y;
return term1 + term2 + term3;
}
long double RungeKutta4::DFDX(
long double x,
long double y,
long double yp)
{
long double x3 = x * x * x;
long double term1 = -2.0 / x3;
long double term2 = -yp / x;
long double term3 = y / x / x;
long double term4 = -2.0 * y * yp;
return term1 + term2 + term3 + term4;
}
void RungeKutta4::ComputeEuler(
int nSteps,
long double x0,
long double x1,
long double y0,
vector<Point>& pts)
{
Point pt0 = {};
pt0.n = 0;
pt0.x = x0;
pt0.y = y0;
double derivative = FX(x0, y0);
pt0.derivative = derivative;
pts.push_back(pt0);
if (nSteps == 1)
return;
long double h = (x1 - x0) / nSteps;
long double xn = x0;
long double yn = y0;
for (int n = 1; n <= nSteps; n++)
{
xn = x0 + n * h;
long double yn = y0 + h * FX(xn, y0);
derivative = FX(xn, y0);
pt0.n = n;
pt0.derivative = derivative;
pt0.x = xn;
pt0.y = yn;
pts.push_back(pt0);
y0 = yn;
}
}
void RungeKutta4::ComputeRK2(
int nSteps,
long double x0,
long double x1,
long double y0,
vector<Point>& pts)
{
Point pt0 = {};
pt0.n = 0;
pt0.x = x0;
pt0.y = y0;
double derivative = FX(x0, y0);
pt0.derivative = derivative;
pts.push_back(pt0);
if (nSteps == 1)
return;
long double h = (x1 - x0) / nSteps;
long double xn = x0;
long double yn = y0;
for (int n = 1; n <= nSteps; n++)
{
long double k1 = h * FX(xn, yn);
long double k2 = h * FX(xn + h, yn + k1);
yn += 0.5 * (k1 + k2);
xn = x0 + n * h;
derivative = FX(xn, yn);
pt0.n = n;
pt0.derivative = derivative;
pt0.x = xn;
pt0.y = yn;
pts.push_back(pt0);
}
}
void RungeKutta4::ComputeRK4(
int nSteps,
long double x0,
long double x1,
long double y0,
vector<Point>& pts)
{
Point pt0 = {};
pt0.n = 0;
pt0.x = x0;
pt0.y = y0;
double derivative = FX(x0, y0);
pt0.derivative = derivative;
pts.push_back(pt0);
if (nSteps == 1)
return;
long double h = (x1 - x0) / nSteps;
long double xn = x0;
long double yn = y0;
for (int n = 1; n <= nSteps; n++)
{
long double k1 = h * FX(xn, yn);
long double k2 = h * FX(xn + h / 2, yn + k1 / 2);
long double k3 = h * FX(xn + h / 2, yn + k2 / 2);
long double k4 = h * FX(xn + h, yn + k3);
xn = x0 + n * h;
yn = yn + (k1 + k2 + k3 + k4) / 6.0;
derivative = FX(xn, yn);
pt0.n = n;
pt0.derivative = derivative;
pt0.x = xn;
pt0.y = yn;
pts.push_back(pt0);
}
}
void RungeKutta4::ComputeTaylor(
int nSteps,
long double x0,
long double x1,
long double y0,
vector<Point>& pts)
{
Point pt0 = {};
pt0.n = 0;
pt0.x = x0;
pt0.y = y0;
double derivative = FX(x0, y0);
pt0.derivative = derivative;
pts.push_back(pt0);
if (nSteps == 1)
return;
long double h = (x1 - x0) / nSteps;
long double xn = x0;
long double yn = y0;
long double yp = FX(x0, y0);
for (int n = 1; n <= nSteps; n++)
{
xn = x0 + n * h;
yp = FX(xn, y0);
long double yn = y0 + h * (FX(xn, y0) + h * DFDX(xn, y0, yp) / 2);
derivative = FX(xn, y0);
pt0.n = n;
pt0.derivative = derivative;
pt0.x = xn;
pt0.y = yn;
pts.push_back(pt0);
y0 = yn;
}
}
using System;
using System.Windows.Forms;
namespace NumericalDifferentiation
{
public partial class Form1 : Form
{
public Form1()
{
InitializeComponent();
comboBox1.SelectedIndex = 0;
}
private double F1(double x)
{
return x * (x * (x - 2) + 4) - 1;
}
private double F2(double x)
{
return x * Math.Exp(-x) + 1;
}
private double F3(double x)
{
return x * x * Math.Sin(x) + 1;
}
private double DF1(double x)
{
return 3 * x * x - 4 * x + 4;
}
private double DF2(double x)
{
return Math.Exp(-x) - x * Math.Exp(-x);
}
private double DF3(double x)
{
return 2 * x * Math.Sin(x) + x * x * Math.Cos(x);
}
private double DDF1(double x)
{
return 6 * x - 4;
}
private double DDF2(double x)
{
return -2 * Math.Exp(-x) + x * Math.Exp(-x);
}
private double DDF3(double x)
{
return
2 * Math.Sin(x) + 2 * x * Math.Cos(x) +
2 * x * Math.Cos(x) - x * x * Math.Sin(x);
}
private void button1_Click(object sender, EventArgs e)
{
double a = (double)numericUpDown1.Value;
double b = (double)numericUpDown2.Value;
int n = (int)numericUpDown3.Value;
double deriv1c, deriv1f;
double deriv2c, deriv2f;
double x = a, y, u, v;
double error1c, error1f, error2c, error2f;
double h = (b - a) / n;
if (comboBox1.SelectedIndex == 0)
textBox1.Text += "x\t f'(x) e f'(x) 1 f'(x) 2 " +
"error1 error2 f''(x)e f''(x)1 f''(x)2 " +
"error1 error2\r\n";
else if (comboBox1.SelectedIndex == 1)
textBox1.Text += "x\t g'(x) e g'(x) 1 g'(x) 2 " +
"error1 error2 g''(x)e g''(x)1 g''(x)2 " +
"error1 error2\r\n";
else if (comboBox1.SelectedIndex == 2)
textBox1.Text += "x\t h'(x) e h'(x) 1 h'(x) 2 " +
"error1 error2 h''(x)e h''(x)1 h''(x)2 " +
"error1 error2\r\n";
while (x <= b)
{
switch (comboBox1.SelectedIndex)
{
case 0:
y = F1(x);
u = DF1(x);
v = DDF1(x);
deriv1c = Differentiation.Derivative1sc(x, h, F1);
deriv1f = Differentiation.Derivative1sf(x, h, F1);
deriv2c = Differentiation.Derivative2sc(x, h, F1);
deriv2f = Differentiation.Derivative2sf(x, h, F1);
error1c = Math.Abs(deriv1c - u);
error1f = Math.Abs(deriv1f - u);
error2c = Math.Abs(deriv2c - v);
error2f = Math.Abs(deriv2f - v);
textBox1.Text +=
x.ToString("F5").PadLeft(8) + " " +
u.ToString("F5").PadLeft(8) + " " +
deriv1c.ToString("F5").PadLeft(8) + " " +
deriv1f.ToString("F5").PadLeft(8) + " " +
error1c.ToString("F5").PadLeft(8) + " " +
error1f.ToString("F5").PadLeft(8) + " " +
v.ToString("F5").PadLeft(8) + " " +
deriv2c.ToString("F5").PadLeft(8) + " " +
deriv2f.ToString("F5").PadLeft(8) + " " +
error2c.ToString("F5").PadLeft(8) + " " +
error2f.ToString("F5").PadLeft(8) + "\r\n";
break;
case 1:
y = F2(x);
u = DF2(x);
v = DDF2(x);
deriv1c = Differentiation.Derivative1sc(x, h, F2);
deriv1f = Differentiation.Derivative1sf(x, h, F2);
deriv2c = Differentiation.Derivative2sc(x, h, F2);
deriv2f = Differentiation.Derivative2sf(x, h, F2);
error1c = Math.Abs(deriv1c - u);
error1f = Math.Abs(deriv1f - u);
error2c = Math.Abs(deriv2c - v);
error2f = Math.Abs(deriv2f - v);
textBox1.Text +=
x.ToString("F5").PadLeft(8) + " " +
u.ToString("F5").PadLeft(8) + " " +
deriv1c.ToString("F5").PadLeft(8) + " " +
deriv1f.ToString("F5").PadLeft(8) + " " +
error1c.ToString("F5").PadLeft(8) + " " +
error1f.ToString("F5").PadLeft(8) + " " +
v.ToString("F5").PadLeft(8) + " " +
deriv2c.ToString("F5").PadLeft(8) + " " +
deriv2f.ToString("F5").PadLeft(8) + " " +
error2c.ToString("F5").PadLeft(8) + " " +
error2f.ToString("F5").PadLeft(8) + "\r\n";
break;
case 2:
y = F3(x);
u = DF3(x);
v = DDF3(x);
deriv1c = Differentiation.Derivative1sc(x, h, F3);
deriv1f = Differentiation.Derivative1sf(x, h, F3);
deriv2c = Differentiation.Derivative2sc(x, h, F3);
deriv2f = Differentiation.Derivative2sf(x, h, F3);
error1c = Math.Abs(deriv1c - u);
error1f = Math.Abs(deriv1f - u);
error2c = Math.Abs(deriv2c - v);
error2f = Math.Abs(deriv2f - v);
textBox1.Text +=
x.ToString("F5").PadLeft(8) + " " +
u.ToString("F5").PadLeft(8) + " " +
deriv1c.ToString("F5").PadLeft(8) + " " +
deriv1f.ToString("F5").PadLeft(8) + " " +
error1c.ToString("F5").PadLeft(8) + " " +
error1f.ToString("F5").PadLeft(8) + " " +
v.ToString("F5").PadLeft(8) + " " +
deriv2c.ToString("F5").PadLeft(8) + " " +
deriv2f.ToString("F5").PadLeft(8) + " " +
error2c.ToString("F5").PadLeft(8) + " " +
error2f.ToString("F5").PadLeft(8) + "\r\n";
break;
default:
break;
}
x += h;
}
textBox1.Text += comboBox1.SelectedItem;
}
private void button2_Click(object sender, EventArgs e)
{
textBox1.Text = string.Empty;
}
}
}
// Derivative formulas translated from FORTRAN 77
// source code found in "Elementary Numerical Analysis:
// An Algorithmic Approach" by S. D. Conte and Carl de
// Boor. Translator: James Pate Williams, Jr. (c)
// August 18 - 19, 2023
using System;
namespace NumericalDifferentiation
{
// Numerically compute first and second
// ordinary derivatives
class Differentiation
{
static public double Derivative1sc(
double a, double h, Func<double, double> f)
{
return (f(a + h) - f(a - h)) / (h + h);
}
static public double Derivative1sf(
double a, double h, Func<double, double> f)
{
return (-3.0 * f(a) + 4.0 * f(a + h) - f(a + h + h)) / (h + h);
}
static public double Derivative2sc(
double a, double h, Func<double, double> f)
{
return (f(a - h) - 2.0 * f(a) + f(a + h)) / (h + h);
}
static public double Derivative2sf(
double a, double h, Func<double, double> f)
{
return (f(a) - 2.0 * f(a + h) + f(a + h + h)) / (h + h);
}
}
}
using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using System.Windows.Forms;
namespace AdamsMoultonMethod1
{
public partial class Form1 : Form
{
public Form1()
{
InitializeComponent();
}
private double Solution(double x)
{
return Math.Exp(x) - 1.0 - x;
}
private void button1_Click(object sender, EventArgs e)
{
int numberSteps = (int)numericUpDown1.Value;
AdamsMoulton.AdamsMoultonMethod(
numberSteps, Solution, textBox1);
}
private void button2_Click(object sender, EventArgs e)
{
textBox1.Text = string.Empty;
}
}
}
// Adams-Moulton Method translated from FORTRAN 77
// source code found in "Elementary Numerical Analysis:
// An Algorithmic Approach" by S. D. Conte and Carl de
// Boor. Translator: James Pate Williams, Jr. (c)
// August 18, 2023
using System;
using System.Windows.Forms;
namespace AdamsMoultonMethod1
{
public struct Point
{
public double x, y;
}
class AdamsMoulton
{
static private string FormatNumber(double x)
{
return x.ToString("E11").PadLeft(12);
}
static private string FormatLine(
int n, double x, double y, double delta, double error)
{
return
n.ToString("D2") + " " +
FormatNumber(x) + " " +
FormatNumber(y) + " " +
FormatNumber(delta) + " " +
FormatNumber(error) + "\r\n";
}
static public void AdamsMoultonMethod(
int numberSteps, Func<double, double> solution,
TextBox textBox)
{
double h = 1.0 / numberSteps;
double xn = 0.0, yn = 0.0, yofxn;
double xBegin = 0.0, yBegin = 0.0;
double fnPredict, ynPredict;
double delta = 0.0, error = 0.0;
double[] f = new double[numberSteps];
int n = 0;
textBox.Text += "n\txn\tyn\tdelta\terror\r\n";
textBox.Text += FormatLine(n, xn, yn, delta, error);
f[1] = xBegin + yBegin;
for (n = 1; n <= 3; n++)
{
xn = xBegin + n * h;
yn = solution(xn);
f[n + 1] = xn + yn;
textBox.Text += FormatLine(n, xn, yn, delta, error);
}
for (n = 4; n <= numberSteps; n++)
{
ynPredict =
yn + h / 24.0 * (55.0 * f[4] - 59.0 * f[3] +
37.0 * f[2] - 9.0 * f[1]);
xn = xBegin + n * h;
fnPredict = xn + ynPredict;
yn += h / 24.0 * (9.0 * fnPredict + 19.0 * f[4]
- 5.0 * f[3] + f[2]);
delta = (yn - ynPredict) / 14.0;
f[1] = f[2];
f[2] = f[3];
f[3] = f[4];
f[4] = xn + yn;
yofxn = solution(xn);
error = yn - yofxn;
textBox.Text += FormatLine(n, xn, yn, delta, error);
}
}
}
}
Below are three screenshots of two methods of calculating the determinant of a matrix, namely the Bareiss Algorithm and Gaussian Elimination:
using System;
using System.Diagnostics;
using System.Windows.Forms;
namespace MatrixInverseComparison
{
public partial class MainForm : Form
{
public MainForm()
{
InitializeComponent();
}
static private string FormatNumber(double x)
{
string result = string.Empty;
if (x > 0)
result += x.ToString("F5").PadLeft(10);
else
result += x.ToString("F5").PadLeft(10);
return result;
}
private void MultiplyPrintMatricies(
double[,] A, double[,] B, int n)
{
double[,] I = new double[n, n];
textBox1.Text += "Matrix Product:\r\n";
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
double sum = 0.0;
for (int k = 0; k < n; k++)
sum += A[i, k] * B[k, j];
I[i, j] = sum;
}
}
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
textBox1.Text += FormatNumber(I[i, j]) + " ";
}
textBox1.Text += "\r\n";
}
}
private void PrintMatrix(string title, double[,] A, int n)
{
textBox1.Text += title + "\r\n";
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
textBox1.Text += FormatNumber(A[i, j]) + " ";
}
textBox1.Text += "\r\n";
}
}
private void Compute(double[,] MI, int n)
{
double determinantGE = 1;
double[] b = new double[n];
double[] x = new double[n];
double[,] MB1 = new double[n, n];
double[,] MB2 = new double[n, n];
double[,] MG = new double[n, n];
double[,] MS = new double[n, n];
double[,] IG = new double[n, n];
double[,] IS = new double[n, n];
int[] pivot = new int[n];
Stopwatch sw = new();
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
MB1[i, j] = MB2[i, j] = MG[i, j] = MS[i, j] = MI[i, j];
}
PrintMatrix("Initial Matrix: ", MI, n);
sw.Start();
int flag = DirectMethods.Factor(MG, n, pivot);
if (flag != 0)
{
for (int i = 0; i < n; i++)
determinantGE *= MG[i, i];
determinantGE *= flag;
}
sw.Stop();
PrintMatrix("Gaussian Elimination Final:", MG, n);
for (int i = 0; i < n; i++)
b[i] = 0;
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
b[j] = 1.0;
DirectMethods.Substitute(MG, b, x, n, pivot);
for (int k = 0; k < n; k++)
IG[k, j] = x[k];
b[j] = 0.0;
}
}
PrintMatrix("Gaussian Elimination Inverse:", IG, n);
textBox1.Text += "Determinant: " +
FormatNumber(determinantGE) + "\r\n";
MultiplyPrintMatricies(IG, MI, n);
textBox1.Text += "Runtime (MS) = " +
sw.ElapsedMilliseconds + "\r\n";
sw.Start();
double determinant1 = BareissAlgorithm.Determinant1(MB1, n);
sw.Stop();
PrintMatrix("Bareiss Algorithm Final 1:", MB1, n);
textBox1.Text += "Determinant: " +
FormatNumber(determinant1) + "\r\n";
textBox1.Text += "Runtime (MS) = " +
sw.ElapsedMilliseconds + "\r\n";
sw.Start();
double determinant2 = BareissAlgorithm.Determinant2(MB2, n);
sw.Stop();
PrintMatrix("Bareiss Algorithm Final 2:", MB2, n);
textBox1.Text += "Determinant: " +
FormatNumber(determinant2) + "\r\n";
textBox1.Text += "Runtime (MS) = " +
sw.ElapsedMilliseconds + "\r\n";
}
private void button1_Click(object sender, EventArgs e)
{
int n = (int)numericUpDown1.Value;
int seed = (int)numericUpDown2.Value;
double[,] A = new double[n, n];
Random random = new Random(seed);
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
double q = n * random.NextDouble();
while (q == 0.0)
q = n * random.NextDouble();
A[i, j] = q;
}
}
Compute(A, n);
}
private void button2_Click(object sender, EventArgs e)
{
textBox1.Text = string.Empty;
}
}
}
using System;
namespace MatrixInverseComparison
{
public class DirectMethods
{
// Substitute and Factor translated from FORTRAN 77
// source code found in "Elementary Numerical Analysis:
// An Algorithmic Approach" by S. D. Conte and Carl de
// Boor. Translator: James Pate Williams, Jr. (c)
// August 14 - 17, 2023
static public void Substitute(
double[,] w,
double[] b,
double[] x,
int n,
int[] pivot)
{
double sum;
int i, j, n1 = n - 1;
if (n == 1)
{
x[0] = b[0] / w[0, 0];
return;
}
// forward substitution
x[0] = b[pivot[0]];
for (i = 0; i < n; i++)
{
for (j = 0, sum = 0.0; j < i; j++)
sum += w[i, j] * x[j];
x[i] = b[pivot[i]] - sum;
}
// backward substitution
x[n1] /= w[n1, n1];
for (i = n - 2; i >= 0; i--)
{
for (j = i + 1, sum = 0.0; j < n; j++)
sum += w[i, j] * x[j];
x[i] = (x[i] - sum) / w[i, i];
}
}
// Factor returns +1 if an even number of exchanges
// Factor returns -1 if an odd number of exchanges
// Factor retrurn 0 if matrix is singular
static public int Factor(
double[,] w, int n, int[] pivot)
// returns 0 if matrix is singular
{
double awikod, col_max, ratio, row_max, temp;
double[] d = new double[n];
int flag = 1, i, i_star, j, k;
for (i = 0; i < n; i++)
{
pivot[i] = i;
row_max = 0;
for (j = 0; j < n; j++)
row_max = Math.Max(row_max, Math.Abs(w[i, j]));
if (row_max == 0)
{
flag = 0;
row_max = 1;
}
d[i] = row_max;
}
if (n <= 1)
return flag;
// factorization
for (k = 0; k < n - 1; k++)
{
// determine pivot row the row i_star
col_max = Math.Abs(w[k, k]) / d[k];
i_star = k;
for (i = k + 1; i < n; i++)
{
awikod = Math.Abs(w[i, k]) / d[i];
if (awikod > col_max)
{
col_max = awikod;
i_star = i;
}
}
if (col_max == 0)
flag = 0;
else
{
if (i_star > k)
{
// make k the pivot row by
// interchanging with i_star
flag *= -1;
i = pivot[i_star];
pivot[i_star] = pivot[k];
pivot[k] = i;
temp = d[i_star];
d[i_star] = d[k];
d[k] = temp;
for (j = 0; j < n; j++)
{
temp = w[i_star, j];
w[i_star, j] = w[k, j];
w[k, j] = temp;
}
}
// eliminate x[k]
for (i = k + 1; i < n; i++)
{
w[i, k] /= w[k, k];
ratio = w[i, k];
for (j = k + 1; j < n; j++)
w[i, j] -= ratio * w[k, j];
}
}
}
if (w[n - 1, n - 1] == 0)
flag = 0;
return flag;
}
}
}
namespace MatrixInverseComparison
{
// One implementation is based on https://en.wikipedia.org/wiki/Bareiss_algorithm
// Another perhaps better implementation is found on the following webpage
// https://cs.stackexchange.com/questions/124759/determinant-calculation-bareiss-vs-gauss-algorithm
class BareissAlgorithm
{
static public double Determinant1(double[,] M, int n)
{
double M00;
for (int k = 0; k < n; k++)
{
if (k - 1 == -1)
M00 = 1;
else
M00 = M[k - 1, k - 1];
for (int i = k + 1; i < n; i++)
{
for (int j = k + 1; j < n; j++)
{
M[i, j] = (
M[i, j] * M[k, k] -
M[i, k] * M[k, j]) / M00;
}
}
}
return M[n - 1, n - 1];
}
static public double Determinant2(double[,] A, int dim)
{
if (dim <= 0)
{
return 0;
}
double sign = 1;
for (int k = 0; k < dim - 1; k++)
{
//Pivot - row swap needed
if (A[k, k] == 0)
{
int m;
for (m = k + 1; m < dim; m++)
{
if (A[m, k] != 0)
{
double[] tempRow = new double[dim];
for (int i = 0; i < dim; i++)
tempRow[i] = A[m, i];
for (int i = 0; i < dim; i++)
A[m, i] = A[k, i];
for (int i = 0; i < dim; i++)
A[k, i] = tempRow[i];
sign = -sign;
break;
}
}
//No entries != 0 found in column k -> det = 0
if (m == dim)
{
return 0;
}
}
//Apply formula
for (int i = k + 1; i < dim; i++)
{
for (int j = k + 1; j < dim; j++)
{
A[i, j] = A[k, k] * A[i, j] - A[i, k] * A[k, j];
if (k >= 1)
{
A[i, j] /= A[k - 1, k - 1];
}
}
}
}
return sign * A[dim - 1, dim - 1];
}
}
}
Numerical Explorations 