// 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
Tensor Calculus Examples I © Monday, May, 4, 2026 by James Pate Williams, Jr.
Blog Entry © Sunday, April 26, 2026, by James Pate Williams, Jr. Butterfly Curve
Butterfly curve (transcendental) – Wikipedia
// ButterflyEquation.cpp : Defines the entry point for the application.
//
#include "framework.h"
#include "ButterflyEquation.h"
#include <vector>
#define MAX_LOADSTRING 100
typedef struct tagPoint3d
{
double t, x, y;
} Point3d, * PPoint3d;
// Global Variables:
HINSTANCE hInst; // current instance
WCHAR szTitle[MAX_LOADSTRING]; // The title bar text
WCHAR szWindowClass[MAX_LOADSTRING]; // the main window class name
std::vector<Point3d> points;
// Forward declarations of functions included in this code module:
ATOM MyRegisterClass(HINSTANCE hInstance);
BOOL InitInstance(HINSTANCE, int);
LRESULT CALLBACK WndProc(HWND, UINT, WPARAM, LPARAM);
INT_PTR CALLBACK About(HWND, UINT, WPARAM, LPARAM);
int APIENTRY wWinMain(_In_ HINSTANCE hInstance,
_In_opt_ HINSTANCE hPrevInstance,
_In_ LPWSTR lpCmdLine,
_In_ int nCmdShow)
{
UNREFERENCED_PARAMETER(hPrevInstance);
UNREFERENCED_PARAMETER(lpCmdLine);
// TODO: Place code here.
// Initialize global strings
LoadStringW(hInstance, IDS_APP_TITLE, szTitle, MAX_LOADSTRING);
LoadStringW(hInstance, IDC_BUTTERFLYEQUATION, szWindowClass, MAX_LOADSTRING);
MyRegisterClass(hInstance);
// Perform application initialization:
if (!InitInstance (hInstance, nCmdShow))
{
return FALSE;
}
HACCEL hAccelTable = LoadAccelerators(hInstance, MAKEINTRESOURCE(IDC_BUTTERFLYEQUATION));
MSG msg;
// Main message loop:
while (GetMessage(&msg, nullptr, 0, 0))
{
if (!TranslateAccelerator(msg.hwnd, hAccelTable, &msg))
{
TranslateMessage(&msg);
DispatchMessage(&msg);
}
}
return (int) msg.wParam;
}
//
// FUNCTION: MyRegisterClass()
//
// PURPOSE: Registers the window class.
//
ATOM MyRegisterClass(HINSTANCE hInstance)
{
WNDCLASSEXW wcex = { 0 };
wcex.cbSize = sizeof(WNDCLASSEX);
wcex.style = CS_HREDRAW | CS_VREDRAW;
wcex.lpfnWndProc = WndProc;
wcex.cbClsExtra = 0;
wcex.cbWndExtra = 0;
wcex.hInstance = hInstance;
wcex.hIcon = LoadIcon(hInstance, MAKEINTRESOURCE(IDI_BUTTERFLYEQUATION));
wcex.hCursor = LoadCursor(nullptr, IDC_ARROW);
wcex.hbrBackground = (HBRUSH)(COLOR_WINDOW+1);
wcex.lpszMenuName = MAKEINTRESOURCEW(IDC_BUTTERFLYEQUATION);
wcex.lpszClassName = szWindowClass;
wcex.hIconSm = LoadIcon(wcex.hInstance, MAKEINTRESOURCE(IDI_SMALL));
return RegisterClassExW(&wcex);
}
//
// FUNCTION: InitInstance(HINSTANCE, int)
//
// PURPOSE: Saves instance handle and creates main window
//
// COMMENTS:
//
// In this function, we save the instance handle in a global variable and
// create and display the main program window.
//
BOOL InitInstance(HINSTANCE hInstance, int nCmdShow)
{
hInst = hInstance; // Store instance handle in our global variable
HWND hWnd = CreateWindowW(szWindowClass, szTitle, WS_OVERLAPPEDWINDOW,
CW_USEDEFAULT, 0, CW_USEDEFAULT, 0, nullptr, nullptr, hInstance, nullptr);
if (!hWnd)
{
return FALSE;
}
ShowWindow(hWnd, nCmdShow);
UpdateWindow(hWnd);
return TRUE;
}
void CreateButterflyGraphPoints()
{
double p = 4.0 * atan(1.0);
double h = 12.0 * p / 1024.0;
double t = 0.0;
for (int i = 0; i <= 1024; i++)
{
Point3d pt = { 0 };
double c = 2.0 * cos(4.0 * t);
double d = pow(sin(t / 12.0), 5.0);
double x = sin(t) * (exp(cos(t)) - c - d);
double y = cos(t) * (exp(cos(t)) - c - d);
pt.t = t;
pt.x = x;
pt.y = y;
points.push_back(pt);
t += h;
}
}
static void FindMinMax(
double& xMin, double& xMax,
double& yMin, double& yMax)
{
// uses global 2D double points structure
xMin = yMin = DBL_MAX;
xMax = yMax = DBL_MIN;
for (size_t i = 0; i < points.size(); i++)
{
Point3d pt = points[i];
double x = pt.x;
double y = pt.y;
if (x < xMin)
xMin = x;
if (x > xMax)
xMax = x;
if (y < yMin)
yMin = y;
if (y > yMax)
yMax = y;
}
}
static void DrawFormattedText(HDC hdc, char text[], RECT rect)
{
// Draw the text with formatting options
DrawTextA(hdc, text, -1, &rect, DT_SINGLELINE | DT_NOCLIP);
}
//
// FUNCTION: WndProc(HWND, UINT, WPARAM, LPARAM)
//
// PURPOSE: Processes messages for the main window.
//
// WM_COMMAND - process the application menu
// WM_PAINT - Paint the main window
// WM_DESTROY - post a quit message and return
//
//
LRESULT CALLBACK WndProc(HWND hWnd, UINT message, WPARAM wParam, LPARAM lParam)
{
switch (message)
{
case WM_CREATE:
CreateButterflyGraphPoints();
break;
case WM_COMMAND:
{
int wmId = LOWORD(wParam);
// Parse the menu selections:
switch (wmId)
{
case IDM_ABOUT:
DialogBox(hInst, MAKEINTRESOURCE(IDD_ABOUTBOX), hWnd, About);
break;
case IDM_EXIT:
DestroyWindow(hWnd);
break;
default:
return DefWindowProc(hWnd, message, wParam, lParam);
}
}
break;
case WM_PAINT:
{
PAINTSTRUCT ps;
HDC hdc = BeginPaint(hWnd, &ps);
double h = 0, pi = 0, plm = 0, theta = 0;
double xMax = 0, xMin = 0, yMax = 0, yMin = 0;
FindMinMax(xMin, xMax, yMin, yMax);
float xSpan = (float)(xMax - xMin);
float ySpan = (float)(yMax - yMin);
RECT rect = { };
GetClientRect(hWnd, &rect);
float width = (float)(rect.right - rect.left + 1);
float height = (float)(rect.bottom - rect.top - 32 + 1);
float sx0 = 2.0f * width / 16.0f;
float sx1 = 14.0f * width / 16.0f;
float sy0 = 2.0f * height / 16.0f;
float sy1 = 14.0f * height / 16.0f;
float deltaX = xSpan / 8.0f;
float deltaY = ySpan / 8.0f;
float xSlope = (sx1 - sx0) / xSpan;
float xInter = (float)(sx0 - xSlope * xMin);
float ySlope = (sy0 - sy1) / ySpan;
float yInter = (float)(sy0 - ySlope * yMax);
float px = 0, py = 0, sx = 0, sy = 0;
POINT wPt = { };
int i = 0;
float x = (float)xMin;
float y = (float)yMax;
px = x;
py = y;
sx = xSlope * px + xInter;
sy = ySlope * py + yInter;
MoveToEx(hdc, (int)sx, (int)sy0, &wPt);
char buffer[128] = { };
while (i <= 8)
{
sx = xSlope * x + xInter;
wPt.x = wPt.y = 0;
MoveToEx(hdc, (int)sx, (int)sy0, &wPt);
LineTo(hdc, (int)sx, (int)sy1);
sprintf_s(buffer, "%5.4lf", x);
SIZE size = { };
GetTextExtentPoint32A(
hdc,
buffer,
(int)strlen(buffer),
&size);
RECT textRect = { };
textRect.left = (long)(sx - size.cx / 2.0f);
textRect.right = (long)(sx + size.cx / 2.0f);
textRect.top = (long)sy1;
textRect.bottom = (long)(sy1 + size.cy / 2.0f);
DrawFormattedText(hdc, buffer, textRect);
x += deltaX;
i++;
}
i = 0;
y = (float)yMin;
while (i <= 8)
{
sy = ySlope * y + yInter;
wPt.x = wPt.y = 0;
MoveToEx(hdc, (int)sx0, (int)sy, &wPt);
LineTo(hdc, (int)sx, (int)sy);
if (i != 0)
{
sprintf_s(buffer, "%+5.3lf", y);
SIZE size = { };
GetTextExtentPoint32A(
hdc,
buffer,
(int)strlen(buffer),
&size);
RECT textRect = { };
textRect.left = (long)(sx0 - size.cx - size.cx / 5.0f);
textRect.right = (long)(sx0 - size.cx / 2.0f);
textRect.top = (long)(sy - size.cy / 2.0f);
textRect.bottom = (long)(sy + size.cy / 2.0f);
DrawFormattedText(hdc, buffer, textRect);
}
y += deltaY;
i++;
}
HGDIOBJ bPenNew = NULL;
HGDIOBJ hPenOld = NULL;
bPenNew = CreatePen(PS_SOLID, 2, RGB(0, 0, 255));
hPenOld = SelectObject(hdc, bPenNew);
px = (float)points[0].x;
py = (float)points[0].y;
sx = xSlope * px + xInter;
sy = ySlope * py + yInter;
wPt.x = wPt.y = 0;
MoveToEx(hdc, (int)sx, (int)sy, &wPt);
for (size_t j = 1; j < points.size(); j++)
{
px = (float)points[j].x;
py = (float)points[j].y;
sx = xSlope * px + xInter;
sy = ySlope * py + yInter;
LineTo(hdc, (int)sx, (int)sy);
}
SelectObject(hdc, hPenOld);
DeleteObject(bPenNew);
EndPaint(hWnd, &ps);
}
break;
case WM_DESTROY:
PostQuitMessage(0);
break;
default:
return DefWindowProc(hWnd, message, wParam, lParam);
}
return 0;
}
// Message handler for about box.
INT_PTR CALLBACK About(HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
UNREFERENCED_PARAMETER(lParam);
switch (message)
{
case WM_INITDIALOG:
return (INT_PTR)TRUE;
case WM_COMMAND:
if (LOWORD(wParam) == IDOK || LOWORD(wParam) == IDCANCEL)
{
EndDialog(hDlg, LOWORD(wParam));
return (INT_PTR)TRUE;
}
break;
}
return (INT_PTR)FALSE;
}
Cello Horror © Wednesday, April 22, 2026, Rough cut MP3 using a Gibson SG Guitar and an E-Bow
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.
Numerical Explorations 