Numerical Explorations

// Integration1dTests.cpp
// Copyright Saturday, October 3, 2026
// James Pate Williams, Jr. BA, BS, MSwE, PhD
// One-dimensional integration tests
// Our canonical integrands are x^3 exp(x) for
// x = 0 to x = 1 and x^3 sin(x) for x = 0 to x = pi 

#include <iomanip>
#include <iostream>
using namespace std;

const int steps = 1024;

static double f1(double x)
{
	return x * x * x * exp(x);
}

static double f2(double x)
{
	return x * x * x * sin(x);
}

static double g1()
{
	return -2.0 * exp(1.0) + 6;
}

static double g2()
{
	double pi = 4.0 * atan(1.0);

	return -pi * pi * pi * cos(pi) + 6.0 * pi * cos(pi);
}

double TrapezoidalRule(
	int n, double a, double b, double(*fx)(double))
{
	double h = (b - a) / n;
	double s = 0.5 * (fx(a) + fx(b));
	double x = a + h;

	for (int i = 1; i < n; i++)
	{
		s += fx(x);
		x += h;
	}

	return h * s;
}

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 PrintResult(double x)
{
	cout << setw(15) << fixed;
	cout << right << setprecision(12);
	cout << setfill(' ');
	cout << x;
}

static void PrintError(double x, double exact)
{
	cout << '\t';
	cout << setw(15) << fixed;
	cout << right << setprecision(12);
	cout << setfill(' ');
	cout << 100.0 * fabs(x - exact) /
		fabs(exact) << endl;
}

int main()
{
	double gAct1 = g1();
	double gAct2 = g2();
	double pi = 4.0 * atan(1.0);
	double trap1 = TrapezoidalRule(steps, 0.0, 1.0, f1);
	double simp1 = SimpsonsRule(steps, 0.0, 1.0, f1);
	double trap2 = TrapezoidalRule(steps, 0.0,pi, f2);
	double simp2 = SimpsonsRule(steps, 0.0, pi, f2);

	cout << "Functions are: " << endl;
	cout << "f1(x) = x ^ 3 * exp(x)" << endl;
	cout << "f2(x) = x ^ 3 * sin(x)" << endl;
	cout << "Trapezoidal and Simpson Use 1024 Steps" << endl;
	cout << "         Value           % Error";
	cout << endl;
	cout << "trap1 = ";
	PrintResult(trap1);
	PrintError(trap1, gAct1);
	cout << "simp1 = ";
	PrintResult(simp1);
	PrintError(simp1, gAct1);
	cout << "true1 = ";
	PrintResult(gAct1);
	cout << endl;
	cout << "trap2 = ";
	PrintResult(trap2);
	PrintError(trap2, gAct2);
	cout << "simp2 = ";
	PrintResult(simp2);
	PrintError(simp2, gAct2);
	cout << "true2 = ";
	PrintResult(gAct2);
	cout << endl;
	return 0;
}

// Cardioid.cpp : Defines the entry point for the application.
// Copyright (c) Sunday, September 13, 2026
// by James Pate Williams, Jr. BA, BS, MSwE, PhD
// https://en.wikipedia.org/wiki/Cardioid

#include "framework.h"
#include "Cardioid.h"
#include <float.h>
#include <stdio.h>
#include <vector>

#define MAX_LOADSTRING 100

typedef struct tagPoint2d
{
    double x, y;
} Point2d, * PPoint2d;

// 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<Point2d> 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_CARDIOID, szWindowClass, MAX_LOADSTRING);
    MyRegisterClass(hInstance);

    // Perform application initialization:
    if (!InitInstance (hInstance, nCmdShow))
    {
        return FALSE;
    }

    HACCEL hAccelTable = LoadAccelerators(hInstance, MAKEINTRESOURCE(IDC_CARDIOID));

    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_CARDIOID));
    wcex.hCursor        = LoadCursor(nullptr, IDC_ARROW);
    wcex.hbrBackground  = (HBRUSH)(COLOR_WINDOW+1);
    wcex.lpszMenuName   = MAKEINTRESOURCEW(IDC_CARDIOID);
    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;
}

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++)
    {
        Point2d 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);
}

static void CreateCurve(double a, int n)
{
    double pi = 4.0 * atan(1.0), phi = 0.0;

    while (phi < 2.0 * pi)
    {
        double phi1 = 1.0 - cos(phi);
        double y = 2.0 * a * phi1 * cos(phi);
        double x = 2.0 * a * phi1 * sin(phi);
        Point2d pt = { 0, 0 };

        pt.x = x;
        pt.y = y;
        points.push_back(pt);
        phi += 0.001;
    }
}

//
//  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_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:
    {
        CreateCurve(1.0, 256);
        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;
        PAINTSTRUCT ps;
        POINT wPt = { };
        HDC hdc = BeginPaint(hWnd, &ps);
        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);

        return (INT_PTR)FALSE;
    }
    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;
}