#!/usr/local/bin/perl

# ----------------------------------------------------------------
# John Kerl
# kerl.john.r@gmail.com
# 2005-01-11
#
# This is an example of numerical differentiation in one variable:
# f(x) = sin(x)/x.  The estimated first and second derivatives are
# compared against analytically known results.
#
# Output:
#
# -6.28319 -0.00000 -0.15993 -0.15915 7.7346e-04 -0.05065 -0.05066 7.0666e-06
# -6.25177 -0.00502 -0.16139 -0.16068 7.0577e-04 -0.04637 -0.04638 6.8285e-06
# -6.22035 -0.01009 -0.16270 -0.16207 6.3680e-04 -0.04201 -0.04201 6.5836e-06
# -6.18894 -0.01521 -0.16389 -0.16332 5.6661e-04 -0.03757 -0.03757 6.3320e-06
# -6.15752 -0.02035 -0.16492 -0.16443 4.9524e-04 -0.03305 -0.03305 6.0738e-06
# -6.12611 -0.02554 -0.16582 -0.16539 4.2274e-04 -0.02846 -0.02846 5.8093e-06
# -6.09469 -0.03075 -0.16656 -0.16622 3.4917e-04 -0.02379 -0.02380 5.5387e-06
# -6.06327 -0.03598 -0.16716 -0.16689 2.7457e-04 -0.01907 -0.01907 5.2620e-06
# -6.03186 -0.04123 -0.16761 -0.16741 1.9899e-04 -0.01428 -0.01428 4.9796e-06
# -6.00044 -0.04650 -0.16791 -0.16779 1.2249e-04 -0.00942 -0.00943 4.6917e-06
#  ...
# ----------------------------------------------------------------

$pi  = 4.0 * atan2(1,1);
$xlo = -2.0 * $pi;
$xhi =  2.0 * $pi;

$nx  = 400;
$nx  = $ARGV[0] if (@ARGV);
$dx  = ($xhi - $xlo) / $nx;

# Alternatively:
#$dx  = 0.001;

for ($x = $xlo; $x < $xhi; $x += $dx) {
	$y = f($x);

	$estimated1 = (f($x + 0.5*$dx) - f($x - 0.5*$dx)) / $dx;
	$actual1    = fprime($x);

	$estimated2 = (f($x+$dx) - 2.0 * $y + f($x-$dx)) / $dx / $dx;
	$actual2    = fdoubleprime($x);

	printf "%10.5f %10.5f %10.5f %10.5f %9.4e %10.5f %10.5f %9.4e\n",
		$x, $y,
		$estimated1, $actual1, abs($estimated1 - $actual1),
		$estimated2, $actual2, abs($estimated2 - $actual2);
}

# ----------------------------------------------------------------
sub f {
	my ($x) = @_;
	if ($x == 0.0) {
		return 1.0;
	}
	else {
		return sin($x) / $x;
	}
}

# ----------------------------------------------------------------
sub fprime {
	my ($x) = @_;
	if ($x == 0.0) {
		return 0.0;
	}
	else {
		return ($x * cos($x) - sin($x)) / ($x * $x);
	}
}

# ----------------------------------------------------------------
sub fdoubleprime {
	my ($x) = @_;
	if (abs($x) < 1.0e-12) {
		return -1.0 / 3.0;
	}
	else {
		return (2 * sin($x) - 2 * $x * cos($x) - $x*$x*sin($x))
			/ ($x * $x * $x);
	}
}