72 lines
2.2 KiB
Python
72 lines
2.2 KiB
Python
import numpy as np
|
|
|
|
|
|
def step(a, b, n):
|
|
return (b - a) / n
|
|
|
|
|
|
# this is basically linspace
|
|
def interpolation_nodes(a, b, n):
|
|
h = step(a, b, n - 1)
|
|
return [a + i * h for i in range(n)]
|
|
|
|
|
|
def composite_trapezoidal(f, a, b, n):
|
|
h = step(a, b, n)
|
|
x = interpolation_nodes(a, b, n + 1)
|
|
return h / 2 * sum([f(x[i]) + f(x[i + 1]) for i in range(n)])
|
|
|
|
|
|
def error_trapezoidal(f, a, b, n, expected):
|
|
return expected - composite_trapezoidal(f, a, b, n)
|
|
|
|
|
|
def format_error(f, a, b, n, expected):
|
|
return f"error for {f} from {a} to {b} with {n} subintervals is {error_trapezoidal(f, a, b, n, expected)}."
|
|
|
|
|
|
def find_subinterval_count(f, a, b, expected, error_threshold):
|
|
n = 1
|
|
e = lambda n: error_trapezoidal(f, a, b, n, expected)
|
|
while np.abs(e(n)) > error_threshold:
|
|
n += 1
|
|
return n
|
|
|
|
|
|
def estimate_convergence_rate(f, a, b, expected):
|
|
m = 20
|
|
errors = [np.log(abs(error_trapezoidal(f, a, b, n, expected))) for n in range(2, m)]
|
|
steps = [np.log(step(a, b, n)) for n in range(2, m)]
|
|
coeffs = np.polyfit(steps, errors, 1) # log-log regression line
|
|
p = coeffs[0]
|
|
C = np.exp(coeffs[1])
|
|
return p, C
|
|
|
|
|
|
def adaptive_trapezoidal(f, a, b, tol):
|
|
T1 = (b - a) / 2 * (f(a) + f(b))
|
|
c = (a + b) / 2
|
|
T2 = (b - a) / 4 * (f(a) + 2 * f(c) + f(b))
|
|
|
|
E1 = 4 / 3 * (T2 - T1)
|
|
|
|
if abs(E1) <= tol:
|
|
E2 = 1 / 3 * (T2 - T1)
|
|
return T2 + E2, 1
|
|
else:
|
|
left_val, left_count = adaptive_trapezoidal(f, a, c, tol * 0.7)
|
|
right_val, right_count = adaptive_trapezoidal(f, c, b, tol * 0.7)
|
|
return left_val + right_val, left_count + right_count
|
|
|
|
|
|
print(format_error(np.exp, 0, 1, 16, np.e - 1))
|
|
print(f"subintervals: {find_subinterval_count(np.exp, 0, 1, np.e - 1, 1e-5)}")
|
|
print(f"exp p: {estimate_convergence_rate(np.exp, 0, 1, np.e - 1)[0]:.5f}")
|
|
print(f"sqrt p: {estimate_convergence_rate(np.sqrt, 0, 1, 2/3)[0]:.5f}")
|
|
print(f"adaptive estimate exp: {adaptive_trapezoidal(np.exp, 0, 1, 1e-5)[0]:.5f}")
|
|
print(f"adaptive estimate sqrt: {adaptive_trapezoidal(np.sqrt, 0, 1, 1e-5)[0]:.5f}")
|
|
print(f"adaptive n subintervals exp: {adaptive_trapezoidal(np.exp, 0, 1, 1e-5)[1]:.0f}")
|
|
print(
|
|
f"adaptive n subintervals sqrt: {adaptive_trapezoidal(np.sqrt, 0, 1, 1e-5)[1]:.0f}"
|
|
)
|