11회 실험 · ε과 δ

11회 실험 · ε과 δ#

이 실험

짝 서술

11회 · 극한과 연속

확인하는 것

\(\epsilon\)이 먼저 주어지고 \(\delta\)가 뒤에 정해진다는 순서를, \(\delta\)를 실제로 계산하여 확인한다

서술 파트의 정의 18은 \(\epsilon\)을 받아 \(\delta\)를 내놓으라고 요구한다. 이 요구는 코드로 옮기기 좋다. \(\epsilon\)이 함수의 인자가 되고 \(\delta\)가 반환값이 된다. 8회 실험에서 required_n(f, a, b, eps)가 같은 자리에 있었다.

여기서는 \(f(x) = x^2\) 하나만 다룬다. 서술 파트 예제 7가 \(p=3\), \(\epsilon=0.01\)에서 \(\delta = 1/700\)을 골랐는데, 그 값이 최선인지 아닌지를 확인하는 것이 첫 번째 목적이다.

먼저 이분법으로 통하는 \(\delta\)를 직접 찾는다. 그다음 같은 값을 sympy로 정확히 푼다. 수치는 값 하나를 주고 기호 계산은 \(\epsilon\)과 \(p\)의 관계식을 준다. 마지막으로 그 관계식을 \(p\)의 함수로 그려, 서술 파트 2절 흔한 오해가 말한 \(\delta\)가 점에 의존한다는 사실을 눈으로 본다.

서술 파트의 내용

여기서 확인하는 방법

(11.1)의 \(0 < \lvert x-p \rvert\) 조건

표본에서 \(x = p\)를 빼고 검사한다

(11.9) — \(\delta = \min\{1,\ \epsilon/7\}\)이 최선이 아니라는 점

이분법으로 통하는 가장 큰 \(\delta\)를 찾아 비교한다

\(\delta\)가 \(\epsilon\)과 \(p\) 두 값에 의존한다는 것

\(\delta\)를 \(\epsilon\)과 \(p\)의 식으로 푼다

2절 흔한 오해의 \(1/700\)과 \(1/2100\)

식이 주는 \(\delta\)를 \(p\)의 함수로 그린다

0. 준비#

이 회차에서 쓸 것만 불러온다. 색은 상수로 고정한다. 검정은 정확한 값이고 초록은 기호 계산에서 얻은 근사식이다.

Hide code cell source

 
# Setup: imports, plotting defaults, and the colors used in this lecture
import numpy as np
import matplotlib.pyplot as plt
import sympy as sp

plt.rcParams["figure.figsize"] = (11.0, 4.2)
plt.rcParams["figure.dpi"] = 110
plt.rcParams["axes.grid"] = True
plt.rcParams["grid.alpha"] = 0.3
plt.rcParams["font.size"] = 11

C_FUNC = "#000000"
C_SYM = "#2ca02c"

1. 통하는 δ를 이분법으로 찾는다#

정의 18이 요구하는 것을 그대로 코드로 옮긴다. \(\delta\) 하나를 받아 그것이 통하는지 판정하는 함수가 먼저 필요하다.

판정은 \(0 < |x-p| < \delta\)인 모든 \(x\)에서 \(|f(x)-f(p)| < \epsilon\)인지 보는 것이다. 모든 \(x\)를 볼 수는 없으므로 촘촘한 표본으로 대신한다. 여기서 한 줄이 정의를 그대로 옮긴다 — 표본에서 \(x = p\)를 빼는 줄이 (11.1)의 왼쪽 부등식 \(0 < |x-p|\)에 대응한다.

# Does a given delta work? The line that drops x = p is the condition 0 < |x - p|
def f(x):
    return x ** 2


def delta_works(p, eps, delta, n_samples=2001):
    """True if |f(x) - f(p)| < eps for every sampled x with 0 < |x - p| < delta."""
    xs = np.linspace(p - delta, p + delta, n_samples)
    xs = xs[np.abs(xs - p) > 0.0]
    return bool(np.max(np.abs(f(xs) - f(p))) < eps)


print("p = 3, eps = 0.01, delta = 1/700  :", delta_works(3.0, 0.01, 1.0 / 700))
print("p = 3, eps = 0.01, delta = 0.01   :", delta_works(3.0, 0.01, 0.01))

p = 3, eps = 0.01, delta = 1/700  : True
p = 3, eps = 0.01, delta = 0.01   : False

서술 파트가 고른 \(1/700\)은 통하고, 그보다 큰 \(0.01\)은 통하지 않는다. 그렇다면 통하는 것과 통하지 않는 것의 경계가 어디인지 물을 수 있다.

경계는 이분법으로 찾는다. 통하는 \(\delta\)를 하나 확보한 뒤, 통하지 않는 값과의 사이를 반씩 좁힌다.

# The largest delta that still works, found by halving and then bisecting
def find_delta_bisection(p, eps, delta_start=1.0, n_steps=60):
    """Return the largest delta that passes delta_works, up to bisection error."""
    hi = delta_start
    while not delta_works(p, eps, hi):
        hi /= 2.0
    lo = hi
    hi = lo * 2.0
    for _ in range(n_steps):
        mid = 0.5 * (lo + hi)
        if delta_works(p, eps, mid):
            lo = mid
        else:
            hi = mid
    return lo


eps = 0.01
for p in (1.0, 3.0, 5.0):
    print(f"p = {p}: delta = {find_delta_bisection(p, eps):.10f}")

p = 1.0: delta = 0.0049875621
p = 3.0: delta = 0.0016662040
p = 5.0: delta = 0.0009999000

\(p=3\)에서 \(0.0016662040\)이 나온다. 서술 파트가 고른 \(1/700 = 0.0014285714\)보다 크다.

두 값이 다른 것은 문제가 아니다. (11.8)에서 \(|x+3|\)을 \(7\) 아래로 눌렀는데 실제로 \(x\)가 \(3\) 근처이면 그 값은 \(6\)에 가까우므로, 서술 파트의 \(\delta\)는 필요 이상으로 좁게 잡힌 것이다. 정의 18은 통하는 \(\delta\)가 하나 있으면 된다고 할 뿐 가장 큰 것을 요구하지 않는다. 작은 쪽을 골라도 정의는 충족된다.

세 줄을 함께 보면 \(p\)가 커질수록 \(\delta\)가 작아진다. 이 관계를 다음 절에서 식으로 얻는다.

2. 같은 δ를 식으로 푼다#

이분법은 수치 하나를 준다. \(p\)나 \(\epsilon\)을 바꾸면 처음부터 다시 돌려야 한다. 관계식을 얻으면 그럴 필요가 없다.

\(f(x) = x^2\)은 \(p > 0\) 근처에서 증가하므로, \(|x^2 - p^2| < \epsilon\)이 깨지는 곳은 구간의 양 끝이다. 두 끝에서 각각 등식을 풀어 어느 쪽이 먼저 깨지는지 본다.

여기서 sp.limit을 쓰지 않는다. 극한값을 구하는 것이 목적이 아니라 정의가 요구하는 \(\delta\)를 구하는 것이 목적이기 때문이다.

# Solve for the delta that makes the eps band tight at each end
p, eps_sym, delta = sp.symbols("p eps delta", positive=True)

right = sp.solve(sp.Eq((p + delta) ** 2 - p ** 2, eps_sym), delta)
left = sp.solve(sp.Eq(p ** 2 - (p - delta) ** 2, eps_sym), delta)

print("right end :", right)
print("left end  :", left)

at_3 = {p: 3, eps_sym: sp.Rational(1, 100)}
print()
print("right end at p = 3 :", float(right[0].subs(at_3)))
print("left end at p = 3  :", [float(e.subs(at_3)) for e in left])

right end : [-p + sqrt(eps + p**2)]
left end  : [p - sqrt(-eps + p**2), p + sqrt(-eps + p**2)]

right end at p = 3 : 0.0016662039607268763
left end at p = 3  : [0.0016671298870100271, 5.99833287011299]

오른쪽 끝은 \(\sqrt{p^2+\epsilon} - p\) 하나이고, 왼쪽 끝은 두 근을 준다. 왼쪽의 두 근 중 \(5.998\ldots\)는 \(x\)가 \(p\)를 지나 반대편으로 넘어간 해이므로 버리고 작은 쪽 \(0.0016671299\)를 택한다.

두 값을 비교하면 오른쪽이 더 작다. \(x^2\)은 오른쪽으로 갈 때 더 빨리 커지므로 허용 띠를 먼저 벗어난다. 따라서 둘 중 작은 쪽인 오른쪽 값을 \(\delta\)로 택해야 양쪽 모두 안전하다.

# The right end is the binding one, so it is the delta
delta_expr = right[0]
print("delta =", delta_expr)
print("latex :", sp.latex(delta_expr))

series = sp.expand(sp.series(delta_expr, eps_sym, 0, 3).removeO())
print("series in eps :", series)

delta = -p + sqrt(eps + p**2)
latex : - p + \sqrt{eps + p^{2}}
series in eps : -eps**2/(8*p**3) + eps/(2*p)

\(\delta = \sqrt{p^2+\epsilon} - p\)이고, \(\epsilon\)이 작을 때의 급수 전개는

\[ \delta = \frac{\epsilon}{2p} - \frac{\epsilon^2}{8p^3} + \cdots \]

이다. 첫 항이 \(\epsilon/(2p)\)라는 것이 이 절의 결론이다. \(\epsilon\)에 비례하고 \(p\)에 반비례한다.

식이 이분법의 값과 맞는지 대조한다.

# Compare the formula with the bisection search
header = f"{'p':>5} {'delta (formula)':>18} {'delta (bisection)':>19} {'difference':>13}"
print(header)
print("-" * len(header))
for p_val in (1, 3, 5, 10):
    formula = float(delta_expr.subs({p: p_val, eps_sym: sp.Rational(1, 100)}))
    numeric = find_delta_bisection(float(p_val), 0.01)
    print(f"{p_val:>5} {formula:>18.10f} {numeric:>19.10f} {abs(formula - numeric):>13.2e}")

    p    delta (formula)   delta (bisection)    difference
----------------------------------------------------------
    1       0.0049875621        0.0049875621      3.30e-17
    3       0.0016662040        0.0016662040      4.34e-17
    5       0.0009999000        0.0009999000      1.69e-16
   10       0.0004999875        0.0004999875      7.82e-16

차이가 \(10^{-16}\) 수준이다. 이분법이 찾은 경계와 식이 주는 경계가 같은 값이다.

두 방법이 주는 것은 다르다. 이분법은 주어진 \(p\)와 \(\epsilon\) 하나에 대한 수치를 주고, 기호 계산은 모든 \(p\)와 \(\epsilon\)에 대한 관계식을 준다. 관계식이 있어야 다음 절처럼 \(p\)를 움직이며 볼 수 있다.

3. δ는 점에 따라 달라진다#

서술 파트 2절의 흔한 오해 상자가 \(p=3\)에서 \(1/700\)이고 \(p=10\)에서 \(1/2100\)이라고 적었다. 앞 절의 식으로 이 의존성을 \(p\) 전체에 대하여 그린다.

기호식을 수치 함수로 바꾸는 데 sp.lambdify를 쓴다. 이 한 줄이 기호 계산에서 수치 계산으로 내려오는 다리이다.

# Turn the symbolic delta into a numeric function and plot it against p
delta_fn = sp.lambdify((p, eps_sym), delta_expr, "numpy")

eps = 0.01
ps = np.linspace(0.5, 10.0, 60)
exact = delta_fn(ps, eps)
approx = eps / (2 * ps)

fig, ax = plt.subplots()
ax.plot(ps, exact, color=C_FUNC, linewidth=1.8, label="exact delta")
ax.plot(ps, approx, color=C_SYM, linestyle="--", linewidth=1.6,
        label="approximation eps/(2p)")
ax.set_xlabel("point p")
ax.set_ylabel("required delta")
ax.set_title(f"eps = {eps} is fixed; the delta it needs shrinks as p grows")
ax.legend()
fig.tight_layout()
plt.show()

print(f"p = 0.5 : delta = {float(delta_fn(0.5, eps)):.10f}")
print(f"p = 10  : delta = {float(delta_fn(10.0, eps)):.10f}")
print(f"ratio   = {float(delta_fn(0.5, eps)) / float(delta_fn(10.0, eps)):.7f}")
../_images/b1c6f2053e03c9f5d90400fc9da05493b4ba9916299f289759c26e2f62d4b15d.png 
p = 0.5 : delta = 0.0099019514
p = 10  : delta = 0.0004999875
ratio   = 19.8043978

같은 \(\epsilon = 0.01\)을 지키는 데 \(p = 0.5\)에서는 \(0.0099019514\)가 필요하고 \(p = 10\)에서는 \(0.0004999875\)가 필요하다. 약 \(19.8\)배 차이이다.

근사식 \(\epsilon/(2p)\)로 계산하면 \(p\)가 \(20\)배 커졌으므로 정확히 \(20\)배여야 한다. \(19.8\)은 그보다 조금 작은데, 급수 전개의 둘째 항 \(-\epsilon^2/(8p^3)\)이 \(p\)가 작을 때 더 크게 빠지기 때문이다. 그림에서 두 선이 왼쪽 끝에서 아주 조금 벌어지는 것이 그것이다.

허용 오차를 고정해도 필요한 \(\delta\)는 점마다 다르다. 이것은 정의가 허용하는 일이다. (11.2)에서 \(\delta\)를 고르는 자리가 \(\epsilon\) 뒤이고 \(p\)는 그 전에 이미 고정되어 있기 때문이다.

마치며#

정의 18이 요구하는 작업을 두 방법으로 해 보았다. 이분법은 \(\epsilon\)과 \(p\) 하나에 대한 \(\delta\)를 수치로 찾고, 기호 계산은 \(\delta = \sqrt{p^2+\epsilon}-p\)라는 관계식을 준다. 두 결과가 \(10^{-16}\)까지 일치하였다.

서술 파트가 고른 \(1/700\)은 이분법이 찾은 \(0.0016662040\)보다 작았다. (11.8)에서 \(|x+3|\)을 넉넉하게 눌렀기 때문이며, 정의는 통하는 \(\delta\) 하나만 요구하므로 작은 쪽을 골라도 된다.

\(\delta\)를 \(p\)의 함수로 그리면 오른쪽으로 갈수록 낮아진다. \(p=0.5\)와 \(p=10\) 사이에서 약 \(19.8\)배 차이가 났다. \(\epsilon\)을 고정해도 \(\delta\)는 점마다 다르다.

서술 파트의 내용

이 노트북의 코드

(11.1)의 \(0 < \lvert x-p \rvert\) 조건

1절 delta_works의 xs = xs[np.abs(xs - p) > 0.0]

\(\epsilon\)이 먼저이고 \(\delta\)가 나중이라는 순서

1절 find_delta_bisection(p, eps)의 서명 — eps가 인자이고 delta가 반환값이다

(11.9)의 \(\delta\)가 최선이 아니라는 점

1절에서 \(1/700\)과 \(0.0016662040\)의 대조

\(\delta\)가 \(\epsilon\)과 \(p\) 두 값에 의존한다는 것

2절 sp.solve가 준 \(\sqrt{p^2+\epsilon}-p\)

2절 흔한 오해의 \(p=3\)과 \(p=10\) 대조

3절 delta_fn으로 그린 곡선과 \(19.8\)배

더 해 볼 것#

  1. \(f(x) = x^3\)으로 바꾸어 같은 방식으로 \(\delta\)를 식으로 풀어 보자. sp.solve가 주는 근이 몇 개인지, 그중 어느 것을 택해야 하는지 확인하고, 급수 전개의 첫 항이 무엇이 되는지 보자.

  2. \(\epsilon\)을 \(0.01\)에서 \(0.0001\)로 줄여 3절의 그림을 다시 그려 보자. 정확식과 근사식의 차이가 커지는가 작아지는가. 급수 전개의 둘째 항이 \(\epsilon^2\)에 비례한다는 사실과 맞는지 확인해 보자.

  3. 3절의 그림에서 모든 \(p\)에 대하여 통하는 하나의 \(\delta\)를 고를 수 있는가. \(p\)의 범위를 \([0.5, 10]\)으로 제한하면 어떠한가. 범위를 \((0, \infty)\) 전체로 넓히면 어떠한가.

짝이 되는 서술

11회 · 극한과 연속