Category: Numerical Analysis

Yet Another Revisitation of Reproducing Ordnance Pamphlet 770 by James Pate Williams, Jr.

This is another attempt to reproduce the United States Navy’s Ordnance Pamphlet 770: https://eugeneleeslover.com/USN-GUNS-AND-RANGE-TABLES/OP-770-1.html which contains ballistic tables for the battleship USS Iowa (BB-61) artillery (16-inch/50 caliber) and is dated October 1941. My C# Windows desktop application is capable of calculating the elevation from range table which has the columns range in yards, angle of elevation in degrees and minutes, positive angle of fall in degrees and minutes, time of flight in seconds, apogee also called summit in feet, striking velocity in feet per second, and energy in foot pound force. Three corrections can be applied to the trajectory: trunnion height in feet, acceleration of gravity correction, and the curvature of the Earth correction (Vincenty calculation). The first image below is the ballistic settings interface. The second image is the uncorrected table. The third image is the application of a trunnion height of 32 feet. The fourth image is the curvature of the Earth correction. The fifth image is the trunnion height of 32 feet and Vincenty corrections. It is to be noticed that the striking velocity and kinetic energy are the only non-monotonically increasing or decreasing data fields.

Parabolic Cylinder Functions by James Pate Williams, Jr.

Solution to the differential equation: d2y/dx2 = -0.25 * x * x * y – a * y which is valid for all real and complex numbers. We examine the real solutions to the second order ordinary differential equations. See “Handbook of Mathematical Functions” by Milton Abramowitz and Irene A. Stegun, Chapter 19. Parabolic Cylinder functions. Page 686 Equation 19.2.7 Recurrence formula and 19.2.5 Series solutions.

Win32 C Orthogonal Polynomials

I learned about the Laguerre and Legendre polynomials when I first read “Introduction to Quantum Mechanics” by Pauling and Wilson way back in the early 1970s. I later learned about the Chebyshev and other orthogonal polynomials. Beginning on March 30, 2015, I created yet another application to graph various orthogonal polynomials in C#. A few days ago, I wrote a Win32 C application to graph the Chebyshev, Laguerre, and Legendre orthogonal polynomials. At a later date, I will probably add Hermite and Jacobi polynomials.


Win32 C Romberg Extrapolation by James Pate Williams, Jr.

// Romberg.cpp : Defines the entry point for the application.
// https://learn.microsoft.com/en-us/windows/win32/controls/create-a-simple-combo-box

#include "stdafx.h"
#include <CommCtrl.h>
#include <math.h>
#include <stdio.h>
#include "Romberg.h"

#define MAX_LOADSTRING 100

// Global Variables:
HINSTANCE hInst;                                // current instance
WCHAR szTitle[MAX_LOADSTRING];                  // The title bar text
WCHAR szWindowClass[MAX_LOADSTRING];            // the main window class name

// 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_PTR CALLBACK    RombergDialog(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_ROMBERG, szWindowClass, MAX_LOADSTRING);
    MyRegisterClass(hInstance);

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

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

    MSG msg;

    // Main message loop:
    while (GetMessage(&msg, NULL, 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;

    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_ROMBERG));
    wcex.hCursor        = LoadCursor(NULL, IDC_ARROW);
    wcex.hbrBackground  = (HBRUSH)(COLOR_WINDOW+1);
    wcex.lpszMenuName   = MAKEINTRESOURCEW(IDC_ROMBERG);
    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, NULL, NULL, hInstance, NULL);

   if (!hWnd)
   {
      return FALSE;
   }

   ShowWindow(hWnd, nCmdShow);
   UpdateWindow(hWnd);

   return TRUE;
}

//
//  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 ID_OPTION_ROMBERGINTEGRATION:
				DialogBox(hInst, MAKEINTRESOURCE(IDD_DIALOG1), hWnd, RombergDialog);
				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);
            // TODO: Add any drawing code that uses hdc here...
            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;
}

double f0(double x)
{
	return exp(-x * x);
}

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

double f2(double x)
{
	double pi = 4.0 * atan(1.0);
	return sin(101.0 * pi * x);
}

double f3(double x)
{
	double pi = 4.0 * atan(1.0);
	return 1.0 + sin(10.0 * pi * x);
}

double f4(double x)
{
	return fabs(x - 1.0 / 3.0);
}

double f5(double x)
{
	return sqrt(x);
}

double f6(double x)
{
	double result = 1.0;

	if (x != 0.0)
		result = sin(x) / x;

	return result;
}

/* https://en.wikipedia.org/wiki/Romberg%27s_method#Implementation */

char table[8192];
double R1[2048], R2[2048];

void print_row(size_t i, double* R)
{
	char buffer[128];
	sprintf_s(buffer, 128, "R[%2zu] = ", i);
	strcat_s(table, 8192, buffer);

	for (size_t j = 0; j <= i; ++j)
	{
		sprintf_s(buffer, 128, "%.*lf ", 16, R[j]);
		strcat_s(table, 8192, buffer);
	}

	strcat_s(table, 8192, "\r\n");
}

double RombergExtrapolation(
	double(*f)(double),
	double a, double b,
	size_t max_steps,
	double acc)
{
	double *Rp = &R1[0], *Rc = &R2[0];
	double h = (b - a);
	Rp[0] = (f(a) + f(b))*h*.5;

	print_row(0, Rp);

	for (size_t i = 1; i < max_steps; ++i)
	{
		h /= 2.;
		double c = 0;
		size_t ep = 1 << (i - 1);
		
		for (size_t j = 1; j <= ep; ++j)
		{
			c += f(a + (2 * j - 1)*h);
		}

		Rc[0] = h*c + .5*Rp[0];

		for (size_t j = 1; j <= i; ++j)
		{
			double n_k = pow(4, j);
			Rc[j] = (n_k*Rc[j - 1] - Rp[j - 1]) / (n_k - 1);
		}

		print_row(i, Rc);

		if (i > 1 && fabs(Rp[i - 1] - Rc[i]) < acc)
		{
			return Rc[i];
		}

		double *rt = Rp;
		Rp = Rc;
		Rc = rt;
	}

	return Rp[max_steps - 1];
}

INT_PTR CALLBACK RombergDialog(HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
	UNREFERENCED_PARAMETER(lParam);
	WORD word = LOWORD(wParam);
	static TCHAR functions[7][128] =
	{
		TEXT("f(x) = exp(-x * x)"),
		TEXT("f(x) = x * x"),
		TEXT("f(x) = sin(101 * pi * x)"),
		TEXT("f(x) = 1 + sin(10 * pi * x)"),
		TEXT("f(x) = abs(x - 1.0 / 3.0)"),
		TEXT("f(x) = sqrt(x)"),
		TEXT("f(x) = sin(x) / x or 1 if x = 0")
	};
	static TCHAR A[1024];
	static HWND cbWindow;
	static int k = 0;
	static int itemIndex;

	switch (message)
	{
	case WM_INITDIALOG:
		cbWindow = GetDlgItem(hDlg, IDC_COMBO1);
		memset(&A, 0, sizeof(A));
		for (k = 0; k < 6; k++)
		{
			// Add string to combobox.
			SendMessage(cbWindow, CB_ADDSTRING, (WPARAM)k, (LPARAM)functions[k]);
		}
		// Send the CB_SETCURSEL message to display an initial item 
		// in the selection field  
		SendMessage(cbWindow, CB_SETCURSEL, (WPARAM)0, (LPARAM)0);
		return (INT_PTR)TRUE;
	case WM_COMMAND:
		if (word == IDOK || word == IDCANCEL)
		{
			EndDialog(hDlg, word);
			return (INT_PTR)TRUE;
		}
		if (HIWORD(wParam) == CBN_SELCHANGE)
			// If the user makes a selection from the list:
			//   Send CB_GETCURSEL message to get the index of the selected list item.
			//   Send CB_GETLBTEXT message to get the item.
			//   Display the item in a messagebox.
		{
			itemIndex = SendMessage((HWND)lParam, (UINT)CB_GETCURSEL,
				(WPARAM)0, (LPARAM)0);
			TCHAR listItem[256];
			(TCHAR)SendMessage((HWND)lParam, (UINT)CB_GETLBTEXT,
				(WPARAM)itemIndex, (LPARAM)listItem);
		}
		else if (word == IDC_BUTTON1)
		{
			char buffer[128];
			GetDlgItemTextA(hDlg, IDC_EDIT1, buffer, 128);
			double a = atof(buffer);
			GetDlgItemTextA(hDlg, IDC_EDIT2, buffer, 128);
			double b = atof(buffer);
			double integral = 0.0;
			GetDlgItemTextA(hDlg, IDC_EDIT4, buffer, 128);
			double acc = atof(buffer);
			int max_steps = GetDlgItemInt(hDlg, IDC_EDIT5, FALSE, FALSE);

			switch (itemIndex)
			{
			case 0:
				integral = RombergExtrapolation(f0, a, b, max_steps, acc);
				break;
			case 1:
				integral = RombergExtrapolation(f1, a, b, max_steps, acc);
				break;
			case 2:
				integral = RombergExtrapolation(f2, a, b, max_steps, acc);
				break;
			case 3:
				integral = RombergExtrapolation(f3, a, b, max_steps, acc);
				break;
			case 4:
				integral = RombergExtrapolation(f4, a, b, max_steps, acc);
				break;
			case 5:
				integral = RombergExtrapolation(f5, a, b, max_steps, acc);
				break;
			case 6:
				integral = RombergExtrapolation(f6, a, b, max_steps, acc);
				break;
			}

			sprintf_s(buffer, 128, "Integral = %.*lf\r\n", 16, integral);
			strcat_s(table, 8192, buffer);
			SetDlgItemTextA(hDlg, IDC_EDIT7, table);
		}

		else if (word == IDC_BUTTON3)
		{
			table[0] = '\0';
			SetDlgItemTextA(hDlg, IDC_EDIT7, table);
		}
	}
	return (INT_PTR)FALSE;
}

Win32 C Real Double Precision and Integer Matrix Add and Subtract by James Pate Williams, Jr.

void DblAdd(int l, double* a, double* b, double* c)
{
	_asm
	{
		mov eax, a; load a matrix address
		mov ebx, b; load a matrix address
		mov edi, c; load result address
		mov edx, l; load matrix length
		finit;
	ll:	fld qword ptr[eax];
		fld qword ptr[ebx];
		fadd;
		fstp qword ptr[edi];
		fwait;
		add edi, 8;
		add eax, 8;
		add ebx, 8;
		dec edx;
		jnz ll;
	}
}

void DblSub(int l, double* a, double* b, double* c)
{
	_asm
	{
		mov eax, a; load a matrix address
		mov ebx, b; load a matrix address
		mov edi, c; load result address
		mov edx, l; load matrix length
		finit;
	ll:	fld qword ptr[eax];
		fld qword ptr[ebx];
		fsub;
		fstp qword ptr[edi];
		fwait;
		add edi, 8;
		add eax, 8;
		add ebx, 8;
		dec edx;
		jnz ll;
	}
}

Eigenfunctions of a Nonlinear Second Order Ordinary Differential Equation Initial Value Eigenvalue Problem by James Pate Williams, Jr.

Solutions of a nonlinear second order ordinary differential equation initial value eigenvalue problem:

y”(x) = x * x * y(x) + n * n * y(x) * y(x) for all x in [0, 1]
y(0) = 0
y'(0) = 1
n in [0, 1, 2, …]

I graphed the first five eigenfunctions.

Graphs for n = 0 to n = 4