1b bisection ported
This commit is contained in:
1 parent
6e35861154
commit
90e0d7638b
3 files changed
+127
-42
No files matched your search
@@ -0,0 +1 @@
|
|||||||
|
__pycache__
|
||||||
@@ -0,0 +1,119 @@
|
|||||||
|
"bisection tailored to fit lab 1b"
|
||||||
|
|
||||||
|
from typing import Callable
|
||||||
|
from numbers import Real
|
||||||
|
from math import isclose, fabs
|
||||||
|
|
||||||
|
|
||||||
|
# this is the already existing bis, see git blame
|
||||||
|
def bis(
|
||||||
|
f: Callable[[Real], Real],
|
||||||
|
low: Real,
|
||||||
|
high: Real,
|
||||||
|
*,
|
||||||
|
rel_tol: Real = 1e-9,
|
||||||
|
abs_tol: Real = 0.0,
|
||||||
|
max_loop: int = 100,
|
||||||
|
) -> tuple[Real, int]:
|
||||||
|
"""real
|
||||||
|
|
||||||
|
rel_tol and abs_tol passes to `math.isclose`
|
||||||
|
"""
|
||||||
|
|
||||||
|
def ldaprx(a, b):
|
||||||
|
return isclose(a, b, rel_tol=rel_tol, abs_tol=abs_tol)
|
||||||
|
|
||||||
|
steps = 0
|
||||||
|
fl = f(low)
|
||||||
|
fh = f(high)
|
||||||
|
# raise ValueError(f"at {low =} f(x) is {fl}, at {high =} f(x) is {fh} which has the same sign")
|
||||||
|
|
||||||
|
if ldaprx(0, fl):
|
||||||
|
return (low, 0)
|
||||||
|
if ldaprx(0, fh):
|
||||||
|
return (high, 0)
|
||||||
|
|
||||||
|
mid = (low + high) / 2
|
||||||
|
|
||||||
|
while (steps := steps + 1) and (max_loop := max_loop - 1) + 1 > 0:
|
||||||
|
match f(mid):
|
||||||
|
case x if ldaprx(0, x):
|
||||||
|
return (mid, steps)
|
||||||
|
case x if x < 0:
|
||||||
|
low = mid
|
||||||
|
case x if x > 0:
|
||||||
|
high = mid
|
||||||
|
mid = (low + high) / 2
|
||||||
|
return (mid, steps)
|
||||||
|
|
||||||
|
|
||||||
|
def bis1b(
|
||||||
|
f: Callable[..., Real],
|
||||||
|
bounds: tuple[Real, Real],
|
||||||
|
rel_tol: Real = 1e-6,
|
||||||
|
maxit: int = 50,
|
||||||
|
fa: tuple = None,
|
||||||
|
fkw: dict = None,
|
||||||
|
) -> tuple[Real, Real, Real, int]:
|
||||||
|
"""bisect: root location zeroes
|
||||||
|
|
||||||
|
uses bisection method to find the root of func
|
||||||
|
|
||||||
|
input:
|
||||||
|
func: name of function
|
||||||
|
bounds: lower and upper guesses
|
||||||
|
rel_tol: desired relative error (default = 0.0001%)
|
||||||
|
maxit: maximum allowable iterations (default = 50)
|
||||||
|
fa, fkw: additional parameters used by func
|
||||||
|
output:
|
||||||
|
root: real root
|
||||||
|
fx: function value at root
|
||||||
|
ea: approximate relative error. not in percent because that is bullshit
|
||||||
|
iter: number of iterations
|
||||||
|
|
||||||
|
function should take one real variable x (Real) + farg and fkwarg like so:
|
||||||
|
`f(x, *fa, **fkw)`
|
||||||
|
"""
|
||||||
|
if fa is None:
|
||||||
|
fa = ()
|
||||||
|
if fkw is None:
|
||||||
|
fkw = {}
|
||||||
|
|
||||||
|
def fi(x):
|
||||||
|
return f(x, *fa, **fkw)
|
||||||
|
|
||||||
|
if len(bounds) != 2:
|
||||||
|
raise ValueError("bounds should have two values?")
|
||||||
|
xl, xu = bounds
|
||||||
|
# shit naming competition
|
||||||
|
|
||||||
|
if fi(xl) * fi(xu) > 0: # maybe incorrect
|
||||||
|
raise ValueError("no sign change")
|
||||||
|
|
||||||
|
it = 0
|
||||||
|
xr = xl
|
||||||
|
ea = 1
|
||||||
|
|
||||||
|
while True:
|
||||||
|
xr_old = xr
|
||||||
|
xr = (xl + xu) / 2
|
||||||
|
it += 1
|
||||||
|
|
||||||
|
if xr != 0:
|
||||||
|
ea = fabs((xr - xr_old) / xr) # approx
|
||||||
|
|
||||||
|
cur_iter = fi(xl) * fi(xu)
|
||||||
|
match cur_iter:
|
||||||
|
case x if x > 0:
|
||||||
|
xu = xr
|
||||||
|
case x if x < 0:
|
||||||
|
xl = xr
|
||||||
|
case 0:
|
||||||
|
ea = 0
|
||||||
|
|
||||||
|
if ea <= rel_tol or it >= maxit:
|
||||||
|
break
|
||||||
|
root = xr
|
||||||
|
fx = fi(xr)
|
||||||
|
|
||||||
|
return (root, fx, ea, it)
|
||||||
+7
-42
@@ -192,7 +192,7 @@
|
|||||||
},
|
},
|
||||||
{
|
{
|
||||||
"cell_type": "code",
|
"cell_type": "code",
|
||||||
"execution_count": 18,
|
"execution_count": 3,
|
||||||
"metadata": {},
|
"metadata": {},
|
||||||
"outputs": [],
|
"outputs": [],
|
||||||
"source": [
|
"source": [
|
||||||
@@ -208,7 +208,7 @@
|
|||||||
},
|
},
|
||||||
{
|
{
|
||||||
"cell_type": "code",
|
"cell_type": "code",
|
||||||
"execution_count": 19,
|
"execution_count": 4,
|
||||||
"metadata": {},
|
"metadata": {},
|
||||||
"outputs": [
|
"outputs": [
|
||||||
{
|
{
|
||||||
@@ -217,7 +217,7 @@
|
|||||||
"np.float64(1.5213797068045676)"
|
"np.float64(1.5213797068045676)"
|
||||||
]
|
]
|
||||||
},
|
},
|
||||||
"execution_count": 19,
|
"execution_count": 4,
|
||||||
"metadata": {},
|
"metadata": {},
|
||||||
"output_type": "execute_result"
|
"output_type": "execute_result"
|
||||||
}
|
}
|
||||||
@@ -258,51 +258,16 @@
|
|||||||
},
|
},
|
||||||
{
|
{
|
||||||
"cell_type": "code",
|
"cell_type": "code",
|
||||||
"execution_count": 20,
|
"execution_count": 1,
|
||||||
"metadata": {},
|
"metadata": {},
|
||||||
"outputs": [],
|
"outputs": [],
|
||||||
"source": [
|
"source": [
|
||||||
"from typing import Callable\n",
|
"from bis import bis # to new file"
|
||||||
"from numbers import Real\n",
|
|
||||||
"from math import isclose\n",
|
|
||||||
"\n",
|
|
||||||
"\n",
|
|
||||||
"def bis(\n",
|
|
||||||
" f: Callable[[Real], Real],\n",
|
|
||||||
" low: Real,\n",
|
|
||||||
" high: Real,\n",
|
|
||||||
" *,\n",
|
|
||||||
" rel_tol: Real = 1e-9,\n",
|
|
||||||
" abs_tol: Real = 0,\n",
|
|
||||||
" max_loop: int = 100,\n",
|
|
||||||
") -> tuple[Real, int]:\n",
|
|
||||||
" ldaprx = lambda a, b: isclose(a, b, rel_tol=rel_tol, abs_tol=abs_tol)\n",
|
|
||||||
" steps = 0\n",
|
|
||||||
" fl = f(low)\n",
|
|
||||||
" fh = f(high)\n",
|
|
||||||
" # raise ValueError(f\"at {low =} f(x) is {fl}, at {high =} f(x) is {fh} which has the same sign\")\n",
|
|
||||||
"\n",
|
|
||||||
" if ldaprx(0, fl):\n",
|
|
||||||
" return (low, 0)\n",
|
|
||||||
" elif ldaprx(0, fh):\n",
|
|
||||||
" return (high, 0)\n",
|
|
||||||
" else:\n",
|
|
||||||
" mid = (low + high) / 2\n",
|
|
||||||
" while (steps := steps + 1) and (max_loop := max_loop - 1) + 1 > 0:\n",
|
|
||||||
" match f(mid):\n",
|
|
||||||
" case x if ldaprx(0, x):\n",
|
|
||||||
" return (mid, steps)\n",
|
|
||||||
" case x if x < 0:\n",
|
|
||||||
" low = mid\n",
|
|
||||||
" case x if x > 0:\n",
|
|
||||||
" high = mid\n",
|
|
||||||
" mid = (low + high) / 2\n",
|
|
||||||
" return (mid, steps)"
|
|
||||||
]
|
]
|
||||||
},
|
},
|
||||||
{
|
{
|
||||||
"cell_type": "code",
|
"cell_type": "code",
|
||||||
"execution_count": 21,
|
"execution_count": 5,
|
||||||
"metadata": {},
|
"metadata": {},
|
||||||
"outputs": [
|
"outputs": [
|
||||||
{
|
{
|
||||||
@@ -311,7 +276,7 @@
|
|||||||
"(1.5213797068045676, 52)"
|
"(1.5213797068045676, 52)"
|
||||||
]
|
]
|
||||||
},
|
},
|
||||||
"execution_count": 21,
|
"execution_count": 5,
|
||||||
"metadata": {},
|
"metadata": {},
|
||||||
"output_type": "execute_result"
|
"output_type": "execute_result"
|
||||||
}
|
}
|
||||||
|
|||||||
Reference in new issue
Block a user