Category: Computer Science

Comparison of Euler Method, Runge-Kutta 2, Runge-Kutta 4, and Taylor Series Solutions of a Single Ordinary Differential Equation

#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;
	}
}

Numerical Differentiation in C# by James Pate Williams, Jr. Translated from FORTRAN 77 Formulas

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);
        }
    }
}

Adams-Moulton Predictor Corrector Method to Solve a One-Dimensional Ordinary Differential Equation Translated from FORTRAN 77 Code by James Pate Williams, Jr.

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);
            }
        }
    }
}

Determinant Calculations: Bareiss Algorithm and Gaussian Elimination by James Pate Williams, Jr.

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];
        }
    }
}

Microsoft C# Constraint Satisfaction Dynamic Link Library (DLL) and Test Application

Constraint Satisfaction DLL Source Code Files
Source Code File 	Lines of Code
Arc.cs	42
Common.cs	214
Label.cs	62
NQPArcConsisitencies.cs	148
NQPBacktrack.cs	51
NQPSosicGu.cs	234
Vertex.cs	18
YakooGCPWCSA.cs	160
Total	929

Source Code File	Lines of Code
GCPGraphForm.cs	173
MainForm.cs	42
NQPGraphForm.cs	137
NQPTableForm.cs	223
Total	575

Source
1.	Foundations of Constraint Satisfaction by Edward Tsang

Sosic and Gu Fast N-Queens Problem Solver Implemented by James Pate Williams, Jr.

I implemented the Sosic and Gu fast N-Queens solver first in C++ and later in C#. The five screen shots below were created by latest C# application. None of my other algorithms can compete with this method.

More N-Queens Puzzle Results by James Pate Williams, Jr.

N-Queens Tests				
8 Queens	Microseconds		50 Samples	
Algorithm	Minimum	Average	Maximum	Std Dev
Arc-Consistency 3	1846	4148	14220	2294
Arc-Consistency 4	3850	8987	27903	5022
Backjump	175	1057	3873	883
Backmark	198	737	3323	596
Backtracking	171	789	3132	652
Hill-Climber	1991	2097	2913	157
				
9 Queens	Microseconds		50 Samples	
Algorithm	Minimum	Average	Maximum	Std Dev
Arc-Consistency 3	2628	5785	15748	3202
Arc-Consistency 4	7492	18057	56587	11114
Backjump	208	1410	6184	1325
Backmark	214	866	4468	880
Backtracking	210	1269	6156	1186
Hill-Climber	1309	1398	1953	102
				
10 Queens	Microseconds		50 Samples	
Algorithm	Minimum	Average	Maximum	Std Dev
Arc-Consistency 3	3820	10680	30708	6963
Arc-Consistency 4	12624	46645	148221	31212
Backjump	248	3970	16736	3805
Backmark	283	1504	7613	1433
Backtracking	227	2857	12882	3022
Hill-Climber	24575	25608	33563	2037

My hill-climber results are pretty disappointing. I have gotten some decent runs of back-marking for n = 50 and n = 60 queens. I show the n = 60 queens experiment below.

==Menu==
**Instance Submenu**
1 Single Instance Test
2 Multiple Instances Test
3 Exit
Enter an option '1' to '3': 2
**Algorithm Submenu**
1 Arc-Consistency 3
2 Arc-Consistency 4
3 Backjump
4 Backmark
5 Chronological Bactracking
6 Evolutionary Hill-Climber
7 Exit
Enter an algorithm ('1' to '7'): 4
Enter the number of queens for option '2': 60
Enter the number of samples to be analyzed: 10
Minimum Runtime in Microseconds: 4819
Average Runtime in Microseconds: 10419330
Maximum Runtime in Micorseconds: 45439444
Standard Sample Deviation: 16665688

Notice the misspelling of microseconds.

Six C++ Algorithms to Solve the N-Queens Puzzle (Problem) by James Pate Williams, Jr.

Most of the six algorithms were developed back in the years 1999 to 2001 when I was a Master of Software Engineering and Doctoral candidate at Auburn University. Professor Gerry V. Dozier was my Master research advisor. The algorithms developed were in alphabetical order:

  • Arc-Consistency 3
  • Arc-Consistency 4
  • Backjump
  • Backmark
  • Chronological Backtracking
  • Evolutionary Hill-Climber (My algorithm)

The first five algorithms were stated in pseudocode in the excellent treatise Foundations of Constraint Satisfaction by Edward Tsang. Below are six runs of my test C++ application (Win32 Console App) using 8 queens and 50 samples.