14회 실험 · 존재와 값

14회 실험 · 존재와 값#

이 실험

짝 서술

14회 · 중간값 정리

확인하는 것

이분법이 좁히는 구간과 sympy·scipy가 주는 값을 대조한다. 조건이 빠지면 결론이 무너지는 두 장면을 그린다

정리 39는 \(f(c)=0\)인 \(c\)가 존재한다고만 말한다. 그 값이 얼마인지는 말하지 않는다. 서술 파트 2절의 흔한 오해 상자가 지적한 것이 이 구분이며, 이 노트북은 그 구분을 코드로 갈라 놓는다.

존재는 정리가 주고 값은 알고리즘이 준다. 이분법을 스무 번 돌리면 근이 든 구간이 백만분의 일 아래로 좁아지지만, 그 절차 어디에도 “근이 있다”는 근거가 없다. 근거는 완비성 공리이다.

그다음 sympy와 scipy로 같은 값을 구해 대조한다. sympy는 실근이 몇 개인지까지 알려 주는데, 그것은 정리가 말하지 않은 별개의 정보이다.

마지막으로 서술 파트 3절의 반례 두 개를 그린다. 각각 증명의 어느 단계에서 막히는지 함께 출력한다.

서술 파트의 내용

여기서 확인하는 방법

정리 39의 가정

1절이 양 끝의 부호가 다름을 먼저 출력한다

존재와 값의 구분

2절이 존재(정리)와 값(알고리즘)을 나란히 놓는다

근이 하나뿐인 것은 아니다

2절이 실근 개수와 차수를 함께 출력한다

3절의 반례 두 개

3절이 두 그림을 그리고 무너지는 단계를 출력한다

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
from scipy.optimize import brentq

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_CEX = "#9467bd"
C_SYM = "#2ca02c"
C_BAND = "#808080"

1. 이분법이 구간을 좁힌다#

서술 파트 2절에서 \(f(x) = x^3 - x - 1\)을 \([1,2]\)에서 보았다. \(f(1) = -1\)이고 \(f(2) = 5\)이므로 정리 39의 가정이 성립한다.

가정부터 확인하고 시작한다. 부호가 다르다는 것이 정리가 요구하는 전부이다.

# Bisection narrows an interval that must contain a zero
def f(t):
    return t ** 3 - t - 1


left, right = 1.0, 2.0

print(f"f({left}) = {f(left)}")
print(f"f({right}) = {f(right)}")
print("signs differ :", f(left) * f(right) < 0)

f(1.0) = -1.0
f(2.0) = 5.0
signs differ : True

이제 구간을 반씩 줄인다. 중점에서의 함숫값 부호를 보고 근이 든 쪽을 남긴다.

# Twenty halvings: the width is exactly halved each time
header = f"{'step':>5} {'left':>14} {'right':>14} {'width':>14}"
print(header)
print("-" * len(header))

for step in range(1, 21):
    mid = 0.5 * (left + right)
    if f(left) * f(mid) <= 0.0:
        right = mid
    else:
        left = mid
    if step <= 3 or step >= 18:
        print(f"{step:>5} {left:>14.8f} {right:>14.8f} {right - left:>14.6e}")

print()
print(f"final interval : [{left:.8f}, {right:.8f}]")
print(f"width          : {right - left:.6e}")

 step           left          right          width
--------------------------------------------------
    1     1.00000000     1.50000000   5.000000e-01
    2     1.25000000     1.50000000   2.500000e-01
    3     1.25000000     1.37500000   1.250000e-01
   18     1.32471466     1.32471848   3.814697e-06
   19     1.32471657     1.32471848   1.907349e-06
   20     1.32471752     1.32471848   9.536743e-07

final interval : [1.32471752, 1.32471848]
width          : 9.536743e-07

폭이 매 단계 정확히 절반이 된다. 스무 단계 뒤에 \(2^{-20} = 9.536743\times10^{-7}\)이다.

그런데 이 절차는 근을 찾는 방법이지 근의 존재를 보장하는 정리가 아니다. 구간이 계속 좁아진다는 것과 그 구간들이 어떤 한 점을 남긴다는 것은 다른 말이다. 뒤쪽을 말하려면 수열의 극한이나 중첩구간에 관한 도구가 있어야 하는데, PART I에는 그런 것이 없다.

존재의 근거는 공리 1이다. 이분법이 하는 일은 이미 존재하는 것의 위치를 좁히는 것뿐이다.

2. 존재와 값은 다른 문제이다#

정리가 준 것과 계산이 준 것을 나란히 놓는다. sympy는 다항식의 실근을 정확히 다룰 수 있다.

# The theorem gives existence; sympy gives the value and the count
t = sp.Symbol("t", real=True)
poly = t ** 3 - t - 1
roots = sp.real_roots(poly)

print("number of real roots :", len(roots))
print("degree of the polynomial :", sp.degree(poly, t))
print()
sp.pprint(roots[0])
print()
print("root to 20 digits :", roots[0].evalf(20))

number of real roots : 1
degree of the polynomial : 3

       ⎛ 3           ⎞
CRootOf⎝x  - x - 1, 0⎠

root to 20 digits : 1.3247179572447460260

실근이 하나인데 차수는 셋이다. 두 수가 다르다는 것에 유의한다. 정리 39는 근이 하나 있다고만 말했고, 하나뿐이라는 것은 정리가 준 정보가 아니다. sympy가 계산으로 알려 준 별개의 사실이다.

scipy로 같은 값을 수치로 구해 대조한다.

# scipy needs the same hypothesis: opposite signs at the two ends
numeric = brentq(lambda x: x ** 3 - x - 1, 1.0, 2.0)
exact = float(roots[0].evalf(20))

print("brentq     :", repr(numeric))
print("sympy      :", repr(exact))
print("difference :", f"{abs(exact - numeric):.3e}")
print()
print("root inside the interval from section 1 :", left <= numeric <= right)

brentq     : 1.324717957244746
sympy      : 1.324717957244746
difference : 0.000e+00

root inside the interval from section 1 : True

세 방법이 같은 값을 가리킨다. 1절의 이분법이 남긴 구간 안에 그 값이 들어 있다.

brentq가 요구하는 조건에 유의한다. 양 끝에서 부호가 달라야 동작하며, 이것은 정리 39의 가정과 같은 조건이다. 수치 알고리즘이 정리의 가정을 그대로 입력으로 받는 셈이다.

3. 조건이 빠지면 무너진다#

서술 파트 3절의 반례 두 개를 그린다. 왼쪽은 연속이 아닌 경우이고 오른쪽은 정의역이 구간이 아닌 경우이다.

# Two counterexamples, each breaking a different step of the proof
fig, axes = plt.subplots(1, 2, figsize=(11.0, 4.0))

ax = axes[0]
xs_left = np.linspace(-1.0, 0.0, 200)[:-1]
xs_right = np.linspace(0.0, 1.0, 200)
ax.plot(xs_left, -np.ones_like(xs_left), color=C_CEX, linewidth=1.8)
ax.plot(xs_right, np.ones_like(xs_right), color=C_CEX, linewidth=1.8)
ax.plot([0.0], [1.0], marker="o", color=C_CEX)
ax.plot([0.0], [-1.0], marker="o", markerfacecolor="white", markeredgecolor=C_CEX)
ax.axhline(0.0, color=C_BAND, linestyle="--", linewidth=1.0)
ax.set_title("f(-1) = -1, f(1) = 1 but f never equals 0")
ax.set_xlabel("x")
ax.set_ylabel("f(x)")
ax.set_ylim(-1.6, 1.6)

ax = axes[1]
piece_a = np.linspace(-2.0, -1.0, 100)
piece_b = np.linspace(1.0, 2.0, 100)
ax.plot(piece_a, piece_a, color=C_CEX, linewidth=1.8)
ax.plot(piece_b, piece_b, color=C_CEX, linewidth=1.8)
ax.axhline(0.0, color=C_BAND, linestyle="--", linewidth=1.0)
ax.set_title("domain [-2,-1] U [1,2] : f(x) = x never equals 0")
ax.set_xlabel("x")
ax.set_ylabel("f(x)")
ax.set_xlim(-2.4, 2.4)

fig.tight_layout()
plt.show()

print("left  : c = sup S = 0, f(c) = 1 > 0 -> case 1 needs sign preservation at c,")
print("        but f is not left-continuous there")
print("right : c = sup S = -1, f(c) = -1 < 0 -> case 2 needs a point y with")
print("        c < y < c + delta inside the domain, and there is none")
../_images/28a66a2881ce22fa2372ea08be234f4ed73984dbb2125539abbea97e4c550308.png 
left  : c = sup S = 0, f(c) = 1 > 0 -> case 1 needs sign preservation at c,
        but f is not left-continuous there
right : c = sup S = -1, f(c) = -1 < 0 -> case 2 needs a point y with
        c < y < c + delta inside the domain, and there is none

왼쪽 그림에서 \(x=0\)의 채운 점이 \(f(0)=1\)이고 빈 점이 왼쪽에서 다가갈 때의 값 \(-1\)이다. 두 점이 떨어져 있다는 것이 좌연속이 아니라는 뜻이며, 그래서 따름정리 8를 \(c=0\)에 적용할 수 없다.

오른쪽 그림에서 두 선분 사이가 비어 있다. \(c = \sup S = -1\)에서 오른쪽으로 조금 나간 점 \(y\)를 골라야 하는데, \(-1\)과 \(1\) 사이에 정의역의 점이 하나도 없다.

같은 정리가 두 곳에서 서로 다른 이유로 막힌다. 왼쪽은 부호 보존을 쓸 수 없어서이고, 오른쪽은 고른 점이 정의역 밖이어서이다.

마치며#

이분법을 스무 번 돌려 근이 든 구간을 \(10^{-6}\) 아래로 좁혔다. 그러나 그 절차는 위치를 좁힐 뿐이고, 근이 있다는 것은 공리 1가 보장한다. 존재와 값은 다른 문제이다.

sympy가 실근이 하나이고 차수가 셋이라고 알려 주었다. 정리 39는 하나 있다고만 말했으므로, 개수는 정리 밖의 정보이다. brentq가 양 끝의 부호가 다를 것을 요구하는 것은 정리의 가정과 같다.

반례 두 개를 그렸다. 하나는 좌연속이 아니어서 부호 보존을 쓸 수 없고, 다른 하나는 고른 점이 정의역 밖이다. 조건마다 무너지는 자리가 다르다.

서술 파트의 내용

이 노트북의 코드

정리 39의 가정

1절 f(left) * f(right) < 0

존재 증명이며 위치를 알려주지 않는다

1절에서 위치를 좁힌 것은 스무 번의 함숫값 계산이다

근이 하나뿐인 것은 아니다

2절 len(roots)와 sp.degree의 대비

값을 구하는 것은 별개의 작업이다

2절 brentq와 sp.real_roots의 대조

3절 반례 — 연속이 아니면

3절 왼쪽 그림의 채운 점과 빈 점

3절 반례 — 정의역이 구간이 아니면

3절 오른쪽 그림의 두 선분 사이

더 해 볼 것#

  1. 1절의 초기 구간을 \([1, 10]\)으로 바꾸어 스무 단계를 돌려 보자. 폭이 얼마가 되는가. \(10^{-6}\) 아래로 내려가려면 몇 단계가 더 필요한가.

  2. \(f(x) = x^3 - 4x\)를 \([-3, 3]\)에서 이분법에 넣어 보자. 서술 파트에서 본 대로 근이 셋인데 절차는 하나만 준다. 어느 근으로 가는가. 초기 구간을 바꾸면 다른 근으로 갈 수 있는가.

  3. brentq에 부호가 같은 두 끝점을 주면 무슨 일이 일어나는가. 그 오류 메시지가 정리 39의 어느 가정에 해당하는지 적어 보자.

짝이 되는 서술

14회 · 중간값 정리