17회 실험 · 되돌리는 함수#
이 실험
짝 서술 |
|
확인하는 것 |
역함수를 이분법으로 직접 만들어 접선 기울기가 서로 역수임을 본다. (17.3)이 실제로 오류를 내는 장면을 만든다. 연쇄법칙과 유리수 지수 공식을 정의와 대조한다 |
이 노트북의 중심은 2절이다.
서술 파트가 (17.3)을 두고 “항등식이 아니다”라고 했다. 말로만 하면 지나치기 쉬운 지적인데, 코드는 그것을 오류로 만든다. \(g\)가 \(p\) 근처에서 상수인 함수를 만들어 (17.3)의 오른쪽을 계산시키면 ZeroDivisionError가 난다. 같은 자리에서 보조정리 9의 \(F\)는 멀쩡한 값을 준다.
1절은 역함수를 손으로 만든다. 서술 파트 1절이 \(f^{-1}\)의 존재를 정리로 얻었지만 그 정리는 값을 알려 주지 않는다. 값을 얻으려면 따로 계산해야 하고, 그 계산이 이분법이다. 14회 실험에서 쓴 것과 같은 절차인데 이번에는 근이 아니라 역함숫값을 찾는다.
서술 파트의 내용 |
여기서 확인하는 방법 |
|---|---|
1절이 이분법으로 \(f^{-1}\)의 값을 구한다 |
|
(17.8)의 역수 관계 |
1절이 두 기울기를 곱해 \(1\)이 나오는지 본다 |
(17.3)이 항등식이 아니다 |
2절이 오류를 실제로 받는다 |
보조정리 9의 \(F\) |
2절이 같은 자리에서 값을 얻는다 |
3절이 전개한 결과와 공식을 대조한다 |
|
4절이 이분법 역함수의 수치 미분과 맞춘다 |
0. 준비#
이 회차에서 쓸 것만 불러온다. 상합·하합을 그리지 않으므로 그쪽 색 상수는 두지 않는다.
1. 역함수를 직접 만든다#
정리 52와 정리 53는 \(f^{-1}\)이 있다고 말한다. 그 값이 얼마인지는 말하지 않는다.
값을 구하려면 \(f(x) = y\)를 푸는 절차가 따로 필요하다. 이분법이 그 절차이며, 그 절차가 옳게 도는 근거가 정확히 정리의 가정이다. \(f\)가 순증가하므로 중점에서의 함숫값을 \(y\)와 견주면 근이 어느 쪽에 있는지 알 수 있다.
# Bisection for the inverse: it works because f is continuous and strictly increasing
def inverse_by_bisection(f, y, lo, hi, rounds=60):
for _ in range(rounds):
mid = 0.5 * (lo + hi)
if f(mid) < y:
lo = mid
else:
hi = mid
return 0.5 * (lo + hi)
cube = lambda t: t ** 3
print(f"{'y':>8} {'inverse value':>18} {'cubed back':>14}")
print("-" * 42)
for y in (-8.0, -1.0, 0.125, 1.0, 8.0):
x = inverse_by_bisection(cube, y, -3.0, 3.0)
print(f"{y:>8} {x:>18.12f} {x ** 3:>14.10f}")
y inverse value cubed back
------------------------------------------
-8.0 -2.000000000000 -8.0000000000
-1.0 -1.000000000000 -1.0000000000
0.125 0.500000000000 0.1250000000
1.0 1.000000000000 1.0000000000
8.0 2.000000000000 8.0000000000
마지막 열이 첫 열로 돌아온다. (17.1)의 \(f(f^{-1}(y)) = y\)가 이것이다.
\(f\)가 순증가하지 않으면 이 절차가 무너진다. \(f(x)=x^2\)을 \([-3,3]\)에서 넣으면 f(mid) < y 비교가 뜻을 잃는다. 왼쪽에도 오른쪽에도 답이 있기 때문이다. 정의 33이 가정에 들어간 이유가 그것이다.
두 기울기가 서로 역수이다#
(17.8)은 \((f^{-1})'(f(p)) = 1/f'(p)\)라고 한다. 두 값을 곱하면 \(1\)이어야 한다.
\(f^{-1}\)의 도함수값은 공식을 쓰지 않고 차분몫으로 구한다. 그래야 대조가 된다.
# f' from the formula, the inverse's derivative from a difference quotient
def numeric_derivative(f, point, step=1e-5):
return (f(point + step) - f(point - step)) / (2 * step)
inverse_cube = lambda v: inverse_by_bisection(cube, v, -4.0, 4.0)
print(f"{'p':>6} {'q = f(p)':>10} {'f prime at p':>14} "
f"{'inv prime at q':>16} {'product':>10}")
print("-" * 60)
for p in (0.5, 1.0, 2.0):
q = p ** 3
slope = 3 * p ** 2
inverse_slope = numeric_derivative(inverse_cube, q)
print(f"{p:>6} {q:>10.3f} {slope:>14.6f} "
f"{inverse_slope:>16.6f} {slope * inverse_slope:>10.6f}")
p q = f(p) f prime at p inv prime at q product
------------------------------------------------------------
0.5 0.125 0.750000 1.333333 1.000000
1.0 1.000 3.000000 0.333333 1.000000
2.0 8.000 12.000000 0.083333 1.000000
마지막 열이 모두 \(1\)이다.
그림으로 보면 왜 그런지가 드러난다. 두 그래프가 직선 \(y=x\)에 대하여 대칭이고, 대칭시키면 기울기가 뒤집힌다.
두 점선의 기울기가 \(3p^2\)과 \(1/(3p^2)\)이다. \(p=1.2\)에서 각각 \(4.32\)와 약 \(0.231\)이다.
2. 쪼개기가 무너지는 자리#
서술 파트 2절이 지적한 것을 코드로 확인한다. \(g(x)=g(p)\)인 \(x\)에서 (17.3)의 오른쪽이 아무것도 뜻하지 않는다.
\(p\) 근처에서 통째로 상수인 \(g\)를 만든다. \(|x|\le1\)에서 \(0\)이고 바깥에서 \((|x|-1)^3\)인 함수이다. 이 함수는 \(x=0\)에서 미분 가능하고 \(g'(0)=0\)이다.
# g is identically zero on [-1,1]: every x near p = 0 has g(x) = g(p)
def g_flat(x):
size = abs(x)
return 0.0 if size <= 1.0 else (size - 1.0) ** 3
f_outer = lambda y: y ** 2
point = 0.0
value_at_point = g_flat(point)
print("g(p) =", value_at_point)
print("g(x) for x near p :", [g_flat(x) for x in (0.5, 0.1, 0.01)])
g(p) = 0.0
g(x) for x near p : [0.0, 0.0, 0.0]
세 값이 모두 \(g(p)\)와 같다. 이제 (17.3)의 오른쪽 첫 분수를 계산시킨다.
# The first factor of (17.3) divides by g(x) - g(p), which is zero here
for x in (0.5, 0.1, 0.01):
try:
first = ((f_outer(g_flat(x)) - f_outer(value_at_point))
/ (g_flat(x) - value_at_point))
print(f"x = {x:<6} first factor = {first}")
except ZeroDivisionError as error:
print(f"x = {x:<6} ZeroDivisionError: {error}")
x = 0.5 ZeroDivisionError: float division by zero
x = 0.1 ZeroDivisionError: float division by zero
x = 0.01 ZeroDivisionError: float division by zero
세 줄 모두 오류이다. \(p\)에 가까워질수록 나아지지도 않는다. 이 \(g\)는 \(p\) 근처 어디서나 상수이므로 (17.3)은 \(p\)의 어떤 근방에서도 쓸 수 없다.
이제 보조정리 9의 \(F\)를 만들어 같은 자리를 계산한다. \(f(y)=y^2\)이므로 \(f'(y)=2y\)이고 \(q=g(p)=0\)에서 \(f'(q)=0\)이다.
# The extended difference quotient (17.4) has a value at y = q as well
derivative_at_q = 2 * value_at_point
def F(y):
if y == value_at_point:
return derivative_at_q
return (f_outer(y) - f_outer(value_at_point)) / (y - value_at_point)
print(f"{'x':>8} {'F(g(x))':>12} {'(17.7) right side':>20} "
f"{'direct quotient':>18}")
print("-" * 62)
for x in (0.5, 0.1, 0.01):
by_17_7 = F(g_flat(x)) * (g_flat(x) - value_at_point) / (x - point)
direct = (f_outer(g_flat(x)) - f_outer(value_at_point)) / (x - point)
print(f"{x:>8} {F(g_flat(x)):>12.6f} {by_17_7:>20.6f} {direct:>18.6f}")
x F(g(x)) (17.7) right side direct quotient
--------------------------------------------------------------
0.5 0.000000 0.000000 0.000000
0.1 0.000000 0.000000 0.000000
0.01 0.000000 0.000000 0.000000
오류가 없다. 그리고 (17.7)의 오른쪽과 직접 계산한 차분몫이 같은 값을 준다.
값이 모두 \(0\)인 것도 맞다. \(g\)가 \(p\) 근처에서 상수이므로 \(f(g(x))\)도 상수이고 차분몫이 \(0\)이다. 연쇄법칙이 주는 값도 \(f'(g(p))\,g'(p) = 0 \cdot 0 = 0\)이다.
오류가 나는 것이 요점이다
이 절의 셀 하나는 일부러 예외를 잡는다. ZeroDivisionError가 서술 파트의 지적을 그대로 옮긴 것이기 때문이다.
증명에서 (17.3)을 쓰면 이 오류에 해당하는 자리가 조용히 지나간다. 종이 위에서는 \(0\)으로 나누어도 아무 일이 일어나지 않는다. 코드는 멈춘다.
16회 실험 1절에서 차분몫에 \(x=p\)를 넣어 같은 오류를 받았다. 정의역 밖을 건드리면 코드가 알려 준다.
3. 연쇄법칙을 정의와 대조한다#
\(g(x)=3x+1\), \(f(y)=y^2\)에서 두 계산을 나란히 놓는다. 하나는 전개해서 미분한 것이고 다른 하나는 (17.6)이다.
# Expand first and differentiate, versus the chain rule
t = sp.Symbol("t", real=True)
outer, inner = t ** 2, 3 * t + 1
composed = sp.expand(outer.subs(t, inner))
by_expanding = sp.expand(sp.diff(composed, t))
by_chain = sp.expand(sp.diff(outer, t).subs(t, inner) * sp.diff(inner, t))
print("composed :", composed)
print("expand then differentiate :", by_expanding)
print("f'(g(x)) * g'(x) :", by_chain)
print("agree :", sp.simplify(by_expanding - by_chain) == 0)
composed : 9*t**2 + 6*t + 1
expand then differentiate : 18*t + 6
f'(g(x)) * g'(x) : 18*t + 6
agree : True
서술 파트의 흔한 오해 상자가 든 잘못된 형태도 계산해 본다.
# The wrong reading: substituting g'(x) into f' instead of g(x)
wrong = sp.diff(outer, t).subs(t, sp.diff(inner, t))
print("f'(g'(x)) :", wrong, " <- a constant, not a function of x")
print("correct :", by_chain)
f'(g'(x)) : 6 <- a constant, not a function of x
correct : 18*t + 6
\(f'(g'(x))\)가 상수 \(6\)이다. 옳은 답은 \(18x+6\)이므로 \(x\)에 따라 변하지도 않는다.
16회 마치며가 남긴 \((x^2+1)^{10}\)을 처리한다.
# The debt from lecture 16: no expansion needed
big = (t ** 2 + 1) ** 10
print("derivative, factored :", sp.factor(sp.diff(big, t)))
print("degree if expanded :", sp.degree(sp.expand(big), t),
"-> its derivative has degree", sp.degree(sp.expand(sp.diff(big, t)), t))
derivative, factored : 20*t*(t**2 + 1)**9
degree if expanded : 20 -> its derivative has degree 19
전개하면 차수 \(20\)인 다항식이고 항이 열한 개이다. 연쇄법칙은 그것을 만들지 않고 \(20x(x^2+1)^9\)을 곧바로 준다.
4. n제곱근과 유리수 지수#
정의 35가 \(y^{1/n}\)을 \(x^n\)의 역함수로 정의했다. 그 정의를 그대로 코드로 옮기면 1절의 이분법이다.
# y**(1/n) is by definition the inverse of x**n on [0, b]
def nth_root(y, n, bound=100.0):
return inverse_by_bisection(lambda s: s ** n, y, 0.0, bound)
print(f"{'y':>8} {'n':>4} {'y**(1/n)':>16} {'raised back':>14}")
print("-" * 46)
for y, n in ((2.0, 2), (9.0, 2), (8.0, 3), (81.0, 4)):
root = nth_root(y, n)
print(f"{y:>8} {n:>4} {root:>16.12f} {root ** n:>14.10f}")
y n y**(1/n) raised back
----------------------------------------------
2.0 2 1.414213562373 2.0000000000
9.0 2 3.000000000000 9.0000000000
8.0 3 2.000000000000 8.0000000000
81.0 4 3.000000000000 81.0000000000
첫 줄이 \(2^{1/2}\)이다. 14회 예제 8가 \((1,2)\) 안에 있다고만 한 그 수이며, 이제 자릿수가 보인다.
(17.10)이 맞는지 확인한다. 왼쪽은 이분법 역함수의 수치 미분이고 오른쪽은 공식이다.
# (17.10): the derivative of the inverse, numerically versus by formula
n = 3
print(f"{'y':>8} {'numeric':>16} {'formula (17.10)':>18} {'difference':>14}")
print("-" * 60)
for y in (1.0, 8.0, 27.0):
numeric = numeric_derivative(lambda v: nth_root(v, n), y, step=1e-6)
formula = 1.0 / (n * nth_root(y, n) ** (n - 1))
print(f"{y:>8} {numeric:>16.10f} {formula:>18.10f} "
f"{abs(numeric - formula):>14.2e}")
y numeric formula (17.10) difference
------------------------------------------------------------
1.0 0.3333333334 0.3333333333 4.66e-11
8.0 0.0833333331 0.0833333333 2.10e-10
27.0 0.0370370370 0.0370370370 4.42e-11
차이가 \(10^{-9}\) 아래이다. 남은 차이는 이분법의 반올림과 수치 미분의 자릿수 손실에서 오며, 16회 실험 2절이 그 손실을 다루었다.
이제 (17.11)을 sympy로 확인한다.
# (17.11): rational exponents, and the value at a point
y = sp.Symbol("y", positive=True)
print(f"{'exponent':>10} {'derivative':>22} {'at y = 64':>14}")
print("-" * 48)
for m, k in ((1, 2), (1, 3), (2, 3), (3, 2)):
expr = y ** sp.Rational(m, k)
derivative = sp.diff(expr, y)
print(f"{f'{m}/{k}':>10} {str(derivative):>22} "
f"{str(derivative.subs(y, 64)):>14}")
exponent derivative at y = 64
------------------------------------------------
1/2 1/(2*sqrt(y)) 1/16
1/3 1/(3*y**(2/3)) 1/48
2/3 2/(3*y**(1/3)) 1/6
3/2 3*sqrt(y)/2 12
첫 줄이 \(\dfrac{1}{2\sqrt y}\)이다. 서술 파트가 (17.10)에서 \(n=2\)로 얻은 것과 같다.
\(y=64\)에서 \(y^{1/2}\)의 도함수값이 \(1/16\)이다. 손으로 확인하면 \(2\sqrt{64} = 16\)이므로 맞다.
# The two routes to the same number at y = 64
print("from (17.10) with n = 2 :", 1.0 / (2 * nth_root(64.0, 2)))
print("from (17.11) with m/n = 1/2 :",
float(sp.diff(y ** sp.Rational(1, 2), y).subs(y, 64)))
from (17.10) with n = 2 : 0.0625
from (17.11) with m/n = 1/2 : 0.0625
마치며#
이분법으로 \(x^3\)의 역함수를 만들고 세제곱해서 원래 값으로 돌아오는 것을 확인했다. 정리는 존재를 주고 이분법이 값을 준다. 14회 실험에서 한 구분이 여기서 다시 나온다.
두 기울기를 곱했더니 \(1\)이 나왔다. \(f^{-1}\)의 도함수값은 공식이 아니라 차분몫으로 구했으므로 (17.8)에 대한 독립적인 확인이다.
(17.3)의 오른쪽을 계산시켰더니 ZeroDivisionError가 났다. \(p\)에 가까운 세 점 모두에서 그랬다. 같은 자리에서 보조정리 9의 \(F\)는 값을 주었고, (17.7)의 오른쪽이 직접 계산한 차분몫과 일치했다. 분모를 없앤 것이 실제로 효과가 있었다.
\(f'(g'(x))\)가 상수 \(6\)을 주는 것도 보았다. 옳은 답 \(18x+6\)과 달리 \(x\)에 따라 변하지도 않는다.
\(2^{1/2}\)의 자릿수를 얻었고, (17.10)의 공식이 수치 미분과 \(10^{-9}\) 아래로 일치했다.
서술 파트의 내용 |
이 노트북의 코드 |
|---|---|
정의 33이 가정인 이유 |
1절 |
1절 |
|
1절 |
|
(17.3)이 항등식이 아니다 |
2절 |
(17.4)의 \(F\) |
2절 |
2절 |
|
3절 |
|
\(f'(g'(p))\) 는 틀렸다 |
3절 |
4절 |
|
4절 |
더 해 볼 것#
1절의
inverse_by_bisection에 \(f(x)=x^2\)과 구간 \([-3,3]\)을 넣어 \(y=4\)를 찾아 보자. 어느 답이 나오는가. 구간을 \([0,3]\)으로 바꾸면 달라지는가. 정리 52가 순증가를 요구하는 이유를 이 결과로 설명해 보자.2절의 \(g\)를 \(|x|\le 0.001\)에서만 상수인 함수로 바꾸어 보자. (17.3)이 무너지는 \(x\)의 범위가 줄어드는가. 그래도 \(p\)의 모든 근방에서 무너지는가.
3절의 \(f\)와 \(g\)를 바꾸어 \(g(y)=y^2\)을 바깥에, \(f(x)=3x+1\)을 안에 놓아 보자. \((f\circ g)'\)과 \((g\circ f)'\)이 같은가. 연쇄법칙에서 순서가 어디에 나타나는가.
4절의
nth_root에 \(y=0\)을 넣고 \(n=3\)으로 (17.10)을 계산해 보면 무슨 일이 일어나는가. 그것이 예제 12와 어떻게 이어지는가.
짝이 되는 서술