Numerical Explorations

// FiniteDifference.cpp : This file contains the 'main' function.
// Program execution begins and ends there.
// Copyright (c) Wednesday, August 9, 2026
// by James Pate Williams, Jr., BA, BS, MSwE, PhD
// Reference: "Elementary Numerical Analysis an Algorithmic Approach
// Third Edition" © 1980 by S. D. Conte and Carl de Boor 

#include <math.h>
#include <stdio.h>

#define MAX_ROWS 8192

double ac[MAX_ROWS];
double fv[MAX_ROWS];
double gv[MAX_ROWS];
double qv[MAX_ROWS];
double xv[MAX_ROWS];
double yv[MAX_ROWS];

double ta[MAX_ROWS];
double tb[MAX_ROWS];
double td[MAX_ROWS];
double tc[MAX_ROWS];
double tx[MAX_ROWS];

int SolveTridiagonal(int n)
{
    for (int k = 2; k <= n; k++)
    {
        if (td[k - 1] == 0)
            return 0;

        double m = ta[k] / td[k - 1];
        td[k] -= m * tc[k - 1];
        tb[k] -= m * tb[k - 1];
    }

    if (td[n] == 0)
        return 0;

    tx[n] = tb[n] / td[n];

    for (int k = n - 1; k >= 1; k--)
        tx[k] = (tb[k] - tc[k] * tx[k + 1]) / td[k];

    return 1;
}

double f1(double x)
{
	return 0.0;
}

double g1(double x)
{
	return -1.0;
}

double q1(double x)
{
	return 0.0;
}

double f2(double x)
{
	return 0.0;
}

double g2(double x)
{
	return 1.0;
}

double q2(double x)
{
	return 0.0;
}

double f3(double x)
{
	return x;
}

double g3(double x)
{
	return 1.0;
}

double q3(double x)
{
	return 2.0 * x;
}

double f4(double x)
{
	return 2.0;
}

double g4(double x)
{
	return 1.0;
}

double q4(double x)
{
	return x;
}

double f5(double x)
{
	return 1.0 / x;
}

double g5(double x)
{
	return 1.0;
}

double q5(double x)
{
	return 0.0;
}

double FD_Solution(
	double a,
	double b,
	double ya,
	double yb,
	double (*f)(double),
	double (*g)(double),
	double (*q)(double),
	int N)
{
	double h = (b - a) / N;

	for (int i = 1; i <= N - 1; i++)
	{
		xv[i] = a + i * h;
		fv[i] = f(xv[i]);
		gv[i] = g(xv[i]);
		qv[i] = q(xv[i]);
	}

	tb[1] = h * h * qv[1] - (1.0 - 0.5 * h * fv[1]) * ya;

	for (int i = 2; i <= N - 1; i++)
		tb[i] = h * h * qv[i];

	for (int i = 1; i <= N - 1; i++)
		td[i] = -2.0 + h * h * gv[i];

	for (int i = 1; i <= N - 2; i++)
		tc[i] = 1.0 + 0.5 * h * fv[i];

	for (int i = 2; i <= N - 1; i++)
		ta[i] = 1.0 - 0.5 * h * fv[i];

	tb[N - 1] = h * h * qv[N - 1] - (1.0 + 0.5 * h * fv[N - 1]) * yb;

	return SolveTridiagonal(N - 1);
}

double IS_Solution(
	double a,
	double b,
	double ya,
	double yb,
	double x)
{
	double A[13] = { 0 };

	A[0] = 1.0;
	A[1] = 0.0;
	A[2] = -0.25;
	A[3] = 0.0;
	A[4] = A[0] / 64.0;
	A[5] = 0.0;
	A[6] = -A[0] / (36.0 * 64.0);
	A[7] = 0.0;
	A[8] = A[0] / (64.0 * 36.0 * 64.0);
	A[9] = 0.0;
	A[10] = -A[8] / 100.0;
	A[11] = 0.0;
	A[12] = -A[10] / 144.0;

	double s = A[12];

	for (int i = 11; i >= 0; i--)
		s = s * x + A[i];

	return s;
}

int main()
{
	double h = 0.05, x0 = 0.0, x1 = 1.0, y0 = 0, y1 = 1.0;
	int N = (int)((x1 - x0) / h);
	double fd1 = FD_Solution(
		x0, x1, y0, y1, f1, g1, q1, N);
	FILE* file = 0;
	int errno = fopen_s(&file, "Chapter9.txt", "w");

	if (errno != 0)
		return -1;

	fprintf_s(file, "Example 9.1\r\n");
	fprintf_s(file, " x\tapproximate\texact\t\tpercent error\r\n");

	for (int i = 1; i < N; i++)
	{
		double exact = sinh(xv[i]) / sinh(1.0);
		double error = 100.0 * (fabs(tx[i] - exact) / fabs(exact));
		fprintf_s(file, "%3.2lf\t%11.10lf\t%11.10lf\t%11.10lf\r\n",
			xv[i], tx[i], exact, error);

	}

	h = 0.25;
	x0 = 0.0, x1 = 1.0, y0 = 0, y1 = 1.0;
	N = (int)((x1 - x0) / h);
	double fd2 = FD_Solution(
		x0, x1, y0, y1, f2, g2, q2, N);

	fprintf_s(file, "Exercise 9.1-1\r\n");
	fprintf_s(file, " x\tapprox\r\n");

	for (int i = 1; i < N; i++)
		fprintf_s(file, "%3.2lf\t%5.4lf\r\n", xv[i], tx[i]);

	h = 0.1;
	x0 = 0.0, x1 = 1.0, y0 = 1.0, y1 = 0.0;
	N = (int)((x1 - x0) / h);
	double fd3 = FD_Solution(
		x0, x1, y0, y1, f3, g3, q3, N);

	fprintf_s(file, "Exercise 9.1-3\r\n");
	fprintf_s(file, " x\t approximate\r\n");

	for (int i = 1; i < N; i++)
		fprintf_s(file, "%3.2lf\t%11.10lf\r\n", xv[i], tx[i]);

	h = 1.0 / 16.0;
	x0 = 0.0, x1 = 1.0, y0 = 0.0, y1 = 1.0;
	N = (int)((x1 - x0) / h);
	double fd4 = FD_Solution(
		x0, x1, y0, y1, f4, g4, q4, N);

	fprintf_s(file, "Exercise 9.1-4\r\n");
	fprintf_s(file, " x\t approximate\r\n");

	for (int i = 1; i < N; i++)
		fprintf_s(file, "%3.2lf\t%11.10lf\r\n", xv[i], tx[i]);

	h = 1.0 / 16.0;
	x0 = 0.0, x1 = 2.0, y0 = 1.0, y1 = 2.238907555e-01;
	N = (int)((x1 - x0) / h);
	double fd5 = FD_Solution(
		x0, x1, y0, y1, f5, g5, q5, N);

	fprintf_s(file, "Bessel Equation of Order Zero\r\n");
	fprintf_s(file, "Approximate Finite Difference Solution\r\n");

	fprintf_s(file, " x\t approximate\r\n");

	for (int i = 1; i < N; i++)
		fprintf_s(file, "%3.2lf\t%11.10lf\r\n", xv[i], tx[i]);

	fprintf_s(file, "Bessel Equation of Order Zero\r\n");
	fprintf_s(file, "Approximate Infinite Series Solution\r\n");

	fprintf_s(file, " x\t approximate\r\n");

	for (int i = 1; i < N; i++)
	{
		double x = x0 + i * h;
		fprintf_s(file, "%3.2lf\t%11.10lf\r\n", x,
			IS_Solution(x0, x1, y0, y1, x));
	}

    return 0;
}