NaturalCubicSplineTest.php 3.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145
  1. <?php
  2. namespace MathPHP\Tests\NumericalAnalysis\Interpolation;
  3. use MathPHP\NumericalAnalysis\Interpolation\NaturalCubicSpline;
  4. use MathPHP\Expression\Polynomial;
  5. class NaturalCubicSplineTest extends \PHPUnit\Framework\TestCase
  6. {
  7. /**
  8. * @test Interpolated piecewise function computes expected values: p(x) = expected
  9. * @dataProvider dataProviderForPiecewiseForPolynomialAgrees
  10. * @param int $x
  11. * @param int $expected
  12. * @throws \Exception
  13. */
  14. public function testPolynomialAgrees(int $x, int $expected)
  15. {
  16. // Given
  17. $points = [[0, 0], [1, 5], [3, 2], [7, 10], [10, -4]];
  18. // And
  19. $p = NaturalCubicSpline::interpolate($points);
  20. // When
  21. $evaluated = $p($x);
  22. // Then
  23. $this->assertEqualsWithDelta($expected, $evaluated, 0.00001);
  24. }
  25. /**
  26. * @return array (x, expected)
  27. */
  28. public function dataProviderForPiecewiseForPolynomialAgrees(): array
  29. {
  30. return [
  31. [0, 0], // p(0) = 0
  32. [1, 5], // p(1) = 5
  33. [3, 2], // p(3) = 2
  34. [7, 10], // p(7) = 10
  35. [10, -4], // p(10) = -4
  36. ];
  37. }
  38. /**
  39. * @test Solve zero error
  40. * @dataProvider dataProviderForSolve
  41. * @param float $x
  42. * @throws \Exception
  43. *
  44. * f(x) = 8x³ -13x² -92x + 96
  45. *
  46. * The error in the Cubic Spline Interpolating Polynomial is proportional
  47. * to the max value of the 4th derivative. Thus, if our input Function
  48. * is a 3rd-degree polynomial, the fourth derivative will be zero, and
  49. * thus we will have zero error.
  50. *
  51. * p(x) agrees with f(x) at x = $_
  52. */
  53. public function testSolveZeroError($x)
  54. {
  55. // Given f(x) = 8x³ -13x² -92x + 96
  56. $f = new Polynomial([8, -13, -92, 96]);
  57. // And
  58. $a = 0;
  59. $b = 10;
  60. $n = 50;
  61. $tol = 0;
  62. $roundoff = 0.0001; // round off error
  63. // And
  64. $p = NaturalCubicSpline::interpolate($f, $a, $b, $n);
  65. $expected = $f($x);
  66. // When
  67. $evaluated = $p($x);
  68. // Then
  69. $this->assertEqualsWithDelta($expected, $evaluated, $tol + $roundoff);
  70. }
  71. /**
  72. * @test Solve non-zero error
  73. * @dataProvider dataProviderForSolve
  74. * @param float $x
  75. * @throws \Exception
  76. *
  77. * f(x) = x⁴ + 8x³ -13x² -92x + 96
  78. *
  79. * The error is bounded by:
  80. * |f(x)-p(x)| = tol <= (1/4!) * h⁴ * max f⁽⁴⁾(x)
  81. * where h = max hᵢ
  82. *
  83. * f'(x) = 4x³ +24x² -26x - 92
  84. * f''(x) = 12x² - 48x - 26
  85. * f'''(x) = 24x - 48
  86. * f⁽⁴⁾(x) = 24
  87. *
  88. * So, tol <= (1/24) * (1/5)⁴ * 24 = (1/5)⁴
  89. *
  90. * p(x) agrees with f(x) at x = $_
  91. */
  92. public function testSolveNonZeroError($x)
  93. {
  94. // Given f(x) = x⁴ + 8x³ -13x² -92x + 96
  95. $f = new Polynomial([1, 8, -13, -92, 96]);
  96. // And
  97. $a = 0;
  98. $b = 10;
  99. $n = 51;
  100. // And
  101. $tol = $tol = 0.2 ** 4; // So, tol <= (1/24) * (1/5)⁴ * 24 = (1/5)⁴
  102. $roundoff = 0.000001; // round off error
  103. // And
  104. $p = NaturalCubicSpline::interpolate($f, $a, $b, $n);
  105. $expected = $f($x);
  106. // When
  107. $evaluated = $p($x);
  108. // Then
  109. $this->assertEqualsWithDelta($expected, $evaluated, $tol + $roundoff);
  110. }
  111. /**
  112. * @return array p(x) agrees with f(x) at x = $_
  113. */
  114. public function dataProviderForSolve(): array
  115. {
  116. return [
  117. [0],
  118. [2],
  119. [4],
  120. [6],
  121. [8],
  122. [10],
  123. [7.32], // not a node
  124. ];
  125. }
  126. }