TrapezoidalRuleTest.php 7.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249
  1. <?php
  2. namespace MathPHP\Tests\NumericalAnalysis\NumericalIntegration;
  3. use MathPHP\Expression\Polynomial;
  4. use MathPHP\NumericalAnalysis\NumericalIntegration\TrapezoidalRule;
  5. class TrapezoidalRuleTest extends \PHPUnit\Framework\TestCase
  6. {
  7. /**
  8. * @test approximate with endpoints (0, 1) and (3, 16)
  9. * @throws \Exception
  10. *
  11. * f(x) = x² + 2x + 1
  12. * Antiderivative F(x) = (1/3)x³ + x² + x
  13. * Indefinite integral over [0, 3] = F(3) - F(0) = 21
  14. *
  15. * h₁, h₂, ... denotes the size on interval 1, 2, ...
  16. * ζ₁, ζ₂, ... denotes the max of the second derivative of f(x) on
  17. * interval 1, 2, ...
  18. * f'(x) = 2x + 2
  19. * f''(x) = 2
  20. * ζ = f''(x) = 2
  21. * Error = Sum(ζ₁h₁³ + ζ₂h₂³ + ...) = 2 * Sum(h₁³ + h₂³ + ...)
  22. *
  23. * Approximate with endpoints: (0, 1) and (3, 16)
  24. * Error = 2 * ((3 - 0)²) = 18
  25. */
  26. public function testApproximateWithEndpoints()
  27. {
  28. // Given
  29. $points = [[0, 1], [3, 16]];
  30. $tol = 18;
  31. $expected = 21;
  32. // When
  33. $x = TrapezoidalRule::approximate($points);
  34. // Then
  35. $this->assertEqualsWithDelta($expected, $x, $tol);
  36. }
  37. /**
  38. * @test approximate with endpoints and one interior point: (0, 1), (1, 4) and (3, 16)
  39. * @throws \Exception
  40. *
  41. * f(x) = x² + 2x + 1
  42. * Antiderivative F(x) = (1/3)x³ + x² + x
  43. * Indefinite integral over [0, 3] = F(3) - F(0) = 21
  44. *
  45. * h₁, h₂, ... denotes the size on interval 1, 2, ...
  46. * ζ₁, ζ₂, ... denotes the max of the second derivative of f(x) on
  47. * interval 1, 2, ...
  48. * f'(x) = 2x + 2
  49. * f''(x) = 2
  50. * ζ = f''(x) = 2
  51. * Error = Sum(ζ₁h₁³ + ζ₂h₂³ + ...) = 2 * Sum(h₁³ + h₂³ + ...)
  52. *
  53. * Approximate with endpoints and interior point: (0, 1), (1, 4) and (3, 16)
  54. * Error = 2 * ((1 - 0)² + (3 - 1)²) = 10
  55. */
  56. public function testApproximateWithEndpointsAndOneInteriorPoint()
  57. {
  58. // Given
  59. $points = [[0, 1], [1, 4], [3, 16]];
  60. $tol = 10;
  61. $expected = 21;
  62. // When
  63. $x = TrapezoidalRule::approximate($points);
  64. // Then
  65. $this->assertEqualsWithDelta($expected, $x, $tol);
  66. }
  67. /**
  68. * @test approximate with endpoints and two interior points: (0, 1), (1, 4), (2, 9) and (3, 16)
  69. * @throws \Exception
  70. *
  71. * f(x) = x² + 2x + 1
  72. * Antiderivative F(x) = (1/3)x³ + x² + x
  73. * Indefinite integral over [0, 3] = F(3) - F(0) = 21
  74. *
  75. * h₁, h₂, ... denotes the size on interval 1, 2, ...
  76. * ζ₁, ζ₂, ... denotes the max of the second derivative of f(x) on
  77. * interval 1, 2, ...
  78. * f'(x) = 2x + 2
  79. * f''(x) = 2
  80. * ζ = f''(x) = 2
  81. * Error = Sum(ζ₁h₁³ + ζ₂h₂³ + ...) = 2 * Sum(h₁³ + h₂³ + ...)
  82. *
  83. * Approximate with endpoints and interior point: (0, 1), (1, 4), (2, 9) and (3, 16)
  84. * Error = 2 * ((1 - 0)² + (2 - 1)² + (3 - 2)²) = 6
  85. */
  86. public function testApproximateWithEndpointsAndTwoInteriorPoints()
  87. {
  88. // Given
  89. $points = [[0, 1], [1, 4], [2, 9], [3, 16]];
  90. $tol = 6;
  91. $expected = 21;
  92. // When
  93. $x = TrapezoidalRule::approximate($points);
  94. // Then
  95. $this->assertEqualsWithDelta($expected, $x, $tol);
  96. }
  97. /**
  98. * @test approximate with endpoints and two interior points not sorted: (0, 1), (1, 4), (2, 9) and (3, 16)
  99. * @throws \Exception
  100. *
  101. * f(x) = x² + 2x + 1
  102. * Antiderivative F(x) = (1/3)x³ + x² + x
  103. * Indefinite integral over [0, 3] = F(3) - F(0) = 21
  104. *
  105. * h₁, h₂, ... denotes the size on interval 1, 2, ...
  106. * ζ₁, ζ₂, ... denotes the max of the second derivative of f(x) on
  107. * interval 1, 2, ...
  108. * f'(x) = 2x + 2
  109. * f''(x) = 2
  110. * ζ = f''(x) = 2
  111. * Error = Sum(ζ₁h₁³ + ζ₂h₂³ + ...) = 2 * Sum(h₁³ + h₂³ + ...)
  112. *
  113. * Approximate with endpoints and interior point: (0, 1), (1, 4), (2, 9) and (3, 16)
  114. * Error = 2 * ((1 - 0)² + (2 - 1)² + (3 - 2)²) = 6
  115. */
  116. public function testApproximateWithEndpointsAndTwoInteriorPointsNotSorted()
  117. {
  118. // Given
  119. $points = [[1, 4], [3, 16], [0, 1], [2, 9]];
  120. $tol = 6;
  121. $expected = 21;
  122. // When
  123. $x = TrapezoidalRule::approximate($points);
  124. // Then
  125. $this->assertEqualsWithDelta($expected, $x, $tol);
  126. }
  127. /**
  128. * @test approximate using callback
  129. * @throws \Exception
  130. *
  131. * f(x) = x² + 2x + 1
  132. * Antiderivative F(x) = (1/3)x³ + x² + x
  133. * Indefinite integral over [0, 3] = F(3) - F(0) = 21
  134. *
  135. * h₁, h₂, ... denotes the size on interval 1, 2, ...
  136. * ζ₁, ζ₂, ... denotes the max of the second derivative of f(x) on
  137. * interval 1, 2, ...
  138. * f'(x) = 2x + 2
  139. * f''(x) = 2
  140. * ζ = f''(x) = 2
  141. * Error = Sum(ζ₁h₁³ + ζ₂h₂³ + ...) = 2 * Sum(h₁³ + h₂³ + ...)
  142. */
  143. public function testApproximateUsingCallback()
  144. {
  145. // Given x² + 2x + 1
  146. $func = $func = function ($x) {
  147. return $x ** 2 + 2 * $x + 1;
  148. };
  149. $start = 0;
  150. $end = 3;
  151. $n = 4;
  152. $tol = 6;
  153. $expected = 21;
  154. // When
  155. $x = TrapezoidalRule::approximate($func, $start, $end, $n);
  156. // Then
  157. $this->assertEqualsWithDelta($expected, $x, $tol);
  158. }
  159. /**
  160. * @test approximate e^x² (http://tutorial.math.lamar.edu/Classes/CalcII/ApproximatingDefIntegrals.aspx)
  161. * @dataProvider dataProviderForEXSquared
  162. * @param int $n
  163. * @param float $expected
  164. * @param float $tol
  165. * @throws \Exception
  166. */
  167. public function testApproximateUsingCallback2(int $n, float $expected, float $tol)
  168. {
  169. // Given e^x²
  170. $func = function ($x) {
  171. return M_E ** ($x ** 2);
  172. };
  173. $start = 0;
  174. $end = 2;
  175. // When
  176. $x = TrapezoidalRule::approximate($func, $start, $end, $n);
  177. // Then
  178. $this->assertEqualsWithDelta($expected, $x, $tol);
  179. }
  180. /**
  181. * http://tutorial.math.lamar.edu/Classes/CalcII/ApproximatingDefIntegrals.aspx
  182. * @return array (n, expected, tol)
  183. */
  184. public function dataProviderForEXSquared(): array
  185. {
  186. return [
  187. [4, 20.64455905, 4.19193129],
  188. [8, 17.5650858, 1.1124580],
  189. [16, 16.7353812, 0.2827535],
  190. [32, 16.5236176, 0.0709898],
  191. [64, 16.4703942, 0.0177665],
  192. [128, 16.4570706, 0.0044428],
  193. ];
  194. }
  195. /**
  196. * @test approximate using polynomial
  197. * @throws \Exception
  198. *
  199. * f(x) = x² + 2x + 1
  200. * Antiderivative F(x) = (1/3)x³ + x² + x
  201. * Indefinite integral over [0, 3] = F(3) - F(0) = 21
  202. *
  203. * h₁, h₂, ... denotes the size on interval 1, 2, ...
  204. * ζ₁, ζ₂, ... denotes the max of the second derivative of f(x) on
  205. * interval 1, 2, ...
  206. * f'(x) = 2x + 2
  207. * f''(x) = 2
  208. * ζ = f''(x) = 2
  209. * Error = Sum(ζ₁h₁³ + ζ₂h₂³ + ...) = 2 * Sum(h₁³ + h₂³ + ...)
  210. */
  211. public function testApproximateUsingPolynomial()
  212. {
  213. // Given x² + 2x + 1
  214. $polynomial = new Polynomial([1, 2, 1]);
  215. $start = 0;
  216. $end = 3;
  217. $n = 4;
  218. $tol = 6;
  219. $expected = 21;
  220. // When
  221. $x = TrapezoidalRule::approximate($polynomial, $start, $end, $n);
  222. // Then
  223. $this->assertEqualsWithDelta($expected, $x, $tol);
  224. }
  225. }