/
Bisection.cs
83 lines (68 loc) · 2.28 KB
/
Bisection.cs
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
/*
* Copyright Lamont Granquist, Sebastien Gaggini and the MechJeb contributors
* SPDX-License-Identifier: LicenseRef-PD-hp OR Unlicense OR CC0-1.0 OR 0BSD OR MIT-0 OR MIT OR LGPL-2.1+
*/
#nullable enable
using System;
using MechJebLib.Utils;
// ReSharper disable CompareOfFloatsByEqualityOperator
namespace MechJebLib.Rootfinding
{
using RootFunc = Func<double, object?, double>;
public static class Bisection
{
/// <summary>
/// Bisection search.
/// </summary>
/// <param name="f">1 dimensional function</param>
/// <param name="a">minimum bounds</param>
/// <param name="b">maximum bounds</param>
/// <param name="o">state object to pass to function</param>
/// <param name="rtol">tolerance</param>
/// <param name="preferLeft"></param>
/// <returns>value for which the function evaluates to zero</returns>
/// <exception cref="ArgumentException">guess does not brack the root</exception>
public static (double c, double fc) Solve(RootFunc f, double a, double b, object? o, double tolX = 0,
bool preferLeft = true)
{
Check.Finite(a);
Check.Finite(b);
Check.NonNegative(tolX);
if (a > b)
{
(a, b) = (b, a);
preferLeft = !preferLeft;
}
double fa = f(a, o);
double fb = f(b, o);
Check.Finite(fa);
Check.Finite(fb);
if (fa == 0)
return (a, fa);
if (fb == 0)
return (b, fb);
if (fa * fb > 0)
throw new ArgumentException("Bisection rootfinding method: guess does not bracket the root");
while (true)
{
double c = (a + b) / 2;
if (c == a || c == b || b - c < tolX)
break;
double fc = f(c, o);
if (fc * fa < 0)
{
b = c;
fb = fc;
}
else
{
a = c;
fa = fc;
}
}
if (preferLeft)
return (a, fa);
return (b, fb);
}
}
}