diff --git a/Math/Solver/Root/Bisection.php b/Math/Solver/Root/Bisection.php new file mode 100644 index 000000000..cc57970b4 --- /dev/null +++ b/Math/Solver/Root/Bisection.php @@ -0,0 +1,57 @@ += 0) { + throw new \Exception("Function values at endpoints must have opposite signs."); + } + + $iteration = 0; + while (($b - $a) / 2 > self::EPSILON && $iteration < $maxIterations) { + $c = ($a + $b) / 2; + + $y = $func($c); + + if ($y === 0.0) { + return $c; + } + + if ($y * $func($a) < 0) { + $b = $c; + } else { + $a = $c; + } + + ++$iteration; + } + + return ($a + $b) / 2; + } +} diff --git a/Math/Solver/Root/Illinois.php b/Math/Solver/Root/Illinois.php new file mode 100644 index 000000000..fac0eeff5 --- /dev/null +++ b/Math/Solver/Root/Illinois.php @@ -0,0 +1,71 @@ += 0) { + throw new \Exception("Function values at endpoints must have opposite signs."); + } + + $c = $b; + $iteration = 0; + $sign = 1; + + while (($y = \abs($func($c))) > self::EPSILON && $iteration < $maxIterations) { + $fa = $func($a); + $fb = $func($b); + + if ($y === 0.0) { + return $c; + } + + // @todo: c might be wrong, could be that if and else must be switched + // @see https://en.wikipedia.org/wiki/Regula_falsi#The_Illinois_algorithm + if ($y * $fa < 0) { + $c = $sign === (int) ($y >= 0) + ? (0.5 * $a * $fb - $b * $fa) / (0.5 * $fb - $fa) + : ($a * $fb - $b * $fa) / ($fb - $fa); + + $b = $c; + } else { + $c = $sign === (int) ($y >= 0) + ? ($a * $fb - 0.5 * $b * $fa) / ($fb - 0.5 * $fa) + : ($a * $fb - $b * $fa) / ($fb - $fa); + + $a = $c; + } + + $sign = (int) ($y > 0); + + ++$iteration; + } + + return ($a + $b) / 2; + } +} diff --git a/Math/Solver/Root/RegulaFalsi.php b/Math/Solver/Root/RegulaFalsi.php new file mode 100644 index 000000000..b987ceb6c --- /dev/null +++ b/Math/Solver/Root/RegulaFalsi.php @@ -0,0 +1,60 @@ += 0) { + throw new \Exception("Function values at endpoints must have opposite signs."); + } + + $c = $b; + $iteration = 0; + + while (($y = \abs($func($c))) > self::EPSILON && $iteration < $maxIterations) { + $fa = $func($a); + $fb = $func($b); + + $c = ($a * $fb - $b * $fa) / ($fb - $fa); + + if ($y === 0.0) { + return $c; + } + + if ($y * $fa < 0) { + $b = $c; + } else { + $a = $c; + } + + ++$iteration; + } + + return ($a + $b) / 2; + } +}