NevillesMethodTest.php 3.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147
  1. <?php
  2. namespace MathPHP\Tests\NumericalAnalysis\Interpolation;
  3. use MathPHP\NumericalAnalysis\Interpolation\NevillesMethod;
  4. class NevillesMethodTest extends \PHPUnit\Framework\TestCase
  5. {
  6. /**
  7. * @test Interpolate
  8. * @dataProvider dataProviderForPolynomialAgrees
  9. * @param int $x
  10. * @param int $expected
  11. * @throws \Exception
  12. */
  13. public function testPolynomialAgrees(int $x, int $expected)
  14. {
  15. // Given
  16. $points = [[0, 0], [1, 5], [3, 2], [7, 10], [10, -4]];
  17. $roundoff = 0.0000001; // round off error
  18. // When
  19. $interpolated = NevillesMethod::interpolate($x, $points);
  20. // Then
  21. $this->assertEqualsWithDelta($expected, $interpolated, $roundoff);
  22. }
  23. /**
  24. * @return array (x, expected)
  25. */
  26. public function dataProviderForPolynomialAgrees(): array
  27. {
  28. return [
  29. [0, 0], // p(0) = 0
  30. [1, 5], // p(1) = 5
  31. [3, 2], // p(3) = 2
  32. [7, 10], // p(7) = 10
  33. [10, -4], // p(10) = -4
  34. ];
  35. }
  36. /**
  37. * @test Solve
  38. * @dataProvider dataProviderForSolve
  39. * @param float $x
  40. * @throws \Exception
  41. *
  42. * f(x) = x⁴ + 8x³ -13x² -92x + 96
  43. *
  44. * Given n points, the error in Nevilles Method (Lagrange polynomials) is proportional
  45. * to the max value of the nth derivative. Thus, if we if interpolate n at
  46. * 6 points, the 5th derivative of our original function f(x) = 0, and so
  47. * our resulting polynomial will have no error.
  48. *
  49. * p(x) agrees with f(x) at x = $_
  50. */
  51. public function testSolve($x)
  52. {
  53. // f(x) = x⁴ + 8x³ -13x² -92x + 96
  54. $f = function ($x) {
  55. return $x ** 4 + 8 * $x ** 3 - 13 * $x ** 2 - 92 * $x + 96;
  56. };
  57. // And
  58. $a = 0;
  59. $b = 10;
  60. $n = 5;
  61. $roundoff = 0.000001; // round off error
  62. // And
  63. $expected = $f($x);
  64. // When
  65. $actual = NevillesMethod::interpolate($x, $f, $a, $b, $n);
  66. // Then
  67. $this->assertEqualsWithDelta($expected, $actual, $roundoff);
  68. }
  69. /**
  70. * @test Solve
  71. * @dataProvider dataProviderForSolve
  72. * @param float $x
  73. * @throws \Exception
  74. *
  75. * f(x) = x⁴ + 8x³ -13x² -92x + 96
  76. *
  77. * The error is bounded by:
  78. * |f(x)-p(x)| = tol <= (max f⁽ⁿ⁺¹⁾(x))*(x-x₀)*...*(x-xn)/(n+1)!
  79. *
  80. * f'(x) = 4x³ +24x² -26x - 92
  81. * f''(x) = 12x² - 48x - 26
  82. * f'''(x) = 24x - 48
  83. * f⁽⁴⁾(x) = 24
  84. *
  85. * o, tol <= 24*(x-x₀)*...*(x-xn)/(4!) = (x-x₀)*...*(x-xn) where
  86. *
  87. * p(x) agrees with f(x) at x = $_
  88. */
  89. public function testSolveNonZeroError($x)
  90. {
  91. // f(x) = x⁴ + 8x³ -13x² -92x + 96
  92. $f = function ($x) {
  93. return $x ** 4 + 8 * $x ** 3 - 13 * $x ** 2 - 92 * $x + 96;
  94. };
  95. // And
  96. $a = 0;
  97. $b = 9;
  98. $n = 4;
  99. $x₀ = 0;
  100. $x₁ = 3;
  101. $x₂ = 6;
  102. $x₃ = 9;
  103. $tol = \abs(($x - $x₀) * ($x - $x₁) * ($x - $x₂) * ($x - $x₃));
  104. $roundoff = 0.0000001; // round off error
  105. // And
  106. $expected = $f($x);
  107. // When
  108. $actual = NevillesMethod::interpolate($x, $f, $a, $b, $n);
  109. // Then
  110. $this->assertEqualsWithDelta($expected, $actual, $tol + $roundoff);
  111. }
  112. /**
  113. * @return array p(x) agrees with f(x) at x = $_
  114. */
  115. public function dataProviderForSolve(): array
  116. {
  117. return [
  118. [0],
  119. [2],
  120. [4],
  121. [6],
  122. [8],
  123. [10],
  124. [-90],
  125. [-99],
  126. ];
  127. }
  128. }