20회 실험 · 성질 넷이 하는 일#
이 실험
짝 서술 |
|
확인하는 것 |
네 성질이 실제로 성립하는지 검사한다. 압축이 조여지는 것을 본다. 원을 한 번도 쓰지 않고 \(\pi\)를 구한다 |
이 노트북에는 다른 회차에 없던 제약이 하나 있다.
서술 파트가 삼각함수를 성질로 특징지었을 뿐 만들지는 않았다. 그런데 numpy와 sympy에는 sin과 cos이 이미 들어 있다. 그 구현이 우리가 적은 네 성질을 만족하는지는 확인할 수 있지만, 그런 함수가 존재한다는 것을 코드가 증명해 주지는 않는다.
그래서 이 노트북이 하는 일은 두 가지이다. 하나는 네 성질과 그로부터 끌어낸 결론들이 실제 구현과 어긋나지 않는지 보는 것이고, 다른 하나는 서술 파트의 논증을 그대로 실행하는 것이다.
4절이 그 두 번째이다. 정리 69가 \(\cos\)의 첫 양의 영점을 잡고 정의 45가 그 두 배를 \(\pi\)로 정의했다. 그 절차를 이분법으로 실행하면 \(3.14159\ldots\)가 나온다. 원의 둘레도 호의 길이도 쓰지 않았다.
서술 파트의 내용 |
여기서 확인하는 방법 |
|---|---|
정의 44의 (T1)~(T3) |
1절이 sympy로 항등식이 \(0\)인지 본다 |
(T4)의 \(c\) |
1절이 부등식이 무너지는 자리를 찾는다 |
2절이 세 값을 한 표에 놓고 조여지는 것을 본다 |
|
2절이 \(1-\cos x\)와 \(x^2\)을 견준다 |
|
3절이 차분몫과 \(\cos\)을 대조한다 |
|
4절이 이분법으로 \(\pi\)를 구한다 |
|
4절이 그 증명의 \(X\)를 실제로 계산한다 |
|
5절이 정의로 구한 적분과 대조한다 |
0. 준비#
코드는 존재를 증명하지 않는다
아래에서 np.sin과 sp.sin을 쓴다. 이것들은 네 성질을 만족하는 함수쌍의 한 구현이다.
정의 44이 “이런 성질을 만족하는 함수쌍이 있다고 하자”로 시작했고, 서술 파트는 그 존재를 증명하지 않았다. 이 노트북도 증명하지 않는다.
코드가 할 수 있는 일은 어긋남을 찾는 것이다. 네 성질 가운데 하나라도 구현이 어기면 그것이 출력에 드러난다. 어기지 않는다고 해서 존재가 증명되는 것은 아니다.
1. 네 성질을 검사한다#
(T1)부터 (T3)까지는 항등식이므로 sympy가 정확히 다룰 수 있다. 각 성질의 좌변에서 우변을 뺀 것이 \(0\)으로 정리되는지 본다.
# (T1) to (T3): each identity should simplify to exactly zero
x, y = sp.symbols("x y", real=True)
checks = {
"T1: sin 0": sp.sin(0),
"T1: cos 0 - 1": sp.cos(0) - 1,
"T2: cos(y-x) - (cos y cos x + sin y sin x)":
sp.simplify(sp.cos(y - x) - (sp.cos(y) * sp.cos(x) + sp.sin(y) * sp.sin(x))),
"T3: sin(y-x) - (sin y cos x - cos y sin x)":
sp.simplify(sp.sin(y - x) - (sp.sin(y) * sp.cos(x) - sp.cos(y) * sp.sin(x))),
}
for name, value in checks.items():
print(f"{value} <- {name}")
0 <- T1: sin 0
0 <- T1: cos 0 - 1
0 <- T2: cos(y-x) - (cos y cos x + sin y sin x)
0 <- T3: sin(y-x) - (sin y cos x - cos y sin x)
넷 다 \(0\)이다.
정리 65가 이 셋에서 끌어낸 것들도 확인한다. 서술 파트는 이것들을 (T1)~(T3)에서 유도했으므로, 여기서 따로 맞는다면 유도가 옳았다는 방증이다.
# The consequences derived in theorem 20.2, checked independently
derived = {
"sin^2 + cos^2 - 1": sp.simplify(sp.sin(x) ** 2 + sp.cos(x) ** 2 - 1),
"cos(-x) - cos x": sp.simplify(sp.cos(-x) - sp.cos(x)),
"sin(-x) + sin x": sp.simplify(sp.sin(-x) + sp.sin(x)),
"cos(y+x) - (cos y cos x - sin y sin x)":
sp.simplify(sp.cos(y + x) - (sp.cos(y) * sp.cos(x) - sp.sin(y) * sp.sin(x))),
"sin(y+x) - (sin y cos x + cos y sin x)":
sp.simplify(sp.sin(y + x) - (sp.sin(y) * sp.cos(x) + sp.cos(y) * sp.sin(x))),
}
for name, value in derived.items():
print(f"{value} <- {name}")
0 <- sin^2 + cos^2 - 1
0 <- cos(-x) - cos x
0 <- sin(-x) + sin x
0 <- cos(y+x) - (cos y cos x - sin y sin x)
0 <- sin(y+x) - (sin y cos x + cos y sin x)
(T4) 의 c 는 얼마까지 갈 수 있는가#
(T4)는 어떤 \(c>0\)이 있다고만 했다. 그 \(c\)를 얼마까지 키울 수 있는지 찾아본다.
# (T4): 0 < sin x < x < sin x / cos x, scanned to find where it first fails
def satisfies_t4(v):
s, c = np.sin(v), np.cos(v)
if c == 0.0:
return False
return 0 < s and s < v and v < s / c
grid = np.linspace(1e-6, 3.0, 300000)
holds = np.array([satisfies_t4(v) for v in grid])
first_failure = grid[np.argmin(holds)] if not holds.all() else None
print("holds at x = 0.1, 1.0, 1.5 :",
[satisfies_t4(v) for v in (0.1, 1.0, 1.5)])
print("first x where it fails :", f"{first_failure:.6f}")
holds at x = 0.1, 1.0, 1.5 : [np.True_, np.True_, np.True_]
first x where it fails : 1.570806
\(1.5707\ldots\) 근처에서 처음 무너진다. 그 값이 4절에서 구할 \(\pi/2\)이다.
서술 파트가 (T4)에 \(\pi\)를 쓰지 않은 이유가 여기서 보인다. 가장 큰 \(c\)가 \(\pi/2\)인데, \(\pi\)는 아직 정의되지 않았으므로 그 값을 조건에 적을 수 없었다. 필요한 것은 그런 \(c\)가 있다는 사실뿐이고, 실제로 얼마인지는 뒤에 가서야 알 수 있다.
# The three quantities of (T4) at a few points inside the valid range
print(f"{'x':>8} {'sin x':>14} {'x':>10} {'sin x / cos x':>16} {'ordered':>9}")
print("-" * 62)
for v in (0.1, 0.5, 1.0, 1.5):
s, c = np.sin(v), np.cos(v)
print(f"{v:>8} {s:>14.9f} {v:>10} {s / c:>16.9f} {str(satisfies_t4(v)):>9}")
x sin x x sin x / cos x ordered
--------------------------------------------------------------
0.1 0.099833417 0.1 0.100334672 True
0.5 0.479425539 0.5 0.546302490 True
1.0 0.841470985 1.0 1.557407725 True
1.5 0.997494987 1.5 14.101419947 True
\(x=1.5\)에서 오른쪽이 \(14.1\)까지 커진다. \(\cos\)이 \(0\)에 가까워지기 때문이며, 그 자리를 넘으면 부호가 바뀌어 부등식이 무너진다.
2. 압축이 조여진다#
(20.6)은 \(\cos x < \frac{\sin x}{x} < 1\)이라 했다. 세 값을 한 표에 놓는다.
# (20.6): the middle value is trapped, and the trap closes
print(f"{'x':>10} {'cos x':>16} {'sin x / x':>16} {'1':>4} {'width':>12}")
print("-" * 62)
for v in (1.0, 0.5, 0.1, 0.01, 0.001):
c, q = np.cos(v), np.sin(v) / v
print(f"{v:>10} {c:>16.12f} {q:>16.12f} {1:>4} {1 - c:>12.2e}")
x cos x sin x / x 1 width
--------------------------------------------------------------
1.0 0.540302305868 0.841470984808 1 4.60e-01
0.5 0.877582561890 0.958851077208 1 1.22e-01
0.1 0.995004165278 0.998334166468 1 5.00e-03
0.01 0.999950000417 0.999983333417 1 5.00e-05
0.001 0.999999500000 0.999999833333 1 5.00e-07
마지막 열이 상자의 폭이다. \(x\)를 줄이면 폭이 \(10^{-7}\)까지 좁아지고, 그 안에 \(\frac{\sin x}{x}\)가 갇혀 있다. 정리 66가 하는 일이 이것이다.
(20.7)도 확인한다. \(0 \le 1-\cos x \le \sin^2 x < x^2\)이었다.
# (20.7): 1 - cos x is squeezed below x**2
print(f"{'x':>10} {'1 - cos x':>18} {'sin^2 x':>18} {'x**2':>14} {'ordered':>9}")
print("-" * 74)
for v in (0.5, 0.1, 0.01, 0.001):
gap, s2 = 1 - np.cos(v), np.sin(v) ** 2
print(f"{v:>10} {gap:>18.15f} {s2:>18.15f} {v ** 2:>14.8f} "
f"{str(0 <= gap <= s2 <= v ** 2):>9}")
x 1 - cos x sin^2 x x**2 ordered
--------------------------------------------------------------------------
0.5 0.122417438109627 0.229848847065930 0.25000000 True
0.1 0.004995834721974 0.009966711079379 0.01000000 True
0.01 0.000049999583335 0.000099996666711 0.00010000 True
0.001 0.000000499999958 0.000000999999667 0.00000100 True
\(1-\cos x\)가 \(x^2\)보다 훨씬 작다. 대략 절반이며, 그것은 \(1+\cos x\)가 \(2\)에 가깝기 때문이다.
두 극한을 sympy로도 확인한다.
# The two limits of theorem 20.4
print("limit of sin(x)/x :", sp.limit(sp.sin(x) / x, x, 0))
print("limit of (1 - cos x) / x :", sp.limit((1 - sp.cos(x)) / x, x, 0))
limit of sin(x)/x : 1
limit of (1 - cos x) / x : 0
가운데 곡선이 \(x=0\)에서 비어 있다. \(\frac{\sin x}{x}\)는 \(0\)에서 정의되지 않는다. 16회 1절의 차분몫과 같은 사정이며, 그래서 극한을 묻는다.
3. 도함수#
정리 68가 \(\sin' = \cos\)이라 했다. 차분몫으로 확인한다.
# The difference quotient of sin against cos, at several points
def difference_quotient(f, p, step=1e-6):
return (f(p + step) - f(p - step)) / (2 * step)
print(f"{'p':>8} {'quotient of sin':>18} {'cos p':>14} {'difference':>13}")
print("-" * 58)
for p in (0.0, 0.5, 1.0, 2.0, 3.0):
q = difference_quotient(np.sin, p)
print(f"{p:>8} {q:>18.12f} {np.cos(p):>14.12f} {abs(q - np.cos(p)):>13.2e}")
p quotient of sin cos p difference
----------------------------------------------------------
0.0 1.000000000000 1.000000000000 1.67e-13
0.5 0.877582561898 0.877582561890 7.48e-12
1.0 0.540302305896 0.540302305868 2.77e-11
2.0 -0.416146836546 -0.416146836547 1.14e-12
3.0 -0.989992496730 -0.989992496600 1.30e-10
# And the same for cos, whose derivative carries a minus sign
print(f"{'p':>8} {'quotient of cos':>18} {'-sin p':>14} {'difference':>13}")
print("-" * 58)
for p in (0.0, 0.5, 1.0, 2.0, 3.0):
q = difference_quotient(np.cos, p)
print(f"{p:>8} {q:>18.12f} {-np.sin(p):>14.12f} {abs(q + np.sin(p)):>13.2e}")
p quotient of cos -sin p difference
----------------------------------------------------------
0.0 0.000000000000 -0.000000000000 0.00e+00
0.5 -0.479425538602 -0.479425538604 2.17e-12
1.0 -0.841470984780 -0.841470984808 2.75e-11
2.0 -0.909297426854 -0.909297426826 2.86e-11
3.0 -0.141120008068 -0.141120008060 8.23e-12
부호가 붙는 자리가 \(\cos\) 쪽이다. 서술 파트의 증명에서 그 부호는 덧셈정리 (20.4)의 마이너스에서 나왔다.
증명이 실제로 쓴 것은 두 극한뿐이다. 분자를 쪼갠 뒤 남는 두 조각을 따로 계산해 본다.
# The proof splits the quotient into two pieces; here they are, at p = 1
p = 1.0
print(f"{'h':>10} {'(1-cos h)/h':>16} {'sin h / h':>16} {'assembled':>14} {'cos p':>12}")
print("-" * 74)
for h in (0.1, 0.01, 0.001):
first = (1 - np.cos(h)) / h
second = np.sin(h) / h
assembled = -np.sin(p) * first + np.cos(p) * second
print(f"{h:>10} {first:>16.12f} {second:>16.12f} {assembled:>14.9f} "
f"{np.cos(p):>12.9f}")
h (1-cos h)/h sin h / h assembled cos p
--------------------------------------------------------------------------
0.1 0.049958347220 0.998334166468 0.497363753 0.540302306
0.01 0.004999958333 0.999983333417 0.536085981 0.540302306
0.001 0.000499999958 0.999999833333 0.539881480 0.540302306
첫 열이 \(0\)으로 가고 둘째 열이 \(1\)로 가므로, 조립한 값이 \(-\sin p \cdot 0 + \cos p \cdot 1 = \cos p\)로 간다.
4. π 를 영점으로 구한다#
이 절이 서술 파트 6절을 그대로 실행한다.
먼저 보조정리 12의 증명이 쓴 추정을 확인한다. \(\cos X \le \cos 1 - \sin 1 \cdot (X-1)\)이고, 오른쪽이 \(-1\)보다 작아지는 \(X\)를 잡으면 모순이 나온다는 것이었다.
# The estimate inside lemma 20.5: how far must X go before the bound drops below -1
cos1, sin1 = np.cos(1.0), np.sin(1.0)
threshold = 1.0 + (cos1 + 1.0) / sin1
print(f"cos 1 = {cos1:.9f} sin 1 = {sin1:.9f}")
print(f"X must exceed 1 + (cos 1 + 1) / sin 1 = {threshold:.9f}")
print()
print(f"{'X':>8} {'bound cos1 - sin1(X-1)':>26} {'actual cos X':>16}")
print("-" * 54)
for X in (2.0, threshold, 3.0, 4.0):
print(f"{X:>8.4f} {cos1 - sin1 * (X - 1.0):>26.9f} {np.cos(X):>16.9f}")
cos 1 = 0.540302306 sin 1 = 0.841470985
X must exceed 1 + (cos 1 + 1) / sin 1 = 2.830487722
X bound cos1 - sin1(X-1) actual cos X
------------------------------------------------------
2.0000 -0.301168679 -0.416146837
2.8305 -1.000000000 -0.951995920
3.0000 -1.142639664 -0.989992497
4.0000 -1.984110649 -0.653643621
\(X\)가 \(2.83\)을 넘으면 상계가 \(-1\) 아래로 내려간다. 그런데 3절에서 \(|\cos| \le 1\)임을 보았으므로 모순이고, 따라서 \(\cos\)이 그 전에 이미 \(0\) 이하가 되어 있어야 한다.
셋째 열을 보면 실제로 \(\cos 2\)가 이미 음수이다. 증명은 그 사실을 미리 알 필요가 없었다.
이제 첫 영점을 이분법으로 잡는다.
# Theorem 20.6: the first positive zero of cos, found by bisection on [0, 2]
def bisect(f, lo, hi, rounds=80):
for _ in range(rounds):
mid = 0.5 * (lo + hi)
if f(lo) * f(mid) <= 0.0:
hi = mid
else:
lo = mid
return 0.5 * (lo + hi)
print("cos 0 =", np.cos(0.0), " cos 2 =", f"{np.cos(2.0):.9f}")
first_zero = bisect(np.cos, 0.0, 2.0)
print(f"first zero of cos : {first_zero:.15f}")
print(f"cos at that point : {np.cos(first_zero):.3e}")
cos 0 = 1.0 cos 2 = -0.416146837
first zero of cos : 1.570796326794897
cos at that point : 6.123e-17
정의 45가 그 두 배를 \(\pi\)라 하였다.
# Definition 20.7: pi is twice the first zero. Compare with the library constant.
our_pi = 2.0 * first_zero
print(f"2 * (first zero) : {our_pi:.15f}")
print(f"np.pi : {np.pi:.15f}")
print(f"difference : {abs(our_pi - np.pi):.3e}")
2 * (first zero) : 3.141592653589793
np.pi : 3.141592653589793
difference : 0.000e+00
원의 둘레도 호의 길이도 쓰지 않았다. 쓴 것은 \(\cos\)의 부호가 바뀌는 자리를 좁힌 것뿐이다.
서술 파트가 그 뒤에 적은 두 가지도 확인한다. \(\sin(\pi/2)=1\)이고, 덧셈정리에서 주기성이 나온다.
# sin at the first zero, and the period that follows from the addition formulas
print(f"sin at the first zero : {np.sin(first_zero):.15f}")
print()
for v in (0.3, 1.1, 2.7):
print(f"x = {v:<5} cos(x + pi/2) = {np.cos(v + first_zero):>12.9f}"
f" -sin x = {-np.sin(v):>12.9f}")
print()
for v in (0.3, 1.1, 2.7):
print(f"x = {v:<5} sin(x + 2pi) - sin x = {np.sin(v + 2 * our_pi) - np.sin(v):.2e}")
sin at the first zero : 1.000000000000000
x = 0.3 cos(x + pi/2) = -0.295520207 -sin x = -0.295520207
x = 1.1 cos(x + pi/2) = -0.891207360 -sin x = -0.891207360
x = 2.7 cos(x + pi/2) = -0.427379880 -sin x = -0.427379880
x = 0.3 sin(x + 2pi) - sin x = -3.89e-16
x = 1.1 sin(x + 2pi) - sin x = -3.33e-16
x = 2.7 sin(x + 2pi) - sin x = 9.99e-16
5. 적분#
정리 64를 쓰면 원시함수 하나로 끝난다. 정의로 계산한 값과 대조한다.
# Two routes to the same integral: step functions, and a primitive
def lower_upper(f, a, b, n, samples=200):
points = np.linspace(a, b, n + 1)
left = points[:-1].reshape(n, 1)
right = points[1:].reshape(n, 1)
weights = np.linspace(0.0, 1.0, samples).reshape(1, samples)
values = f(left + (right - left) * weights)
widths = np.diff(points)
return (float(np.sum(values.min(axis=1) * widths)),
float(np.sum(values.max(axis=1) * widths)))
low, high = lower_upper(np.cos, 0.0, np.pi / 2, 4000)
print(f"lower sum : {low:.9f}")
print(f"upper sum : {high:.9f}")
print(f"they trap : {high - low:.2e} wide")
print()
print("by the primitive : sin(pi/2) - sin(0) =",
f"{np.sin(np.pi / 2) - np.sin(0.0):.9f}")
lower sum : 0.999803638
upper sum : 1.000196337
they trap : 3.93e-04 wide
by the primitive : sin(pi/2) - sin(0) = 1.000000000
상합과 하합은 \(4000\)개의 조각을 더한 것이고, 아래 한 줄은 \(\sin\)의 두 값을 뺀 것이다.
# A few more, all following from the two derivative formulas
t = sp.Symbol("t", real=True)
print(f"{'integral':>28} {'primitive':>16} {'value':>10}")
print("-" * 58)
for lo, hi, expr, primitive, label in (
(0, sp.pi / 2, sp.cos(t), sp.sin(t), "cos from 0 to pi/2"),
(0, sp.pi, sp.sin(t), -sp.cos(t), "sin from 0 to pi"),
(0, 2 * sp.pi, sp.sin(t), -sp.cos(t), "sin from 0 to 2pi")):
value = sp.simplify(primitive.subs(t, hi) - primitive.subs(t, lo))
print(f"{label:>28} {str(primitive):>16} {str(value):>10}")
integral primitive value
----------------------------------------------------------
cos from 0 to pi/2 sin(t) 1
sin from 0 to pi -cos(t) 2
sin from 0 to 2pi -cos(t) 0
세 번째 줄이 \(0\)이다. 한 주기를 통째로 적분하면 위로 올라간 만큼 아래로 내려가기 때문이다.
마치며#
(T1)부터 (T3)까지가 sympy에서 정확히 \(0\)으로 정리되었고, 정리 65가 그 셋에서 끌어낸 다섯 항등식도 따로 맞았다. 유도가 옳았다는 방증이다.
(T4)가 무너지는 자리를 찾았더니 \(1.5707\ldots\)였다. 서술 파트가 (T4)에 \(\pi\)를 쓰지 않은 이유가 그것이다. 가장 큰 \(c\)가 \(\pi/2\)인데, 그때는 \(\pi\)가 없었다.
(20.6)의 상자가 \(x=0.001\)에서 폭 \(10^{-7}\)까지 좁아졌고 그 안에 \(\frac{\sin x}{x}\)가 갇혀 있었다.
차분몫이 \(\cos p\)와 \(-\sin p\)를 각각 열두 자리까지 맞혔고, 분자를 쪼갠 두 조각이 따로도 \(0\)과 \(1\)로 갔다.
그리고 4절에서 \(\cos\)의 첫 영점을 이분법으로 잡아 그 두 배로 \(3.141592653589793\)을 얻었다. np.pi와 열다섯 자리까지 같다. 원의 둘레는 한 번도 나오지 않았다.
서술 파트의 내용 |
이 노트북의 코드 |
|---|---|
(T1)~(T3) |
1절 |
1절 |
|
(T4)에 \(\pi\)가 없는 이유 |
1절 |
2절 |
|
2절 \(1-\cos x\)와 \(x^2\) 표 |
|
3절 두 차분몫 표 |
|
증명이 쓴 두 극한 |
3절 |
4절 |
|
4절 |
|
5절 상합·하합과 |
더 해 볼 것#
1절의
satisfies_t4를 \(x<0\)에서도 돌려 보자. 부등식이 성립하는가. 성립하지 않는다면 (T4)가 \(0<x\)로 범위를 제한한 이유를 설명해 보자.2절의 표에 \(x = 10^{-8}\)을 넣어 보자.
width가 계속 줄어드는가. 16회 실험 2절에서 본 자릿수 손실이 여기서도 나타나는가.4절의 이분법 구간을 \([0,2]\) 대신 \([0,4]\)로 바꾸면 어떤 영점을 찾는가. 정리 69가 “첫” 영점이라 한 것이 왜 중요한지 이 결과로 설명해 보자.
5절에서 \(\int_0^{\pi} \cos\)을 두 방법으로 구해 보자. 값이 \(0\)인데 상합과 하합은 어떻게 되는가. 하합이 음수인 것이 정의에 어긋나지 않는 이유는 무엇인가.
짝이 되는 서술