Back
C
/******************************************************************************
* LICENSE *
******************************************************************************
* This file is part of mitx_mathematics_programming_examples. *
* *
* mitx_mathematics_programming_examples is free software: you can *
* redistribute it and/or modify it under the terms of the GNU General *
* Public License as published by the Free Software Foundation, either *
* version 3 of the License, or (at your option) any later version. *
* *
* mitx_mathematics_programming_examples is distributed in the hope that *
* it will be useful but WITHOUT ANY WARRANTY; without even the implied *
* warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. *
* See the GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with mitx_mathematics_programming_examples. If not, see *
* <https://www.gnu.org/licenses/>. *
******************************************************************************
* Purpose: *
* Calculates the root of a function using bisection. *
******************************************************************************
* Author: Ryan Maguire *
* Date: 2025/03/09 *
******************************************************************************/
/* stdio.h provides the "printf" function, used for printing text. */
#include <stdio.h>
/* Floating-point absolute value function, fabs, provided here. */
#include <math.h>
/* Function pointer notation is a little confusing. Create a typedef for it *
* so we do not need to explicitly use it later. */
typedef double (*function)(double);
/* Computes the root of a function using the bisection method. */
static double bisection_method(function f, double a, double b)
{
/* Tell the algorithm to stop after several iterations to avoid an *
* infinite loop. Double precision numbers have 52 bits in the mantissa, *
* so if |b - a| ~= 1, after 52 iterations of bisection we will get as *
* close as we can to the root. To allow for |b - a| to be larger, halt *
* the algorithm after at most 64 steps. */
const unsigned int maximum_number_of_iterations = 64U;
/* Getting exact roots is hard using floating-point numbers. Allow a *
* tolerance in our computation. This value is double precision epsilon. */
const double epsilon = 2.220446049250313E-16;
/* Variable for keeping track of how many iterations we have performed. */
unsigned int iters;
/* The midpoint for the bisection method. This will update as we iterate.*/
double midpoint;
/* We do not require a < b, nor do we require f(a) < f(b). We only need *
* one of these to evaluate to a negative under f and one to evaluate to *
* positive. We will call the negative entry left and positive one right.*/
double left, right;
/* Evaluate f at the two endpoints to determine which is positive and *
* which is negative. We transform [a, b] to [left, right] by doing this.*/
const double a_eval = f(a);
const double b_eval = f(b);
/* Rare case, f(a) = 0. Return a, no bisection needed. */
if (a_eval == 0.0)
return a;
/* Similarly, if f(b) = 0, then we have already found the root. Return b.*/
if (b_eval == 0.0)
return b;
/* Compare the two evaluations and set left and right accordingly. */
if (a_eval < b_eval)
{
/* If both evaluations are negative, or if both are positive, then *
* the bisection method will not work. Return NaN. */
if (b_eval < 0.0 || a_eval > 0.0)
return (a - a) / (a - a);
/* Otherwise, since f(a) < f(b), set left = a and right = b. */
left = a;
right = b;
}
/* In this case the function starts positive and tends to a negative. */
else
{
/* Same sanity check as before. We need one evaluation to be *
* negative and one to be positive. Abort if both have the same sign.*/
if (a_eval < 0.0 || b_eval > 0.0)
return (a - a) / (a - a);
/* Since f(a) > f(b), set left = b and right = a. */
left = b;
right = a;
}
/* Start the bisection method. Compute the midpoint of a and b. */
midpoint = 0.5 * (a + b);
/* Iteratively divide the range in half to find the root. */
for (iters = 0; iters < maximum_number_of_iterations; ++iters)
{
/* If f(x) is very small, we are close to a root and can break out *
* of this for loop. Check for this. */
const double eval = f(midpoint);
if (fabs(eval) <= epsilon)
break;
/* Apply bisection to get a better approximation for the root. We *
* have f(left) < 0 < f(right). If f(midpoint) < 0, replace the *
* interval [left, right] with [midpoint, right]. Set left to the *
* midpoint and reset the midpoint to be closer to right. */
if (eval < 0.0)
{
left = midpoint;
midpoint = 0.5 * (midpoint + right);
}
/* In the other case, f(midpoint) > 0, we replace right with the *
* midpoint, changing [left, right] into [left, midpoint]. We then *
* set the midpoint to be closer to left. */
else
{
right = midpoint;
midpoint = 0.5 * (left + midpoint);
}
}
/* After n iterations, we are no more than |b - a| / 2^n away from the *
* root of the function. 1 / 2^n goes to zero very quickly, meaning the *
* convergence is very quick. */
return midpoint;
}
/* End of bisection_method. */
/* Main routine used for testing our implementation of the bisection method. */
int main(void)
{
/* pi is somewhere between 3 and 4, and it is a root to sine. */
const double a = 3.0;
const double b = 4.0;
/* Compute pi using bisection. We should get pi = 3.14159..., accurate *
* to about 16 decimals. */
const double pi = bisection_method(sin, a, b);
printf("pi = %.16f\n", pi);
return 0;
}
C3
/******************************************************************************
* LICENSE *
******************************************************************************
* This file is part of mitx_mathematics_programming_examples. *
* *
* mitx_mathematics_programming_examples is free software: you can *
* redistribute it and/or modify it under the terms of the GNU General *
* Public License as published by the Free Software Foundation, either *
* version 3 of the License, or (at your option) any later version. *
* *
* mitx_mathematics_programming_examples is distributed in the hope that *
* it will be useful but WITHOUT ANY WARRANTY; without even the implied *
* warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. *
* See the GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with mitx_mathematics_programming_examples. If not, see *
* <https://www.gnu.org/licenses/>. *
******************************************************************************
* Purpose: *
* Calculates the root of a function using bisection. *
******************************************************************************
* Author: Ryan Maguire *
* Date: 2025/05/27 *
******************************************************************************/
/* std::io provides the "printfn" function, used for printing text. */
import std::io;
/* Floating-point absolute value function, abs, provided here. */
import std::math;
/* Function pointer notation is a little confusing. Create a typedef for it *
* so we do not need to explicitly use it later. */
alias RealFunc = fn double(double);
/* Computes the root of a function using the bisection method. */
fn double bisection_method(RealFunc f, double a, double b)
{
/* Tell the algorithm to stop after several iterations to avoid an *
* infinite loop. Double precision numbers have 52 bits in the mantissa, *
* so if |b - a| ~= 1, after 52 iterations of bisection we will get as *
* close as we can to the root. To allow for |b - a| to be larger, halt *
* the algorithm after at most 64 steps. */
const uint MAXIMUM_NUMBER_OF_ITERATIONS = 64U;
/* Getting exact roots is hard using floating-point numbers. Allow a *
* tolerance in our computation. This value is double precision epsilon. */
const double EPSILON = 2.220446049250313E-16;
/* Variable for keeping track of how many iterations we have performed. */
uint iters;
/* The midpoint for the bisection method. This will update as we iterate.*/
double midpoint;
/* We do not require a < b, nor do we require f(a) < f(b). We only need *
* one of these to evaluate to a negative under f and one to evaluate to *
* positive. We will call the negative entry left and positive one right.*/
double left, right;
/* Evaluate f at the two endpoints to determine which is positive and *
* which is negative. We transform [a, b] to [left, right] by doing this.*/
double a_eval = f(a);
double b_eval = f(b);
/* Rare case, f(a) = 0. Return a, no bisection needed. */
if (a_eval == 0.0) {
return a;
}
/* Similarly, if f(b) = 0, then we have already found the root. Return b.*/
if (b_eval == 0.0) {
return b;
}
/* Compare the two evaluations and set left and right accordingly. */
if (a_eval < b_eval)
{
/* If both evaluations are negative, or if both are positive, then *
* the bisection method will not work. Return NaN. */
if (b_eval < 0.0 || a_eval > 0.0) {
return (a - a) / (a - a);
}
/* Otherwise, since f(a) < f(b), set left = a and right = b. */
left = a;
right = b;
}
/* In this case the function starts positive and tends to a negative. */
else
{
/* Same sanity check as before. We need one evaluation to be *
* negative and one to be positive. Abort if both have the same sign.*/
if (a_eval < 0.0 || b_eval > 0.0) {
return (a - a) / (a - a);
}
/* Since f(a) > f(b), set left = b and right = a. */
left = b;
right = a;
}
/* Start the bisection method. Compute the midpoint of a and b. */
midpoint = 0.5 * (a + b);
/* Iteratively divide the range in half to find the root. */
for (iters = 0; iters < MAXIMUM_NUMBER_OF_ITERATIONS; ++iters)
{
/* If f(x) is very small, we are close to a root and can break out *
* of this for loop. Check for this. */
double eval = f(midpoint);
if (math::abs(eval) <= EPSILON) {
break;
}
/* Apply bisection to get a better approximation for the root. We *
* have f(left) < 0 < f(right). If f(midpoint) < 0, replace the *
* interval [left, right] with [midpoint, right]. Set left to the *
* midpoint and reset the midpoint to be closer to right. */
if (eval < 0.0)
{
left = midpoint;
midpoint = 0.5 * (midpoint + right);
}
/* In the other case, f(midpoint) > 0, we replace right with the *
* midpoint, changing [left, right] into [left, midpoint]. We then *
* set the midpoint to be closer to left. */
else
{
right = midpoint;
midpoint = 0.5 * (left + midpoint);
}
}
/* After n iterations, we are no more than |b - a| / 2^n away from the *
* root of the function. 1 / 2^n goes to zero very quickly, meaning the *
* convergence is very quick. */
return midpoint;
}
/* End of bisection_method. */
/* math::sin is a macro, which we cannot take the address of. C3 allows *
* linking with the C standard library, including libm. Tell the compiler *
* about the C routine so that we may call it from our main function. */
extern fn double sin(double x) @extern("sin");
/* Main routine used for testing our implementation of the bisection method. */
fn void main()
{
/* pi is somewhere between 3 and 4, and it is a root to sine. */
const double A = 3.0;
const double B = 4.0;
/* Compute pi using bisection. We should get pi = 3.14159..., accurate *
* to about 16 decimals. */
double pi = bisection_method(&sin, A, B);
io::printfn("pi = %.16f", pi);
}
C++
/******************************************************************************
* LICENSE *
******************************************************************************
* This file is part of mitx_mathematics_programming_examples. *
* *
* mitx_mathematics_programming_examples is free software: you can *
* redistribute it and/or modify it under the terms of the GNU General *
* Public License as published by the Free Software Foundation, either *
* version 3 of the License, or (at your option) any later version. *
* *
* mitx_mathematics_programming_examples is distributed in the hope that *
* it will be useful but WITHOUT ANY WARRANTY; without even the implied *
* warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. *
* See the GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with mitx_mathematics_programming_examples. If not, see *
* <https://www.gnu.org/licenses/>. *
******************************************************************************
* Purpose: *
* Calculates the root of a function using the bisection method. *
******************************************************************************
* Author: Ryan Maguire *
* Date: 2025/03/09 *
******************************************************************************/
/* stdio provides the "printf" function, used for printing text. */
#include <cstdio>
/* Floating-point absolute value function, fabs, provided here. */
#include <cmath>
/* Function pointer notation is a little confusing. Create a typedef for it *
* so we do not need to explicitly use it later. */
typedef double (*function)(double);
/* Class providing an implementation of the bisection method. */
class Bisection {
/* The error after n iterations is |b - a| / 2^n. Since double has a *
* 52-bit mantissa, if |b - a| ~= 1, then after 52 steps we can halt the *
* program. To allow for |b - a| to be larger, we stop the process after *
* at most 64 iterations. */
static const unsigned int maximum_number_of_iterations = 64U;
/* We want the function visible outside the class. Declare it public. */
public:
/* Computes the root of a function using the bisection method. */
static double root(function f, double a, double b)
{
/* The maximum allowed error. This is double precision epsilon. */
const double epsilon = 2.220446049250313E-16;
/* Variable for keeping track of the number of iterations. */
unsigned int iters;
/* The midpoint for the bisection. This updates as we iterate. */
double midpoint;
/* We do not require a < b, nor f(a) < f(b). We need one of *
* these to evaluate negative under f and one to evaluate to *
* positive. Call the negative entry left and positive one right.*/
double left, right;
/* Evaluate f at the endpoints to determine which is positive *
* and which is negative, transforming [a, b] to [left, right]. */
const double a_eval = f(a);
const double b_eval = f(b);
/* Rare case, f(a) = 0. Return a, no bisection needed. */
if (a_eval == 0.0)
return a;
/* Similarly, if f(b) = 0, then we found the root. Return b. */
if (b_eval == 0.0)
return b;
/* Compare the two evaluations and set the left and right values.*/
if (a_eval < b_eval)
{
/* If both evaluations are negative, or if both are positive,*
* then the bisection method will not work. Return NaN. */
if (b_eval < 0.0 || a_eval > 0.0)
return (a - a) / (a - a);
/* Otherwise, since f(a) < f(b), set left = a and right = b. */
left = a;
right = b;
}
/* In this case the function starts positive and goes negative. */
else
{
/* Same sanity check as before. We need one evaluation to be *
* negative and one to be positive. Abort if the signs agree.*/
if (a_eval < 0.0 || b_eval > 0.0)
return (a - a) / (a - a);
/* Since f(a) > f(b), set left = b and right = a. */
left = b;
right = a;
}
/* Start the bisection method. Compute the midpoint of a and b. */
midpoint = 0.5 * (a + b);
/* Iteratively divide the range in half to find the root. */
for (iters = 0; iters < maximum_number_of_iterations; ++iters)
{
/* If f(x) is very small, we are close to a root and can *
* break out of this for loop. Check for this. */
const double eval = f(midpoint);
if (std::fabs(eval) <= epsilon)
break;
/* Apply bisection to get a better approximation. We have *
* f(left) < 0 < f(right). If f(midpoint) < 0, replace the *
* interval [left, right] with [midpoint, right]. Set left *
* to the midpoint and set midpoint to be closer to right. */
if (eval < 0.0)
{
left = midpoint;
midpoint = 0.5 * (midpoint + right);
}
/* If f(midpoint) > 0, then replace right with the midpoint, *
* changing [left, right] into [left, midpoint]. We then set *
* the midpoint to be closer to left. */
else
{
right = midpoint;
midpoint = 0.5 * (left + midpoint);
}
}
/* After n iterations, we are at most |b - a| / 2^n from the *
* root of the function. 1 / 2^n goes to zero very quickly, *
* meaning the convergence is very quick. */
return midpoint;
}
/* End of root. */
};
/* End of Bisection definition. */
/* Main routine used for testing our implementation of the Bisection method. */
int main(void)
{
/* pi is somewhere between 3 and 4, and it is a root to sine. */
const double a = 3.0;
const double b = 4.0;
/* Compute pi using bisection. We should get pi = 3.14159..., accurate *
* to about 16 decimals. */
const double pi = Bisection::root(std::sin, a, b);
std::printf("pi = %.16f\n", pi);
return 0;
}
C#
/******************************************************************************
* LICENSE *
******************************************************************************
* This file is part of mitx_mathematics_programming_examples. *
* *
* mitx_mathematics_programming_examples is free software: you can *
* redistribute it and/or modify it under the terms of the GNU General *
* Public License as published by the Free Software Foundation, either *
* version 3 of the License, or (at your option) any later version. *
* *
* mitx_mathematics_programming_examples is distributed in the hope that *
* it will be useful but WITHOUT ANY WARRANTY; without even the implied *
* warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. *
* See the GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with mitx_mathematics_programming_examples. If not, see *
* <https://www.gnu.org/licenses/>. *
******************************************************************************
* Purpose: *
* Calculates the root of a function using the bisection method. *
******************************************************************************
* Author: Ryan Maguire *
* Date: 2025/03/28 *
******************************************************************************/
/* Console.WriteLine and Math.Abs are both provided here. */
using System;
/* Create an alias for functions of the form f: R -> R. */
using function = System.Func<double, double>;
/* Class providing an implementation of the bisection method. */
class Bisection {
/* The error after n iterations is |b - a| / 2^n. Since double has a *
* 52-bit mantissa, if |b - a| ~= 1, then after 52 steps we can halt the *
* program. To allow for |b - a| to be larger, we stop the process after *
* at most 64 iterations. */
const uint maximumNumberOfIterations = 64U;
/* The maximum allowed error. This is double precision epsilon. */
const double epsilon = 2.220446049250313E-16;
/* Computes the root of a function using the bisection method. */
static double Root(function f, double a, double b)
{
/* Variable for keeping track of the number of iterations. */
uint iters;
/* The midpoint for the bisection. This updates as we iterate. */
double midpoint;
/* We do not require a < b, nor f(a) < f(b). We need one of *
* these to evaluate negative under f and one to evaluate to *
* positive. Call the negative entry left and positive one right. */
double left, right;
/* Evaluate f at the endpoints to determine which is positive *
* and which is negative, transforming [a, b] to [left, right]. */
double aEval = f(a);
double bEval = f(b);
/* Rare case, f(a) = 0. Return a, no bisection needed. */
if (aEval == 0.0)
return a;
/* Similarly, if f(b) = 0, then we found the root. Return b. */
if (bEval == 0.0)
return b;
/* Compare the two evaluations and set the left and right values. */
if (aEval < bEval)
{
/* If both evaluations are negative, or if both are positive, *
* then the bisection method will not work. Return NaN. */
if (bEval < 0.0 || aEval > 0.0)
return (a - a) / (a - a);
/* Otherwise, since f(a) < f(b), set left = a and right = b. */
left = a;
right = b;
}
/* In this case the function starts positive and goes negative. */
else
{
/* Same sanity check as before. We need one evaluation to be *
* negative and one to be positive. Abort if the signs agree. */
if (aEval < 0.0 || bEval > 0.0)
return (a - a) / (a - a);
/* Since f(a) > f(b), set left = b and right = a. */
left = b;
right = a;
}
/* Start the bisection method. Compute the midpoint of a and b. */
midpoint = 0.5 * (a + b);
/* Iteratively divide the range in half to find the root. */
for (iters = 0; iters < maximumNumberOfIterations; ++iters)
{
/* If f(x) is very small, we are close to a root and can *
* break out of this for loop. Check for this. */
double eval = f(midpoint);
if (Math.Abs(eval) <= epsilon)
break;
/* Apply bisection to get a better approximation. We have *
* f(left) < 0 < f(right). If f(midpoint) < 0, replace the *
* interval [left, right] with [midpoint, right]. Set left *
* to the midpoint and set midpoint to be closer to right. */
if (eval < 0.0)
{
left = midpoint;
midpoint = 0.5 * (midpoint + right);
}
/* If f(midpoint) > 0, then replace right with the midpoint, *
* changing [left, right] into [left, midpoint]. We then set *
* the midpoint to be closer to left. */
else
{
right = midpoint;
midpoint = 0.5 * (left + midpoint);
}
}
/* After n iterations, we are at most |b - a| / 2^n from the *
* root of the function. 1 / 2^n goes to zero very quickly, *
* meaning the convergence is very quick. */
return midpoint;
}
/* End of root. */
/* Main routine used for testing our Bisection implementation. */
static void Main()
{
/* pi is somewhere between 3 and 4, and it is a root to sine. */
const double a = 3.0;
const double b = 4.0;
/* Compute pi using bisection. We should get pi = 3.14159..., *
* accurate to about 16 decimals. */
double pi = Root(Math.Sin, a, b);
Console.WriteLine($"pi = {pi}");
}
}
D
/******************************************************************************
* LICENSE *
******************************************************************************
* This file is part of mitx_mathematics_programming_examples. *
* *
* mitx_mathematics_programming_examples is free software: you can *
* redistribute it and/or modify it under the terms of the GNU General *
* Public License as published by the Free Software Foundation, either *
* version 3 of the License, or (at your option) any later version. *
* *
* mitx_mathematics_programming_examples is distributed in the hope that *
* it will be useful but WITHOUT ANY WARRANTY; without even the implied *
* warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. *
* See the GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with mitx_mathematics_programming_examples. If not, see *
* <https://www.gnu.org/licenses/>. *
******************************************************************************
* Purpose: *
* Calculates the root of a function using bisection. *
******************************************************************************
* Author: Ryan Maguire *
* Date: 2025/03/09 *
******************************************************************************/
/* Create an alias for real-valued functions f: R -> R. */
alias realfunc = double function(double) pure nothrow @nogc @safe;
/* Computes the root of a function using the bisection method. */
static double bisection_method(realfunc f, double a, double b)
{
/* Floating-point absolute value (fabs) provided here. */
import std.math : fabs;
/* Tell the algorithm to stop after several iterations to avoid an *
* infinite loop. Double precision numbers have 52 bits in the mantissa, *
* so if |b - a| ~= 1, after 52 iterations of bisection we will get as *
* close as we can to the root. To allow for |b - a| to be larger, halt *
* the algorithm after at most 64 steps. */
const uint maximum_number_of_iterations = 64U;
/* Getting exact roots is hard using floating-point numbers. Allow a *
* tolerance in our computation. This value is double precision epsilon. */
const double epsilon = 2.220446049250313E-16;
/* Variable for keeping track of how many iterations we have performed. */
uint iters;
/* The midpoint for the bisection method. This will update as we iterate.*/
double midpoint;
/* We do not require a < b, nor do we require f(a) < f(b). We only need *
* one of these to evaluate to a negative under f and one to evaluate to *
* positive. We will call the negative entry left and positive one right.*/
double left, right;
/* Evaluate f at the two endpoints to determine which is positive and *
* which is negative. We transform [a, b] to [left, right] by doing this.*/
const double a_eval = f(a);
const double b_eval = f(b);
/* Rare case, f(a) = 0. Return a, no bisection needed. */
if (a_eval == 0.0)
return a;
/* Similarly, if f(b) = b, then we have already found the root. Return b.*/
if (b_eval == 0.0)
return b;
/* Compare the two evaluations and set left and right accordingly. */
if (a_eval < b_eval)
{
/* If both evaluations are negative, or if both are positive, then *
* the bisection method will not work. Return NaN. */
if (b_eval < 0.0 || a_eval > 0.0)
return (a - a) / (a - a);
/* Otherwise, since f(a) < f(b), set left = a and right = b. */
left = a;
right = b;
}
/* In this case the function starts positive and tends to a negative. */
else
{
/* Same sanity check as before. We need one evaluation to be *
* negative and one to be positive. Abort if both have the same sign.*/
if (a_eval < 0.0 || b_eval > 0.0)
return (a - a) / (a - a);
/* Since f(a) > f(b), set left = b and right = a. */
left = b;
right = a;
}
/* Start the bisection method. Compute the midpoint of a and b. */
midpoint = 0.5 * (a + b);
/* Iteratively divide the range in half to find the root. */
for (iters = 0; iters < maximum_number_of_iterations; ++iters)
{
/* If f(x) is very small, we are close to a root and can break out *
* of this for loop. Check for this. */
const double eval = f(midpoint);
if (fabs(eval) <= epsilon)
break;
/* Apply bisection to get a better approximation for the root. We *
* have f(left) < 0 < f(right). If f(midpoint) < 0, replace the *
* interval [left, right] with [midpoint, right]. Set left to the *
* midpoint and reset the midpoint to be closer to right. */
if (eval < 0.0)
{
left = midpoint;
midpoint = 0.5 * (midpoint + right);
}
/* In the other case, f(midpoint) > 0, we replace right with the *
* midpoint, changing [left, right] into [left, midpoint]. We then *
* set the midpoint to be closer to left. */
else
{
right = midpoint;
midpoint = 0.5 * (left + midpoint);
}
}
/* After n iterations, we are no more than |b - a| / 2^n away from the *
* root of the function. 1 / 2^n goes to zero very quickly, meaning the *
* convergence is very quick. */
return midpoint;
}
/* End of bisection_method. */
/* Main routine used for testing our implementation of the bisection method. */
int main()
{
/* stdio contains "printf", used for printing text. */
import std.stdio : printf;
/* We will compute pi via bisection of sine on the interval [3, 4]. */
import std.math : sin;
/* pi is somewhere between 3 and 4, and it is a root to sine. */
const double a = 3.0;
const double b = 4.0;
/* Compute pi using bisection. We should get pi = 3.14159..., accurate *
* to about 16 decimals. */
const double pi = bisection_method(&sin, a, b);
printf("pi = %.16f\n", pi);
return 0;
}
Fortran
!------------------------------------------------------------------------------!
! LICENSE !
!------------------------------------------------------------------------------!
! This file is part of mitx_mathematics_programming_examples. !
! !
! mitx_mathematics_programming_examples is free software: you can !
! redistribute it and/or modify it under the terms of the GNU General Public !
! License as published by the Free Software Foundation, either version 3 of !
! the License, or (at your option) any later version. !
! !
! mitx_mathematics_programming_examples is distributed in the hope that it !
! will be useful but WITHOUT ANY WARRANTY; without even the implied warranty !
! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !
! GNU General Public License for more details. !
! !
! You should have received a copy of the GNU General Public License !
! along with mitx_mathematics_programming_examples. If not, see !
! <https://www.gnu.org/licenses/>. !
!------------------------------------------------------------------------------!
! Purpose: !
! Calculates roots using the bisection method. !
!------------------------------------------------------------------------------!
! Author: Ryan Maguire !
! Date: 2025/03/09 !
!------------------------------------------------------------------------------!
MODULE BISECTION
IMPLICIT NONE
! Given a continuous function defined on [a, b], the nth iteration of the
! bisection method is at most |b - a| / 2^n away from the root. A double
! precision number has 52 bits in the mantissa, meaning if |b - a| ~= 1,
! then after 52 iterations we are as close to the root as we can get.
! To allow for a larger range, halt the algorithm after 64 iterations.
INTEGER :: MAXIMUM_NUMBER_OF_ITERATIONS = 64
! Maximum allowed error. This is double precision epsilon.
REAL :: EPSILON = 2.220446049250313E-16
CONTAINS
!--------------------------------------------------------------------------!
! Function: !
! BISECTION_METHOD !
! Purpose: !
! Computes roots using the bisection method. !
! Arguments: !
! FUNC (FUNCTION): !
! A function from [A, B] to the real numbers. !
! A (REAL): !
! One of the endpoints for the interval. !
! B (REAL): !
! The other the endpoint for the interval. !
! OUTPUT: !
! ROOT (REAL): !
! A root for FUNC. !
!--------------------------------------------------------------------------!
FUNCTION BISECTION_METHOD(FUNC, A, B)
IMPLICIT NONE
! FUNC is a function f: [A, B] -> Reals. Define this.
INTERFACE
FUNCTION FUNC(X)
IMPLICIT NONE
REAL, INTENT(IN) :: X
REAL :: FUNC
END FUNCTION FUNC
END INTERFACE
! The inputs are positive real numbers.
REAL, INTENT(IN) :: A, B
! The output is also a positive real number, the root of FUNC.
REAL :: BISECTION_METHOD
! Dummy variable for keeping track of how many iterations we've done.
INTEGER :: ITERS
! The midpoint of the interval. The bisection method iteratively cuts
! the interval in half to find the root.
REAL :: MIDPOINT
! The endpoints of the current interval. These are updated with each
! iteration. We initially define LEFT so that FUNC(LEFT) < 0, and
! similarly choose RIGHT to make 0 < FUNC(RIGHT).
REAL :: LEFT, RIGHT
! The evaluations of A and B are used to determine LEFT and RIGHT.
REAL :: A_EVAL, B_EVAL
! At each iteration we compute FUNC(MIDPOINT) and reset the interval
! depending on whether or not the value is positive.
REAL :: EVAL
! Initial setup, find out which evaluation is positive and which is
! negative. If the signs of the two agree we treat this as an error.
A_EVAL = FUNC(A)
B_EVAL = FUNC(B)
! Special case, A is a root. Set the output to A and return.
IF (ABS(A_EVAL) .LE. EPSILON) THEN
BISECTION_METHOD = A
RETURN
END IF
! Similarly for B, if FUNC(B) ~= 0, return B.
IF (ABS(B_EVAL) .LE. EPSILON) THEN
BISECTION_METHOD = B
RETURN
END IF
! If FUNC(A) < 0 < FUNC(B), then LEFT = A and RIGHT = B.
IF (A_EVAL .LT. B_EVAL) THEN
! If FUNC(A) and FUNC(B) have the same sign (both positive or
! both negative), return NaN. Bisection is undefined.
IF ((A_EVAL .GT. 0) .OR. (B_EVAL .LT. 0)) THEN
BISECTION_METHOD = (A - A) / (A - A)
RETURN
END IF
LEFT = A
RIGHT = B
! If FUNC(B) < 0 < FUNC(A), then LEFT = B and RIGHT = A.
ELSE
! Same sanity check as before, make sure the signs are different.
IF ((A_EVAL .LT. 0) .OR. (B_EVAL .GT. 0)) THEN
BISECTION_METHOD = (A - A) / (A - A)
RETURN
END IF
LEFT = B
RIGHT = A
END IF
MIDPOINT = 0.5 * (A + B)
! Iteratively perform the bisection method.
DO ITERS = 0, MAXIMUM_NUMBER_OF_ITERATIONS
! The interval is cut in half based on the sign of FUNC(MIDPOINT).
! Compute this and compare.
EVAL = FUNC(MIDPOINT)
! If MIDPOINT is close to a root we can exit the function.
IF (ABS(EVAL) .LE. EPSILON) THEN
BISECTION_METHOD = MIDPOINT
RETURN
END IF
! Otherwise, divide the range in half. If EVAL is negative we
! replace LEFT with MIDPOINT and move MIDPOINT closer to RIGHT.
IF (EVAL .LT. 0) THEN
LEFT = MIDPOINT
MIDPOINT = 0.5 * (MIDPOINT + RIGHT)
! Likewise, if EVAL is positive, replace RIGHT with MIDPOINT and
! move MIDPOINT closer to LEFT.
ELSE
RIGHT = MIDPOINT
MIDPOINT = 0.5 * (LEFT + MIDPOINT)
END IF
END DO
! Provided |B - A| is not too big, we should now have a very good
! approximation to the root.
BISECTION_METHOD = MIDPOINT
END FUNCTION BISECTION_METHOD
END MODULE BISECTION
! Program for testing our implementation of the bisection method.
PROGRAM MAIN
USE BISECTION
IMPLICIT NONE
! We will compute pi by finding a root to sine on the interval [3, 4].
! The sine function is an intrinsic provided by Fortran.
INTRINSIC SIN
! Pi is somewhere between 3 and 4, use this interval.
REAL :: A = 3.0
REAL :: B = 4.0
! Variable for the output.
REAL :: PI
! Run the routine, compute pi, and print the result.
PI = BISECTION_METHOD(SIN, A, B)
PRINT "(A,F18.16)", "pi = ", PI
END PROGRAM MAIN
Go
/******************************************************************************
* LICENSE *
******************************************************************************
* This file is part of mitx_mathematics_programming_examples. *
* *
* mitx_mathematics_programming_examples is free software: you can *
* redistribute it and/or modify it under the terms of the GNU General *
* Public License as published by the Free Software Foundation, either *
* version 3 of the License, or (at your option) any later version. *
* *
* mitx_mathematics_programming_examples is distributed in the hope that *
* it will be useful but WITHOUT ANY WARRANTY; without even the implied *
* warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. *
* See the GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with mitx_mathematics_programming_examples. If not, see *
* <https://www.gnu.org/licenses/>. *
******************************************************************************
* Purpose: *
* Calculates the root of a function using bisection. *
******************************************************************************
* Author: Ryan Maguire *
* Date: 2025/03/28 *
******************************************************************************/
package main
/* Only standard library imports are needed. */
import (
"fmt" /* Printf provided here, used for printing text to the screen. */
"math" /* Abs, floating-point absolute value function, found here. */
)
/* Function pointer notation is a little confusing. Create a typedef for it *
* so we do not need to explicitly use it later. */
type realfunc func(x float64) float64
/* Computes the root of a function using the bisection method. */
func bisection_method(f realfunc, a float64, b float64) float64 {
/* Tell the algorithm to stop after several iterations to avoid an *
* infinite loop. Double precision numbers have 52 bits in the mantissa, *
* so if |b - a| ~= 1, after 52 iterations of bisection we will get as *
* close as we can to the root. To allow for |b - a| to be larger, halt *
* the algorithm after at most 64 steps. */
const maximum_number_of_iterations uint32 = 64
/* Getting exact roots is hard using floating-point numbers. Allow a *
* tolerance in our computation. This value is double precision epsilon. */
const epsilon float64 = 2.220446049250313E-16
/* Variable for keeping track of how many iterations we have performed. */
var iters uint32
/* The midpoint for the bisection method. This will update as we iterate.*/
var midpoint float64
/* We do not require a < b, nor do we require f(a) < f(b). We only need *
* one of these to evaluate to a negative under f and one to evaluate to *
* positive. We will call the negative entry left and positive one right.*/
var left, right float64
/* Evaluate f at the two endpoints to determine which is positive and *
* which is negative. We transform [a, b] to [left, right] by doing this.*/
var a_eval = f(a)
var b_eval = f(b)
/* Rare case, f(a) = 0. Return a, no bisection needed. */
if a_eval == 0.0 {
return a
}
/* Similarly, if f(b) = 0, then we have already found the root. Return b.*/
if b_eval == 0.0 {
return b
}
/* Compare the two evaluations and set left and right accordingly. */
if a_eval < b_eval {
/* If both evaluations are negative, or if both are positive, then *
* the bisection method will not work. Return NaN. */
if b_eval < 0.0 || a_eval > 0.0 {
return (a - a) / (a - a)
}
/* Otherwise, since f(a) < f(b), set left = a and right = b. */
left = a
right = b
/* In this case the function starts positive and tends to a negative. */
} else {
/* Same sanity check as before. We need one evaluation to be *
* negative and one to be positive. Abort if both have the same sign.*/
if a_eval < 0.0 || b_eval > 0.0 {
return (a - a) / (a - a)
}
/* Since f(a) > f(b), set left = b and right = a. */
left = b
right = a
}
/* Start the bisection method. Compute the midpoint of a and b. */
midpoint = 0.5 * (a + b)
/* Iteratively divide the range in half to find the root. */
for iters = 0; iters < maximum_number_of_iterations; iters += 1 {
/* If f(x) is very small, we are close to a root and can break out *
* of this for loop. Check for this. */
var eval = f(midpoint)
if math.Abs(eval) <= epsilon {
break
}
/* Apply bisection to get a better approximation for the root. We *
* have f(left) < 0 < f(right). If f(midpoint) < 0, replace the *
* interval [left, right] with [midpoint, right]. Set left to the *
* midpoint and reset the midpoint to be closer to right. */
if eval < 0.0 {
left = midpoint
midpoint = 0.5 * (midpoint + right)
/* In the other case, f(midpoint) > 0, we replace right with the *
* midpoint, changing [left, right] into [left, midpoint]. We then *
* set the midpoint to be closer to left. */
} else {
right = midpoint
midpoint = 0.5 * (left + midpoint)
}
}
/* After n iterations, we are no more than |b - a| / 2^n away from the *
* root of the function. 1 / 2^n goes to zero very quickly, meaning the *
* convergence is very quick. */
return midpoint
}
/* End of bisection_method. */
/* Main routine used for testing our implementation of the bisection method. */
func main() {
/* pi is somewhere between 3 and 4, and it is a root to sine. */
const a float64 = 3.0
const b float64 = 4.0
/* Compute pi using bisection. We should get pi = 3.14159..., accurate *
* to about 16 decimals. */
var pi = bisection_method(math.Sin, a, b)
fmt.Printf("pi = %.16f\n", pi)
}
Java
/******************************************************************************
* LICENSE *
******************************************************************************
* This file is part of mitx_mathematics_programming_examples. *
* *
* mitx_mathematics_programming_examples is free software: you can *
* redistribute it and/or modify it under the terms of the GNU General *
* Public License as published by the Free Software Foundation, either *
* version 3 of the License, or (at your option) any later version. *
* *
* mitx_mathematics_programming_examples is distributed in the hope that *
* it will be useful but WITHOUT ANY WARRANTY; without even the implied *
* warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. *
* See the GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with mitx_mathematics_programming_examples. If not, see *
* <https://www.gnu.org/licenses/>. *
******************************************************************************
* Purpose: *
* Calculates the root of a function using the bisection method. *
******************************************************************************
* Author: Ryan Maguire *
* Date: 2025/03/28 *
******************************************************************************/
/* Class providing an implementation of the bisection method. */
final public class Bisection {
/* Interface used for defining functions of the type f: R -> R. */
interface Function
{
double evaluate(double x);
}
/* The error after n iterations is |b - a| / 2^n. Since double has a *
* 52-bit mantissa, if |b - a| ~= 1, then after 52 steps we can halt the *
* program. To allow for |b - a| to be larger, we stop the process after *
* at most 64 iterations. */
static private final int maximumNumberOfIterations = 64;
/* The maximum allowed error. This is double precision epsilon. */
static private final double epsilon = 2.220446049250313E-16;
/* Computes the root of a function using the bisection method. */
static double root(Function f, double a, double b)
{
/* Variable for keeping track of the number of iterations. */
int iters;
/* The midpoint for the bisection. This updates as we iterate. */
double midpoint;
/* We do not require a < b, nor f(a) < f(b). We need one of *
* these to evaluate negative under f and one to evaluate to *
* positive. Call the negative entry left and positive one right. */
double left, right;
/* Evaluate f at the endpoints to determine which is positive *
* and which is negative, transforming [a, b] to [left, right]. */
double aEval = f.evaluate(a);
double bEval = f.evaluate(b);
/* Rare case, f(a) = 0. Return a, no bisection needed. */
if (aEval == 0.0)
return a;
/* Similarly, if f(b) = 0, then we found the root. Return b. */
if (bEval == 0.0)
return b;
/* Compare the two evaluations and set the left and right values. */
if (aEval < bEval)
{
/* If both evaluations are negative, or if both are positive, *
* then the bisection method will not work. Return NaN. */
if (bEval < 0.0 || aEval > 0.0)
return (a - a) / (a - a);
/* Otherwise, since f(a) < f(b), set left = a and right = b. */
left = a;
right = b;
}
/* In this case the function starts positive and goes negative. */
else
{
/* Same sanity check as before. We need one evaluation to be *
* negative and one to be positive. Abort if the signs agree. */
if (aEval < 0.0 || bEval > 0.0)
return (a - a) / (a - a);
/* Since f(a) > f(b), set left = b and right = a. */
left = b;
right = a;
}
/* Start the bisection method. Compute the midpoint of a and b. */
midpoint = 0.5 * (a + b);
/* Iteratively divide the range in half to find the root. */
for (iters = 0; iters < maximumNumberOfIterations; ++iters)
{
/* If f(x) is very small, we are close to a root and can *
* break out of this for loop. Check for this. */
double eval = f.evaluate(midpoint);
if (Math.abs(eval) <= epsilon)
break;
/* Apply bisection to get a better approximation. We have *
* f(left) < 0 < f(right). If f(midpoint) < 0, replace the *
* interval [left, right] with [midpoint, right]. Set left *
* to the midpoint and set midpoint to be closer to right. */
if (eval < 0.0)
{
left = midpoint;
midpoint = 0.5 * (midpoint + right);
}
/* If f(midpoint) > 0, then replace right with the midpoint, *
* changing [left, right] into [left, midpoint]. We then set *
* the midpoint to be closer to left. */
else
{
right = midpoint;
midpoint = 0.5 * (left + midpoint);
}
}
/* After n iterations, we are at most |b - a| / 2^n from the *
* root of the function. 1 / 2^n goes to zero very quickly, *
* meaning the convergence is very quick. */
return midpoint;
}
/* End of root. */
/* Main routine used for testing our Bisection implementation. */
static public void main(String[] args)
{
/* pi is somewhere between 3 and 4, and it is a root to sine. */
final double a = 3.0;
final double b = 4.0;
/* Create a function that evaluates sine using Math.sin. */
Function sin = new Function()
{
/* The evaluator simply calls the standard library function. */
public double evaluate(double x)
{
return Math.sin(x);
}
};
/* Compute pi using bisection. We should get pi = 3.14159..., *
* accurate to about 16 decimals. */
double pi = Bisection.root(sin, a, b);
System.out.printf("pi = %.16f\n", pi);
}
}
Julia
################################################################################
# LICENSE #
################################################################################
# This file is part of mitx_mathematics_programming_examples. #
# #
# mitx_mathematics_programming_examples is free software: you can #
# redistribute it and/or modify it under the terms of the GNU General Public #
# License as published by the Free Software Foundation, either version 3 of #
# the License, or (at your option) any later version. #
# #
# mitx_mathematics_programming_examples is distributed in the hope that it #
# will be useful, but WITHOUT ANY WARRANTY; without even the implied #
# warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #
# GNU General Public License for more details. #
# #
# You should have received a copy of the GNU General Public License #
# along with mitx_mathematics_programming_examples. If not, see #
# <https://www.gnu.org/licenses/>. #
################################################################################
# Purpose: #
# Calculates roots of a function using bisection. #
################################################################################
# Author: Ryan Maguire #
# Date: March 29, 2025. #
################################################################################
################################################################################
# Function: #
# bisection #
# Purpose: #
# Computes the root of a function f between a and b. #
# Arguments: #
# f (function): #
# A function f: R -> R. The root of f is computed. #
# a (real): #
# One of the endpoints of the interval for f. #
# b (real): #
# The other endpoint for f. #
# Output: #
# root (real): #
# A root of the function f between a and b. #
################################################################################
function bisection(f, a, b)
# Tell the algorithm to stop after several iterations to avoid an
# infinite loop. Double precision numbers have 52 bits in the mantissa,
# so if |b - a| ~= 1, after 52 iterations of bisection we will get as
# close as we can to the root. To allow for |b - a| to be larger, halt
# the algorithm after at most 64 steps.
maximum_number_of_iterations = 64
# The maximum allowed error. This is double precision epsilon.
epsilon = 2.220446049250313E-16
# Evaluate f at the two endpoints to determine which is positive and
# which is negative. We transform [a, b] to [left, right] by doing this.
a_eval = f(a)
b_eval = f(b)
# Rare case, f(a) = 0. Return a, no bisection needed.
if a_eval == 0.0
return a
end
# Similarly, if f(b) = 0, then we have already found the root. Return b.
if b_eval == 0.0
return b
end
# Compare the two evaluations and set left and right accordingly.
if a_eval < b_eval
# If both evaluations are negative, or if both are positive, then
# the bisection method will not work. Return NaN.
if b_eval < 0.0 || a_eval > 0.0
return (a - a) / (a - a)
end
# Otherwise, since f(a) < f(b), set left = a and right = b.
left = a
right = b
# In this case the function starts positive and tends to a negative.
else
# Same sanity check as before. We need one evaluation to be
# negative and one to be positive. Abort if both have the same sign.
if a_eval < 0.0 || b_eval > 0.0
return (a - a) / (a - a)
end
# Since f(a) > f(b), set left = b and right = a.
left = b
right = a
end
# Start the bisection method. Compute the midpoint of a and b.
midpoint = 0.5 * (a + b)
# Iteratively apply the bisection method.
for _ = 1:maximum_number_of_iterations
# If f(x) is very small, we are close to a root and can break out
# of this for loop. Check for this.
eval = f(midpoint)
if abs(eval) <= epsilon
break
end
# Apply bisection to get a better approximation for the root. We
# have f(left) < 0 < f(right). If f(midpoint) < 0, replace the
# interval [left, right] with [midpoint, right]. Set left to the
# midpoint and reset the midpoint to be closer to right.
if eval < 0.0
left = midpoint
midpoint = 0.5 * (midpoint + right)
# In the other case, f(midpoint) > 0, we replace right with the
# midpoint, changing [left, right] into [left, midpoint]. We then
# set the midpoint to be closer to left.
else
right = midpoint
midpoint = 0.5 * (left + midpoint)
end
end
# After n iterations, we are no more than |b - a| / 2^n away from the
# root of the function. 1 / 2^n goes to zero very quickly, meaning the
# convergence is very quick.
return midpoint
end
# pi is somewhere between 3 and 4, and it is a root to sine.
a = 3.0
b = 4.0
# Compute pi using bisection. We should get pi = 3.14159..., accurate
# to about 16 decimals.
pi_by_bisection = bisection(sin, a, b)
println("pi = ", pi_by_bisection)
JavaScript
/******************************************************************************
* LICENSE *
******************************************************************************
* This file is part of mitx_mathematics_programming_examples. *
* *
* mitx_mathematics_programming_examples is free software: you can *
* redistribute it and/or modify it under the terms of the GNU General *
* Public License as published by the Free Software Foundation, either *
* version 3 of the License, or (at your option) any later version. *
* *
* mitx_mathematics_programming_examples is distributed in the hope that *
* it will be useful but WITHOUT ANY WARRANTY; without even the implied *
* warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. *
* See the GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with mitx_mathematics_programming_examples. If not, see *
* <https://www.gnu.org/licenses/>. *
******************************************************************************
* Purpose: *
* Calculates the root of a function using bisection. *
******************************************************************************
* Author: Ryan Maguire *
* Date: 2025/05/18 *
******************************************************************************/
/* Computes the root of a function using the bisection method. */
function bisectionMethod(f, a, b) {
/* Tell the algorithm to stop after several iterations to avoid an *
* infinite loop. Double precision numbers have 52 bits in the mantissa, *
* so if |b - a| ~= 1, after 52 iterations of bisection we will get as *
* close as we can to the root. To allow for |b - a| to be larger, halt *
* the algorithm after at most 64 steps. */
const MAXIMUM_NUMBER_OF_ITERATIONS = 64;
/* Getting exact roots is hard using floating-point numbers. Allow a *
* tolerance in our computation. This value is double precision epsilon. */
const EPSILON = 2.220446049250313E-16;
/* The midpoint for the bisection method. This will update as we iterate.*/
let midpoint;
/* Variable for keeping track of the number of iterations performed. */
let iters;
/* We do not require a < b, nor do we require f(a) < f(b). We only need *
* one of these to evaluate to a negative under f and one to evaluate to *
* positive. We will call the negative entry left and positive one right.*/
let left;
let right;
/* Evaluate f at the two endpoints to determine which is positive and *
* which is negative. We transform [a, b] to [left, right] by doing this.*/
const A_EVAL = f(a);
const B_EVAL = f(b);
/* Rare case, f(a) = 0. Return a, no bisection needed. */
if (A_EVAL == 0.0) {
return a;
}
/* Similarly, if f(b) = 0, then we have already found the root. Return b.*/
if (B_EVAL == 0.0) {
return b;
}
/* Compare the two evaluations and set left and right accordingly. */
if (A_EVAL < B_EVAL) {
/* If both evaluations are negative, or if both are positive, then *
* the bisection method will not work. Return NaN. */
if ((B_EVAL < 0.0) || (A_EVAL > 0.0)) {
return (a - a) / (a - a);
}
/* Otherwise, since f(a) < f(b), set left = a and right = b. */
left = a;
right = b;
/* In this case the function starts positive and tends to a negative. */
} else {
/* Same sanity check as before. We need one evaluation to be *
* negative and one to be positive. Abort if both have the same sign.*/
if ((A_EVAL < 0.0) || (B_EVAL > 0.0)) {
return (a - a) / (a - a);
}
/* Since f(a) > f(b), set left = b and right = a. */
left = b;
right = a;
}
/* Start the bisection method. Compute the midpoint of a and b. */
midpoint = 0.5 * (a + b);
/* Iteratively divide the range in half to find the root. */
for (iters = 0; iters < MAXIMUM_NUMBER_OF_ITERATIONS; ++iters) {
/* If f(x) is very small, we are close to a root and can break out *
* of this for loop. Check for this. */
const EVAL = f(midpoint);
if (Math.abs(EVAL) <= EPSILON) {
break;
}
/* Apply bisection to get a better approximation for the root. We *
* have f(left) < 0 < f(right). If f(midpoint) < 0, replace the *
* interval [left, right] with [midpoint, right]. Set left to the *
* midpoint and reset the midpoint to be closer to right. */
if (EVAL < 0.0) {
left = midpoint;
midpoint = 0.5 * (midpoint + right);
/* In the other case, f(midpoint) > 0, we replace right with the *
* midpoint, changing [left, right] into [left, midpoint]. We then *
* set the midpoint to be closer to left. */
} else {
right = midpoint;
midpoint = 0.5 * (left + midpoint);
}
}
/* After n iterations, we are no more than |b - a| / 2^n away from the *
* root of the function. 1 / 2^n goes to zero very quickly, meaning the *
* convergence is very quick. */
return midpoint;
}
/* End of bisectionMethod. */
/* pi is somewhere between 3 and 4, and it is a root to sine. */
const A = 3.0;
const B = 4.0;
/* Compute pi using bisection. We should get pi = 3.14159..., accurate *
* to about 16 decimals. */
const PI = bisectionMethod(Math.sin, A, B);
console.log("pi = " + PI.toFixed(16));
Kotlin
/******************************************************************************
* LICENSE *
******************************************************************************
* This file is part of mitx_mathematics_programming_examples. *
* *
* mitx_mathematics_programming_examples is free software: you can *
* redistribute it and/or modify it under the terms of the GNU General *
* Public License as published by the Free Software Foundation, either *
* version 3 of the License, or (at your option) any later version. *
* *
* mitx_mathematics_programming_examples is distributed in the hope that *
* it will be useful but WITHOUT ANY WARRANTY; without even the implied *
* warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. *
* See the GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with mitx_mathematics_programming_examples. If not, see *
* <https://www.gnu.org/licenses/>. *
******************************************************************************
* Purpose: *
* Calculates the root of a function using bisection. *
******************************************************************************
* Author: Ryan Maguire *
* Date: 2025/05/18 *
******************************************************************************/
import kotlin.math.abs
import kotlin.math.sin
/* Function pointer notation is a little confusing. Create a typedef for it *
* so we do not need to explicitly use it later. */
typealias RealFunc = (Double) -> Double
/* Computes the root of a function using the bisection method. */
fun bisectionMethod(f: RealFunc, a: Double, b: Double): Double {
/* Tell the algorithm to stop after several iterations to avoid an *
* infinite loop. Double precision numbers have 52 bits in the mantissa, *
* so if |b - a| ~= 1, after 52 iterations of bisection we will get as *
* close as we can to the root. To allow for |b - a| to be larger, halt *
* the algorithm after at most 64 steps. */
val maximum_number_of_iterations: Int = 64
/* Getting exact roots is hard using floating-point numbers. Allow a *
* tolerance in our computation. This value is double precision epsilon. */
val epsilon: Double = 2.220446049250313E-16
/* The midpoint for the bisection method. This will update as we iterate.*/
var midpoint: Double
/* We do not require a < b, nor do we require f(a) < f(b). We only need *
* one of these to evaluate to a negative under f and one to evaluate to *
* positive. We will call the negative entry left and positive one right.*/
var left: Double
var right: Double
/* Evaluate f at the two endpoints to determine which is positive and *
* which is negative. We transform [a, b] to [left, right] by doing this.*/
val aEval: Double = f(a)
val bEval: Double = f(b)
/* Rare case, f(a) = 0. Return a, no bisection needed. */
if (aEval == 0.0) {
return a
}
/* Similarly, if f(b) = 0, then we have already found the root. Return b.*/
if (bEval == 0.0) {
return b
}
/* Compare the two evaluations and set left and right accordingly. */
if (aEval < bEval) {
/* If both evaluations are negative, or if both are positive, then *
* the bisection method will not work. Return NaN. */
if ((bEval < 0.0) || (aEval > 0.0)) {
return (a - a) / (a - a)
}
/* Otherwise, since f(a) < f(b), set left = a and right = b. */
left = a
right = b
/* In this case the function starts positive and tends to a negative. */
} else {
/* Same sanity check as before. We need one evaluation to be *
* negative and one to be positive. Abort if both have the same sign.*/
if ((aEval < 0.0) || (bEval > 0.0)) {
return (a - a) / (a - a)
}
/* Since f(a) > f(b), set left = b and right = a. */
left = b
right = a
}
/* Start the bisection method. Compute the midpoint of a and b. */
midpoint = 0.5 * (a + b)
/* Iteratively divide the range in half to find the root. */
for (iters in 0 .. maximum_number_of_iterations) {
/* If f(x) is very small, we are close to a root and can break out *
* of this for loop. Check for this. */
val eval: Double = f(midpoint)
if (abs(eval) <= epsilon) {
break
}
/* Apply bisection to get a better approximation for the root. We *
* have f(left) < 0 < f(right). If f(midpoint) < 0, replace the *
* interval [left, right] with [midpoint, right]. Set left to the *
* midpoint and reset the midpoint to be closer to right. */
if (eval < 0.0) {
left = midpoint
midpoint = 0.5 * (midpoint + right)
/* In the other case, f(midpoint) > 0, we replace right with the *
* midpoint, changing [left, right] into [left, midpoint]. We then *
* set the midpoint to be closer to left. */
} else {
right = midpoint
midpoint = 0.5 * (left + midpoint)
}
}
/* After n iterations, we are no more than |b - a| / 2^n away from the *
* root of the function. 1 / 2^n goes to zero very quickly, meaning the *
* convergence is very quick. */
return midpoint
}
/* End of bisection_method. */
/* Test out our function, compute the square root of 2. */
fun main() {
/* pi is somewhere between 3 and 4, and it is a root to sine. */
val a: Double = 3.0
val b: Double = 4.0
/* Compute pi using bisection. We should get pi = 3.14159..., accurate *
* to about 16 decimals. */
val pi: Double = bisectionMethod(::sin, a, b)
println("pi = $pi")
}
MATLAB
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LICENSE %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of mitx_mathematics_programming_examples. %
% %
% mitx_mathematics_programming_examples is free software: you can %
% redistribute it and/or modify it under the terms of the GNU General Public %
% License as published by the Free Software Foundation, either version 3 of %
% the License, or (at your option) any later version. %
% %
% mitx_mathematics_programming_examples is distributed in the hope that it %
% will be useful, but WITHOUT ANY WARRANTY; without even the implied %
% warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the %
% GNU General Public License for more details. %
% %
% You should have received a copy of the GNU General Public License %
% along with mitx_mathematics_programming_examples. If not, see %
% <https://www.gnu.org/licenses/>. %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Purpose: %
% Calculates roots of a function using bisection. %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Author: Ryan Maguire %
% Date: March 29, 2025. %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function: %
% bisection %
% Purpose: %
% Computes the root of a function f between a and b. %
% Arguments: %
% f (function): %
% A function f: R -> R. The root of f is computed. %
% a (real): %
% One of the endpoints of the interval for f. %
% b (real): %
% The other endpoint for f. %
% Output: %
% root (real): %
% A root of the function f between a and b. %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function root = bisection(f, a, b)
% Tell the algorithm to stop after several iterations to avoid an
% infinite loop. Double precision numbers have 52 bits in the mantissa,
% so if |b - a| ~= 1, after 52 iterations of bisection we will get as
% close as we can to the root. To allow for |b - a| to be larger, halt
% the algorithm after at most 64 steps.
maximum_number_of_iterations = 64;
% The maximum allowed error. This is double precision epsilon.
epsilon = 2.220446049250313E-16;
% Evaluate f at the two endpoints to determine which is positive and
% which is negative. We transform [a, b] to [left, right] by doing this.
a_eval = f(a);
b_eval = f(b);
% Rare case, f(a) = 0. Return a, no bisection needed.
if a_eval == 0.0
root = a;
return;
end
% Similarly, if f(b) = 0, then we have already found the root. Return b.
if b_eval == 0.0
root = b;
return;
end
% Compare the two evaluations and set left and right accordingly.
if a_eval < b_eval
% If both evaluations are negative, or if both are positive, then
% the bisection method will not work. Return NaN.
if b_eval < 0.0 || a_eval > 0.0
root = (a - a) / (a - a);
return;
end
% Otherwise, since f(a) < f(b), set left = a and right = b.
left = a;
right = b;
% In this case the function starts positive and tends to a negative.
else
% Same sanity check as before. We need one evaluation to be
% negative and one to be positive. Abort if both have the same sign.
if a_eval < 0.0 || b_eval > 0.0
root = (a - a) / (a - a);
return;
end
% Since f(a) > f(b), set left = b and right = a.
left = b;
right = a;
end
% Start the bisection method. Compute the midpoint of a and b.
midpoint = 0.5 * (a + b);
% Iteratively apply the bisection method.
for _ = 1:maximum_number_of_iterations
% If f(x) is very small, we are close to a root and can break out
% of this for loop. Check for this.
f_eval = f(midpoint);
if abs(f_eval) <= epsilon
break;
end
% Apply bisection to get a better approximation for the root. We
% have f(left) < 0 < f(right). If f(midpoint) < 0, replace the
% interval [left, right] with [midpoint, right]. Set left to the
% midpoint and reset the midpoint to be closer to right.
if f_eval < 0.0
left = midpoint;
midpoint = 0.5 * (midpoint + right);
% In the other case, f(midpoint) > 0, we replace right with the
% midpoint, changing [left, right] into [left, midpoint]. We then
% set the midpoint to be closer to left.
else
right = midpoint;
midpoint = 0.5 * (left + midpoint);
end
end
% After n iterations, we are no more than |b - a| / 2^n away from the
% root of the function. 1 / 2^n goes to zero very quickly, meaning the
% convergence is very quick.
root = midpoint;
end
% pi is somewhere between 3 and 4, and it is a root to sine.
a = 3.0;
b = 4.0;
% Compute pi using bisection. We should get pi = 3.14159..., accurate
% to about 16 decimals.
pi_by_bisection = bisection(@sin, a, b);
printf("pi = %.16f\n", pi_by_bisection);
Objective-C
/******************************************************************************
* LICENSE *
******************************************************************************
* This file is part of mitx_mathematics_programming_examples. *
* *
* mitx_mathematics_programming_examples is free software: you can *
* redistribute it and/or modify it under the terms of the GNU General *
* Public License as published by the Free Software Foundation, either *
* version 3 of the License, or (at your option) any later version. *
* *
* mitx_mathematics_programming_examples is distributed in the hope that *
* it will be useful but WITHOUT ANY WARRANTY; without even the implied *
* warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. *
* See the GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with mitx_mathematics_programming_examples. If not, see *
* <https://www.gnu.org/licenses/>. *
******************************************************************************
* Purpose: *
* Calculates the root of a function using the bisection method. *
******************************************************************************
* Author: Ryan Maguire *
* Date: 2025/03/29 *
******************************************************************************/
/* NSObject, base class for all Objective-C classes, found here. math.h and *
* stdio.h are also included in this header file. */
#import <Foundation/Foundation.h>
/* Function pointer notation is a little confusing. Create a typedef for it *
* so we do not need to explicitly use it later. */
typedef double (*function)(double);
/* Class providing an implementation of the bisection method. */
@interface Bisection: NSObject
{
/* How many iterations of the Bisection method we do before exiting. */
const unsigned int maximumNumberOfIterations;
/* Allowed tolerance for the root-finding method. */
const double epsilon;
}
/* Given a function f: R -> R and two real numbers a and b, this finds *
* a root of f between a and b using the bisection method. */
+ (double) root: (function)f start: (double)a finish: (double)b;
@end
@implementation Bisection
/* The error after n iterations is |b - a| / 2^n. Since double has a *
* 52-bit mantissa, if |b - a| ~= 1, then after 52 steps we can halt the *
* program. To allow for |b - a| to be larger, we stop the process after *
* at most 64 iterations. */
static const unsigned int maximumNumberOfIterations = 64U;
/* The maximum allowed error. This is double precision epsilon. */
static const double epsilon = 2.220446049250313E-16;
/* Computes the root of a function using the bisection method. */
+ (double) root: (function)f start: (double)a finish: (double)b
{
/* Variable for keeping track of the number of iterations. */
unsigned int iters;
/* The midpoint for the bisection. This updates as we iterate. */
double midpoint;
/* We do not require a < b, nor f(a) < f(b). We need one of *
* these to evaluate negative under f and one to evaluate to *
* positive. Call the negative entry left and positive one right. */
double left, right;
/* Evaluate f at the endpoints to determine which is positive *
* and which is negative, transforming [a, b] to [left, right]. */
const double aEval = f(a);
const double bEval = f(b);
/* Rare case, f(a) = 0. Return a, no bisection needed. */
if (aEval == 0.0)
return a;
/* Similarly, if f(b) = 0, then we found the root. Return b. */
if (bEval == 0.0)
return b;
/* Compare the two evaluations and set the left and right values. */
if (aEval < bEval)
{
/* If both evaluations are negative, or if both are positive, *
* then the bisection method will not work. Return NaN. */
if (bEval < 0.0 || aEval > 0.0)
return (a - a) / (a - a);
/* Otherwise, since f(a) < f(b), set left = a and right = b. */
left = a;
right = b;
}
/* In this case the function starts positive and goes negative. */
else
{
/* Same sanity check as before. We need one evaluation to be *
* negative and one to be positive. Abort if the signs agree. */
if (aEval < 0.0 || bEval > 0.0)
return (a - a) / (a - a);
/* Since f(a) > f(b), set left = b and right = a. */
left = b;
right = a;
}
/* Start the bisection method. Compute the midpoint of a and b. */
midpoint = 0.5 * (a + b);
/* Iteratively divide the range in half to find the root. */
for (iters = 0; iters < maximumNumberOfIterations; ++iters)
{
/* If f(x) is very small, we are close to a root and can *
* break out of this for loop. Check for this. */
const double eval = f(midpoint);
if (fabs(eval) <= epsilon)
break;
/* Apply bisection to get a better approximation. We have *
* f(left) < 0 < f(right). If f(midpoint) < 0, replace the *
* interval [left, right] with [midpoint, right]. Set left *
* to the midpoint and set midpoint to be closer to right. */
if (eval < 0.0)
{
left = midpoint;
midpoint = 0.5 * (midpoint + right);
}
/* If f(midpoint) > 0, then replace right with the midpoint, *
* changing [left, right] into [left, midpoint]. We then set *
* the midpoint to be closer to left. */
else
{
right = midpoint;
midpoint = 0.5 * (left + midpoint);
}
}
/* After n iterations, we are at most |b - a| / 2^n from the *
* root of the function. 1 / 2^n goes to zero very quickly, *
* meaning the convergence is very quick. */
return midpoint;
}
/* End of root. */
@end
/* End of Bisection implementation. */
/* Main routine used for testing our implementation of the Bisection method. */
int main(void)
{
/* pi is somewhere between 3 and 4, and it is a root to sine. */
const double a = 3.0;
const double b = 4.0;
/* Compute pi using bisection. We should get pi = 3.14159..., accurate *
* to about 16 decimals. */
const double pi = [Bisection root: sin start: a finish: b];
printf("pi = %.16f\n", pi);
return 0;
}
Pascal
(******************************************************************************
* LICENSE *
******************************************************************************
* This file is part of mitx_mathematics_programming_examples. *
* *
* mitx_mathematics_programming_examples is free software: you can *
* redistribute it and/or modify it under the terms of the GNU General Public*
* License as published by the Free Software Foundation, either version 3 of *
* the License, or (at your option) any later version. *
* *
* mitx_mathematics_programming_examples is distributed in the hope that it *
* will be useful but WITHOUT ANY WARRANTY; without even the implied warranty*
* of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with mitx_mathematics_programming_examples. If not, see *
* <https://www.gnu.org/licenses/>. *
******************************************************************************
* Purpose: *
* Computes roots using the bisection method. *
******************************************************************************
* Author: Ryan Maguire *
* Date: 2025/04/17 *
******************************************************************************)
PROGRAM Bisection;
(* We define two constants: The allowed tolerance (or epsilon value), *
* and the maximum number of iterations we allow for the bisection method. *)
CONST
(* Tell the algorithm to stop after several iterations to avoid an *
* infinite loop. Double precision numbers have 52 bits in the mantissa, *
* so if |b - a| ~= 1, after 52 iterations of bisection we will get as *
* close as we can to the root. To allow for |b - a| to be larger, halt *
* the algorithm after at most 64 steps. *)
MaximumNumberOfIterations: Integer = 64;
(* Maximum allowed error. This is double precision epsilon. *)
Epsilon: Real = 2.220446049250313E-16;
(* Pi is between 2 and 4 and is a root to sine. *)
LeftValue: Real = 2.0;
RightValue: Real = 4.0;
(* Create an alias for real functions f: R -> R. This makes the syntax of *
* the parameters for the BisectionMethod function a little more readable. *)
TYPE
RealFunc = Function(Const X: Real): Real;
(* The main program only has one variable, pi computed via bisection. *)
VAR
PiByBisection: Real;
(* Sin is a built-in procedure. Create a function Real -> Real that computes *
* sin(x) so that we may pass it as a parameter to BisectionMethod. *)
Function Sine(const X: Real) : Real;
BEGIN
Sine := SIN(X);
END;
(******************************************************************************
* Function: *
* BisectionMethod *
* Purpose: *
* Computes roots using the bisection method. *
* Arguments: *
* F (RealFunc): *
* A function from [A, B] to the real numbers. *
* A (Real): *
* One of the endpoints for the interval. *
* B (Real): *
* The other the endpoint for the interval. *
* OUTPUT: *
* Root (Real): *
* A root for F between A and B. *
******************************************************************************)
Function BisectionMethod(Const F: RealFunc; Const A, B: Real) : Real;
VAR
(* We do not require A < B, nor F(A) < F(B). We will use Left and Right *
* to re-orient the interval so that F(Left) < F(Right). Midpoint will *
* be the center of the interval, and these three values will be updated *
* iteratively as we perform the bisection method. *)
Left, Right, Midpoint: Real;
(* Variables for the evaluation of F at A, B, and Midpoint, respectively.*)
AEval, BEval, Eval: Real;
(* Dummy variable for tracking how many iterations we've performed. *)
Iters: Integer;
LABEL
(* Several spots allow for early returns in the function. We use a *
* single label for GOTO to allow us to break out of the function. *)
Finished;
BEGIN
(* Initial setup, find out which evaluation is positive and which is *
* negative. If the signs of the two agree we treat this as an error. *)
AEval := F(A);
BEval := F(B);
(* Special case, A is a root. Set the output to A and return. *)
IF (ABS(AEval) < Epsilon) THEN
BEGIN
BisectionMethod := A;
GOTO Finished;
END;
(* Similarly for B, if F(B) ~= 0, return B. *)
IF (ABS(BEval) < Epsilon) THEN
BEGIN
BisectionMethod := B;
GOTO Finished;
END;
(* If F(A) < 0 < F(B), then Left = A and Right = B. *)
IF (AEval < BEval) THEN
BEGIN
(* If F(A) and F(B) have the same sign (both positive or
* both negative), return NaN. Bisection is undefined. *)
IF ((AEval > 0) OR (BEval < 0)) THEN
BEGIN
BisectionMethod := (A - A) / (A - A);
GOTO Finished;
END;
Left := A;
Right := B;
(* If F(B) < 0 < F(A), then Left = B and Right = A. *)
END
ELSE
BEGIN
(* Same sanity check as before, make sure the signs are different. *)
IF ((AEval < 0) OR (BEval > 0)) THEN
BEGIN
BisectionMethod := (A - A) / (A - A);
GOTO Finished;
END;
Left := B;
Right := A;
END;
Midpoint := 0.5 * (A + B);
(* Iteratively perform the bisection method. *)
FOR Iters := 0 TO MaximumNumberOfIterations - 1 DO
BEGIN
(* The interval is cut in half based on the sign of F(Midpoint). *
* Compute this and compare. *)
Eval := F(Midpoint);
(* If MIDPOINT is close to a root we can exit the function. *)
IF (ABS(Eval) < Epsilon) THEN BREAK;
(* Otherwise, divide the range in half. If Eval is negative we *
* replace Left with Midpoint and move Midpoint closer to Right. *)
IF (Eval < 0) THEN BEGIN
Left := Midpoint;
Midpoint := 0.5 * (Midpoint + Right);
(* Likewise, if Eval is positive, replace Right with Midpoint and *
* move Midpoint closer to Left. *)
END ELSE BEGIN
Right := Midpoint;
Midpoint := 0.5 * (Left + Midpoint);
END;
END;
(* Provided |B - A| is not too big, we should now have a very good *
* approximation to the root. *)
BisectionMethod := Midpoint;
Finished:
END;
(* Program for testing our implementation of the bisection method. *)
BEGIN
(* Pi is between 2 and 4 and is a root to sine. *)
PiByBisection := BisectionMethod(@Sine, LeftValue, RightValue);
WriteLn('pi = ', PiByBisection:0:16);
END.
IDL
;------------------------------------------------------------------------------;
; LICENSE ;
;------------------------------------------------------------------------------;
; This file is part of mitx_mathematics_programming_examples. ;
; ;
; mitx_mathematics_programming_examples is free software: you can ;
; redistribute it and/or modify it under the terms of the GNU General Public ;
; License as published by the Free Software Foundation, either version 3 of ;
; the License, or (at your option) any later version. ;
; ;
; mitx_mathematics_programming_examples is distributed in the hope that it ;
; will be useful but WITHOUT ANY WARRANTY; without even the implied warranty ;
; of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ;
; GNU General Public License for more details. ;
; ;
; You should have received a copy of the GNU General Public License ;
; along with mitx_mathematics_programming_examples. If not, see ;
; <https://www.gnu.org/licenses/>. ;
;------------------------------------------------------------------------------;
; Purpose: ;
; Calculates roots using the bisection method. ;
;------------------------------------------------------------------------------;
; Author: Ryan Maguire ;
; Date: 2025/03/09 ;
;------------------------------------------------------------------------------;
;------------------------------------------------------------------------------;
; Function: ;
; BISECTION_METHOD ;
; Purpose: ;
; Computes roots using the bisection method. ;
; Arguments: ;
; FUNC (FUNCTION): ;
; A function from [A, B] to the real numbers. ;
; A (REAL): ;
; One of the endpoints for the interval. ;
; B (REAL): ;
; The other the endpoint for the interval. ;
; OUTPUT: ;
; ROOT (REAL): ;
; A root for FUNC. ;
;------------------------------------------------------------------------------;
FUNCTION BISECTION_METHOD, FUNC, A, B
; Tells the compiler that integers should be 32 bits, not 16.
COMPILE_OPT IDL2
; Error checking code.
ON_ERROR, 2
; Given a continuous function defined on [a, b], the nth iteration of the
; bisection method is at most |b - a| / 2^n away from the root. A double
; precision number has 52 bits in the mantissa, meaning if |b - a| ~= 1,
; then after 52 iterations we are as close to the root as we can get.
; To allow for a larger range, halt the algorithm after 64 iterations.
MAXIMUM_NUMBER_OF_ITERATIONS = 64
; Maximum allowed error. This is double precision epsilon.
EPSILON = 2.220446049250313E-16
; Initial setup, find out which evaluation is positive and which is
; negative. If the signs of the two agree we treat this as an error.
A_EVAL = CALL_FUNCTION(FUNC, A)
B_EVAL = CALL_FUNCTION(FUNC, B)
; Special case, A is a root. Set the output to A and return.
IF (ABS(A_EVAL) LE EPSILON) THEN BEGIN
RETURN, A
ENDIF
; Similarly for B, if FUNC(B) ~= 0, return B.
IF (ABS(B_EVAL) LE EPSILON) THEN BEGIN
RETURN, B
ENDIF
; If FUNC(A) < 0 < FUNC(B), then LEFT = A and RIGHT = B.
IF (A_EVAL LT B_EVAL) THEN BEGIN
; If FUNC(A) and FUNC(B) have the same sign (both positive or
; both negative), return NaN. Bisection is undefined.
IF ((A_EVAL GT 0) OR (B_EVAL LT 0)) THEN BEGIN
RETURN, (A - A) / (A - A)
ENDIF
LEFT = A
RIGHT = B
; If FUNC(B) < 0 < FUNC(A), then LEFT = B and RIGHT = A.
ENDIF ELSE BEGIN
; Same sanity check as before, make sure the signs are different.
IF ((A_EVAL LT 0) OR (B_EVAL GT 0)) THEN BEGIN
RETURN, (A - A) / (A - A)
ENDIF
LEFT = B
RIGHT = A
ENDELSE
MIDPOINT = 0.5 * (A + B)
; Iteratively perform the bisection method.
FOR ITERS = 0, MAXIMUM_NUMBER_OF_ITERATIONS DO BEGIN
; The interval is cut in half based on the sign of FUNC(MIDPOINT).
; Compute this and compare.
EVAL = CALL_FUNCTION(FUNC, MIDPOINT)
; If MIDPOINT is close to a root we can exit the function.
IF (ABS(EVAL) LE EPSILON) THEN BEGIN
BREAK
ENDIF
; Otherwise, divide the range in half. If EVAL is negative we
; replace LEFT with MIDPOINT and move MIDPOINT closer to RIGHT.
IF (EVAL LT 0) THEN BEGIN
LEFT = MIDPOINT
MIDPOINT = 0.5 * (MIDPOINT + RIGHT)
; Likewise, if EVAL is positive, replace RIGHT with MIDPOINT and
; move MIDPOINT closer to LEFT.
ENDIF ELSE BEGIN
RIGHT = MIDPOINT
MIDPOINT = 0.5 * (LEFT + MIDPOINT)
ENDELSE
END
; Provided |B - A| is not too big, we should now have a very good
; approximation to the root.
RETURN, MIDPOINT
END
; Program for testing our implementation of the bisection method.
PRO MAIN
; Tells the compiler that integers should be 32 bits, not 16.
COMPILE_OPT IDL2
; Pi is somewhere between 3 and 4, use this interval.
A = DOUBLE(3.0)
B = DOUBLE(4.0)
; Run the routine, compute pi, and print the result.
PI_BY_BISECTION = BISECTION_METHOD("SIN", A, B)
PRINT, PI_BY_BISECTION, FORMAT = 'PI = %18.16f'
END
Python
"""
################################################################################
# LICENSE #
################################################################################
# This file is part of mitx_mathematics_programming_examples. #
# #
# mitx_mathematics_programming_examples is free software: you can #
# redistribute it and/or modify it under the terms of the GNU General Public #
# License as published by the Free Software Foundation, either version 3 of #
# the License, or (at your option) any later version. #
# #
# mitx_mathematics_programming_examples is distributed in the hope that it #
# will be useful, but WITHOUT ANY WARRANTY; without even the implied #
# warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #
# GNU General Public License for more details. #
# #
# You should have received a copy of the GNU General Public License #
# along with mitx_mathematics_programming_examples. If not, see #
# <https://www.gnu.org/licenses/>. #
################################################################################
# Purpose: #
# Calculates roots of a function using bisection. #
################################################################################
# Author: Ryan Maguire #
# Date: April 17, 2025. #
################################################################################
"""
# Pylint doesn't like "x" as a variable name. Disable this warning.
# pylint: disable = invalid-name
# sin provided here.
import math
# Computes the root of a function via bisection.
def bisection(f, a, b):
"""
Function:
bisection
Purpose:
Computes the root of a function f between a and b.
Arguments:
f (function):
A function f: R -> R. The root of f is computed.
a (float):
One of the endpoints of the interval for f.
b (float):
The other endpoint for f.
Output:
root (float):
A root of the function f between a and b.
"""
# Tell the algorithm to stop after several iterations to avoid an
# infinite loop. Double precision numbers have 52 bits in the mantissa,
# so if |b - a| ~= 1, after 52 iterations of bisection we will get as
# close as we can to the root. To allow for |b - a| to be larger, halt
# the algorithm after at most 64 steps.
maximum_number_of_iterations = 64
# The maximum allowed error. This is double precision epsilon.
epsilon = 2.220446049250313E-16
# Evaluate f at the two endpoints to determine which is positive and
# which is negative. We transform [a, b] to [left, right] by doing this.
a_eval = f(a)
b_eval = f(b)
# Rare case, f(a) = 0. Return a, no bisection needed.
if a_eval == 0.0:
return a
# Similarly, if f(b) = 0, then we have already found the root. Return b.
if b_eval == 0.0:
return b
# Compare the two evaluations and set left and right accordingly.
if a_eval < b_eval:
# If both evaluations are negative, or if both are positive, then
# the bisection method will not work. Return NaN.
if b_eval < 0.0 or a_eval > 0.0:
return (a - a) / (a - a)
# Otherwise, since f(a) < f(b), set left = a and right = b.
left = a
right = b
# In this case the function starts positive and tends to a negative.
else:
# Same sanity check as before. We need one evaluation to be
# negative and one to be positive. Abort if both have the same sign.
if a_eval < 0.0 or b_eval > 0.0:
return (a - a) / (a - a)
# Since f(a) > f(b), set left = b and right = a.
left = b
right = a
# Start the bisection method. Compute the midpoint of a and b.
midpoint = 0.5 * (a + b)
# Iteratively apply the bisection method.
for _ in range(maximum_number_of_iterations):
# If f(x) is very small, we are close to a root and can break out
# of this for loop. Check for this.
f_eval = f(midpoint)
if abs(f_eval) <= epsilon:
break
# Apply bisection to get a better approximation for the root. We
# have f(left) < 0 < f(right). If f(midpoint) < 0, replace the
# interval [left, right] with [midpoint, right]. Set left to the
# midpoint and reset the midpoint to be closer to right.
if f_eval < 0.0:
left = midpoint
midpoint = 0.5 * (midpoint + right)
# In the other case, f(midpoint) > 0, we replace right with the
# midpoint, changing [left, right] into [left, midpoint]. We then
# set the midpoint to be closer to left.
else:
right = midpoint
midpoint = 0.5 * (left + midpoint)
# After n iterations, we are no more than |b - a| / 2^n away from the
# root of the function. 1 / 2^n goes to zero very quickly, meaning the
# convergence is very quick.
return midpoint
# pi is somewhere between 3 and 4, and it is a root to sine.
A = 3.0
B = 4.0
# Compute pi using bisection. We should get pi = 3.14159..., accurate
# to about 16 decimals.
pi_by_bisection = bisection(math.sin, A, B)
print(f'pi = {pi_by_bisection}')
R
################################################################################
# LICENSE #
################################################################################
# This file is part of mitx_mathematics_programming_examples. #
# #
# mitx_mathematics_programming_examples is free software: you can #
# redistribute it and/or modify it under the terms of the GNU General Public #
# License as published by the Free Software Foundation, either version 3 of #
# the License, or (at your option) any later version. #
# #
# mitx_mathematics_programming_examples is distributed in the hope that it #
# will be useful, but WITHOUT ANY WARRANTY; without even the implied #
# warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #
# GNU General Public License for more details. #
# #
# You should have received a copy of the GNU General Public License #
# along with mitx_mathematics_programming_examples. If not, see #
# <https://www.gnu.org/licenses/>. #
################################################################################
# Purpose: #
# Calculates roots of a function using bisection. #
################################################################################
# Author: Ryan Maguire #
# Date: May 18, 2025. #
################################################################################
# Computes the root of a function via bisection.
bisection <- function(f, a, b) {
# Tell the algorithm to stop after several iterations to avoid an
# infinite loop. Double precision numbers have 52 bits in the mantissa,
# so if |b - a| ~= 1, after 52 iterations of bisection we will get as
# close as we can to the root. To allow for |b - a| to be larger, halt
# the algorithm after at most 64 steps.
maximum_number_of_iterations <- 64
# The maximum allowed error. This is double precision epsilon.
epsilon <- 2.220446049250313E-16
# Evaluate f at the two endpoints to determine which is positive and
# which is negative. We transform [a, b] to [left, right] by doing this.
a_eval <- f(a)
b_eval <- f(b)
# Rare case, f(a) = 0. Return a, no bisection needed.
if (a_eval == 0.0) {
return(a)
}
# Similarly, if f(b) = 0, then we have already found the root. Return b.
if (b_eval == 0.0) {
return(b)
}
# Compare the two evaluations and set left and right accordingly.
if (a_eval < b_eval) {
# If both evaluations are negative, or if both are positive, then
# the bisection method will not work. Return NaN.
if ((b_eval < 0.0) || (a_eval > 0.0)) {
return((a - a) / (a - a))
}
# Otherwise, since f(a) < f(b), set left = a and right = b.
left <- a
right <- b
# In this case the function starts positive and tends to a negative.
} else {
# Same sanity check as before. We need one evaluation to be
# negative and one to be positive. Abort if both have the same sign.
if ((a_eval < 0.0) || (b_eval > 0.0)) {
return((a - a) / (a - a))
}
# Since f(a) > f(b), set left = b and right = a.
left <- b
right <- a
}
# Start the bisection method. Compute the midpoint of a and b.
midpoint <- 0.5 * (a + b)
# Iteratively apply the bisection method.
for (iters in 1:maximum_number_of_iterations) {
# If f(x) is very small, we are close to a root and can break out
# of this for loop. Check for this.
f_eval <- f(midpoint)
if (abs(f_eval) <= epsilon) {
break
}
# Apply bisection to get a better approximation for the root. We
# have f(left) < 0 < f(right). If f(midpoint) < 0, replace the
# interval [left, right] with [midpoint, right]. Set left to the
# midpoint and reset the midpoint to be closer to right.
if (f_eval < 0.0) {
left <- midpoint
midpoint <- 0.5 * (midpoint + right)
# In the other case, f(midpoint) > 0, we replace right with the
# midpoint, changing [left, right] into [left, midpoint]. We then
# set the midpoint to be closer to left.
} else {
right <- midpoint
midpoint <- 0.5 * (left + midpoint)
}
}
# After n iterations, we are no more than |b - a| / 2^n away from the
# root of the function. 1 / 2^n goes to zero very quickly, meaning the
# convergence is very quick.
return(midpoint)
}
# pi is somewhere between 3 and 4, and it is a root to sine.
A <- 3.0
B <- 4.0
# Compute pi using bisection. We should get pi = 3.14159..., accurate
# to about 16 decimals.
pi_by_bisection <- bisection(sin, A, B)
output = sprintf("pi = %.16f", pi_by_bisection)
message(output)
Rust
/******************************************************************************
* LICENSE *
******************************************************************************
* This file is part of mitx_mathematics_programming_examples. *
* *
* mitx_mathematics_programming_examples is free software: you can *
* redistribute it and/or modify it under the terms of the GNU General *
* Public License as published by the Free Software Foundation, either *
* version 3 of the License, or (at your option) any later version. *
* *
* mitx_mathematics_programming_examples is distributed in the hope that *
* it will be useful but WITHOUT ANY WARRANTY; without even the implied *
* warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. *
* See the GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with mitx_mathematics_programming_examples. If not, see *
* <https://www.gnu.org/licenses/>. *
******************************************************************************
* Purpose: *
* Calculates the root of a function using bisection. *
******************************************************************************
* Author: Ryan Maguire *
* Date: 2025/04/17 *
******************************************************************************/
/* Function pointer notation is a little confusing. Create a typedef for it *
* so we do not need to explicitly use it later. */
type RealFunc = fn(f64) -> f64;
/* Computes the root of a function using the bisection method. */
fn bisection_method(f: RealFunc, a: f64, b: f64) -> f64 {
/* Tell the algorithm to stop after several iterations to avoid an *
* infinite loop. Double precision numbers have 52 bits in the mantissa, *
* so if |b - a| ~= 1, after 52 iterations of bisection we will get as *
* close as we can to the root. To allow for |b - a| to be larger, halt *
* the algorithm after at most 64 steps. */
const MAXIMUM_NUMBER_OF_ITERATIONS: u32 = 64;
/* Getting exact roots is hard using floating-point numbers. Allow a *
* tolerance in our computation. This value is double precision epsilon. */
const EPSILON: f64 = 2.220446049250313E-16;
/* The midpoint for the bisection method. This will update as we iterate.*/
let mut midpoint: f64;
/* We do not require a < b, nor do we require f(a) < f(b). We only need *
* one of these to evaluate to a negative under f and one to evaluate to *
* positive. We will call the negative entry left and positive one right.*/
let mut left: f64;
let mut right: f64;
/* Evaluate f at the two endpoints to determine which is positive and *
* which is negative. We transform [a, b] to [left, right] by doing this.*/
let a_eval: f64 = f(a);
let b_eval: f64 = f(b);
/* Rare case, f(a) = 0. Return a, no bisection needed. */
if a_eval == 0.0 {
return a;
}
/* Similarly, if f(b) = 0, then we have already found the root. Return b.*/
if b_eval == 0.0 {
return b;
}
/* Compare the two evaluations and set left and right accordingly. */
if a_eval < b_eval {
/* If both evaluations are negative, or if both are positive, then *
* the bisection method will not work. Return NaN. */
if b_eval < 0.0 || a_eval > 0.0 {
return (a - a) / (a - a);
}
/* Otherwise, since f(a) < f(b), set left = a and right = b. */
left = a;
right = b;
/* In this case the function starts positive and tends to a negative. */
} else {
/* Same sanity check as before. We need one evaluation to be *
* negative and one to be positive. Abort if both have the same sign.*/
if a_eval < 0.0 || b_eval > 0.0 {
return (a - a) / (a - a);
}
/* Since f(a) > f(b), set left = b and right = a. */
left = b;
right = a;
}
/* Start the bisection method. Compute the midpoint of a and b. */
midpoint = 0.5 * (a + b);
/* Iteratively divide the range in half to find the root. */
for _ in 0 .. MAXIMUM_NUMBER_OF_ITERATIONS {
/* If f(x) is very small, we are close to a root and can break out *
* of this for loop. Check for this. */
let eval: f64 = f(midpoint);
if eval.abs() <= EPSILON {
break;
}
/* Apply bisection to get a better approximation for the root. We *
* have f(left) < 0 < f(right). If f(midpoint) < 0, replace the *
* interval [left, right] with [midpoint, right]. Set left to the *
* midpoint and reset the midpoint to be closer to right. */
if eval < 0.0 {
left = midpoint;
midpoint = 0.5 * (midpoint + right);
/* In the other case, f(midpoint) > 0, we replace right with the *
* midpoint, changing [left, right] into [left, midpoint]. We then *
* set the midpoint to be closer to left. */
} else {
right = midpoint;
midpoint = 0.5 * (left + midpoint);
}
}
/* After n iterations, we are no more than |b - a| / 2^n away from the *
* root of the function. 1 / 2^n goes to zero very quickly, meaning the *
* convergence is very quick. */
return midpoint;
}
/* End of bisection_method. */
/* Main routine used for testing our implementation of the bisection method. */
fn main() {
/* pi is somewhere between 3 and 4, and it is a root to sine. */
const A: f64 = 3.0;
const B: f64 = 4.0;
/* Compute pi using bisection. We should get pi = 3.14159..., accurate *
* to about 16 decimals. */
let pi: f64 = bisection_method(f64::sin, A, B);
println!("pi( = {}", pi);
}
Swift
/******************************************************************************
* LICENSE *
******************************************************************************
* This file is part of mitx_mathematics_programming_examples. *
* *
* mitx_mathematics_programming_examples is free software: you can *
* redistribute it and/or modify it under the terms of the GNU General *
* Public License as published by the Free Software Foundation, either *
* version 3 of the License, or (at your option) any later version. *
* *
* mitx_mathematics_programming_examples is distributed in the hope that *
* it will be useful but WITHOUT ANY WARRANTY; without even the implied *
* warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. *
* See the GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with mitx_mathematics_programming_examples. If not, see *
* <https://www.gnu.org/licenses/>. *
******************************************************************************
* Purpose: *
* Calculates the root of a function using bisection. *
******************************************************************************
* Author: Ryan Maguire *
* Date: 2025/04/17 *
******************************************************************************/
/* fabs (floating-point absolute value) found here. */
import Foundation
/* Function pointer notation is a little confusing. Create a typedef for it *
* so we do not need to explicitly use it later. */
typealias RealFunc = (Double) -> Double
/* Computes the root of a function using the bisection method. */
func bisectionMethod(f: RealFunc, a: Double, b: Double) -> Double {
/* Tell the algorithm to stop after several iterations to avoid an *
* infinite loop. Double precision numbers have 52 bits in the mantissa, *
* so if |b - a| ~= 1, after 52 iterations of bisection we will get as *
* close as we can to the root. To allow for |b - a| to be larger, halt *
* the algorithm after at most 64 steps. */
let maximum_number_of_iterations: UInt32 = 64
/* Getting exact roots is hard using floating-point numbers. Allow a *
* tolerance in our computation. This value is double precision epsilon. */
let epsilon: Double = 2.220446049250313E-16
/* The midpoint for the bisection method. This will update as we iterate.*/
var midpoint: Double
/* We do not require a < b, nor do we require f(a) < f(b). We only need *
* one of these to evaluate to a negative under f and one to evaluate to *
* positive. We will call the negative entry left and positive one right.*/
var left: Double
var right: Double
/* Evaluate f at the two endpoints to determine which is positive and *
* which is negative. We transform [a, b] to [left, right] by doing this.*/
let a_eval: Double = f(a)
let b_eval: Double = f(b)
/* Rare case, f(a) = 0. Return a, no bisection needed. */
if a_eval == 0.0 {
return a
}
/* Similarly, if f(b) = 0, then we have already found the root. Return b.*/
if b_eval == 0.0 {
return b
}
/* Compare the two evaluations and set left and right accordingly. */
if a_eval < b_eval {
/* If both evaluations are negative, or if both are positive, then *
* the bisection method will not work. Return NaN. */
if b_eval < 0.0 || a_eval > 0.0 {
return (a - a) / (a - a)
}
/* Otherwise, since f(a) < f(b), set left = a and right = b. */
left = a
right = b
/* In this case the function starts positive and tends to a negative. */
} else {
/* Same sanity check as before. We need one evaluation to be *
* negative and one to be positive. Abort if both have the same sign.*/
if a_eval < 0.0 || b_eval > 0.0 {
return (a - a) / (a - a)
}
/* Since f(a) > f(b), set left = b and right = a. */
left = b
right = a
}
/* Start the bisection method. Compute the midpoint of a and b. */
midpoint = 0.5 * (a + b)
/* Iteratively divide the range in half to find the root. */
for _ in 0 ..< maximum_number_of_iterations {
/* If f(x) is very small, we are close to a root and can break out *
* of this for loop. Check for this. */
let eval: Double = f(midpoint)
if abs(eval) <= epsilon {
break
}
/* Apply bisection to get a better approximation for the root. We *
* have f(left) < 0 < f(right). If f(midpoint) < 0, replace the *
* interval [left, right] with [midpoint, right]. Set left to the *
* midpoint and reset the midpoint to be closer to right. */
if eval < 0.0 {
left = midpoint
midpoint = 0.5 * (midpoint + right)
/* In the other case, f(midpoint) > 0, we replace right with the *
* midpoint, changing [left, right] into [left, midpoint]. We then *
* set the midpoint to be closer to left. */
} else {
right = midpoint
midpoint = 0.5 * (left + midpoint)
}
}
/* After n iterations, we are no more than |b - a| / 2^n away from the *
* root of the function. 1 / 2^n goes to zero very quickly, meaning the *
* convergence is very quick. */
return midpoint
}
/* End of bisection_method. */
/* pi is somewhere between 3 and 4, and it is a root to sine. */
let a: Double = 3.0
let b: Double = 4.0
/* Compute pi using bisection. We should get pi = 3.14159..., accurate *
* to about 16 decimals. */
let pi: Double = bisectionMethod(f: sin, a: a, b: b)
print("pi = \(pi)")
TypeScript
/******************************************************************************
* LICENSE *
******************************************************************************
* This file is part of mitx_mathematics_programming_examples. *
* *
* mitx_mathematics_programming_examples is free software: you can *
* redistribute it and/or modify it under the terms of the GNU General *
* Public License as published by the Free Software Foundation, either *
* version 3 of the License, or (at your option) any later version. *
* *
* mitx_mathematics_programming_examples is distributed in the hope that *
* it will be useful but WITHOUT ANY WARRANTY; without even the implied *
* warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. *
* See the GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with mitx_mathematics_programming_examples. If not, see *
* <https://www.gnu.org/licenses/>. *
******************************************************************************
* Purpose: *
* Calculates the root of a function using bisection. *
******************************************************************************
* Author: Ryan Maguire *
* Date: 2025/05/18 *
******************************************************************************/
/* Type for a function of the form f: R -> R. */
type realfunc = (x: number) => number;
/* Computes the root of a function using the bisection method. */
function bisectionMethod(f: realfunc, a: number, b: number): number {
/* Tell the algorithm to stop after several iterations to avoid an *
* infinite loop. Double precision numbers have 52 bits in the mantissa, *
* so if |b - a| ~= 1, after 52 iterations of bisection we will get as *
* close as we can to the root. To allow for |b - a| to be larger, halt *
* the algorithm after at most 64 steps. */
const MAXIMUM_NUMBER_OF_ITERATIONS: number = 64;
/* Getting exact roots is hard using floating-point numbers. Allow a *
* tolerance in our computation. This value is double precision epsilon. */
const EPSILON: number = 2.220446049250313E-16;
/* The midpoint for the bisection method. This will update as we iterate.*/
let midpoint: number;
/* Variable for keeping track of the number of iterations performed. */
let iters: number;
/* We do not require a < b, nor do we require f(a) < f(b). We only need *
* one of these to evaluate to a negative under f and one to evaluate to *
* positive. We will call the negative entry left and positive one right.*/
let left: number;
let right: number;
/* Evaluate f at the two endpoints to determine which is positive and *
* which is negative. We transform [a, b] to [left, right] by doing this.*/
const A_EVAL: number = f(a);
const B_EVAL: number = f(b);
/* Rare case, f(a) = 0. Return a, no bisection needed. */
if (A_EVAL == 0.0) {
return a;
}
/* Similarly, if f(b) = 0, then we have already found the root. Return b.*/
if (B_EVAL == 0.0) {
return b;
}
/* Compare the two evaluations and set left and right accordingly. */
if (A_EVAL < B_EVAL) {
/* If both evaluations are negative, or if both are positive, then *
* the bisection method will not work. Return NaN. */
if ((B_EVAL < 0.0) || (A_EVAL > 0.0)) {
return (a - a) / (a - a);
}
/* Otherwise, since f(a) < f(b), set left = a and right = b. */
left = a;
right = b;
/* In this case the function starts positive and tends to a negative. */
} else {
/* Same sanity check as before. We need one evaluation to be *
* negative and one to be positive. Abort if both have the same sign.*/
if ((A_EVAL < 0.0) || (B_EVAL > 0.0)) {
return (a - a) / (a - a);
}
/* Since f(a) > f(b), set left = b and right = a. */
left = b;
right = a;
}
/* Start the bisection method. Compute the midpoint of a and b. */
midpoint = 0.5 * (a + b);
/* Iteratively divide the range in half to find the root. */
for (iters = 0; iters < MAXIMUM_NUMBER_OF_ITERATIONS; ++iters) {
/* If f(x) is very small, we are close to a root and can break out *
* of this for loop. Check for this. */
const EVAL: number = f(midpoint);
if (Math.abs(EVAL) <= EPSILON) {
break;
}
/* Apply bisection to get a better approximation for the root. We *
* have f(left) < 0 < f(right). If f(midpoint) < 0, replace the *
* interval [left, right] with [midpoint, right]. Set left to the *
* midpoint and reset the midpoint to be closer to right. */
if (EVAL < 0.0) {
left = midpoint;
midpoint = 0.5 * (midpoint + right);
/* In the other case, f(midpoint) > 0, we replace right with the *
* midpoint, changing [left, right] into [left, midpoint]. We then *
* set the midpoint to be closer to left. */
} else {
right = midpoint;
midpoint = 0.5 * (left + midpoint);
}
}
/* After n iterations, we are no more than |b - a| / 2^n away from the *
* root of the function. 1 / 2^n goes to zero very quickly, meaning the *
* convergence is very quick. */
return midpoint;
}
/* End of bisectionMethod. */
/* pi is somewhere between 3 and 4, and it is a root to sine. */
const A: number = 3.0;
const B: number = 4.0;
/* Compute pi using bisection. We should get pi = 3.14159..., accurate *
* to about 16 decimals. */
const PI: number = bisectionMethod(Math.sin, A, B);
const OUTPUT: string = "pi = " + PI.toFixed(16);
console.log(OUTPUT);