26회 실험 · 맞춘 만큼만 맞는다#
이 실험
짝 서술 |
|
확인하는 것 |
계수 공식의 \(k!\)이 어디서 나오는지 미분으로 본다. 세 함수의 다항식을 정의대로 만든다. 근사가 좋아지는 것과 나빠지는 것을 수치로 잰다 |
이 노트북은 sp.series를 쓰지 않고 시작한다.
정의 56가 계수를 \(\frac{f^{(k)}(p)}{k!}\)이라 했으므로, \(k\)번 미분하고 \(p\)를 넣고 \(k!\)로 나누면 그것이 곧 테일러 다항식이다. 정의를 그대로 옮기면 된다. 라이브러리의 급수 전개는 답을 맞춰 보는 데만 쓴다.
1절이 정리 82의 증명을 코드로 확인한다. \((x-p)^k\)를 \(j\)번 미분하고 \(x=p\)를 넣으면 \(j=k\)일 때만 살아남고 그 값이 \(k!\)이다. 표 하나가 증명의 세 경우를 그대로 보인다.
4절이 이 회차의 정직한 자리이다. 서술 파트가 “좋아진다, 나빠진다는 관찰일 뿐”이라고 했다. 코드도 관찰만 할 수 있다. 오차를 재고 그림을 그리되, 그것이 증명이 아님을 확인한다.
서술 파트의 내용 |
여기서 확인하는 방법 |
|---|---|
정리 82의 세 경우 |
1절이 \(j<k\), \(j=k\), \(j>k\)를 표로 본다 |
1절이 귀납적 정의를 그대로 옮긴다 |
|
2절이 정의대로 계수를 만든다 |
|
2절이 네 예를 대조한다 |
|
3절이 등비합의 닫힌 형태를 확인한다 |
|
3절이 유일하게 손으로 아는 \(R_n\)을 본다 |
|
5절의 관찰 두 줄 |
4절이 수치로 재고 그림으로 본다 |
0. 준비#
1. k! 이 어디서 나오는가#
정의 54을 그대로 옮긴다. 귀납적 정의이므로 코드도 재귀이다.
# def-26-factorial, verbatim: 0! = 1 and n! = n * (n-1)!
def factorial(n):
return 1 if n == 0 else n * factorial(n - 1)
print("k :", list(range(8)))
print("k! :", [factorial(k) for k in range(8)])
print("sympy :", [int(sp.factorial(k)) for k in range(8)])
k : [0, 1, 2, 3, 4, 5, 6, 7]
k! : [1, 1, 2, 6, 24, 120, 720, 5040]
sympy : [1, 1, 2, 6, 24, 120, 720, 5040]
이제 정리 82의 증명을 확인한다. \((x-p)^k\)를 \(j\)번 미분하고 \(x=p\)를 넣는다.
# The three cases of the proof: j < k, j = k, j > k
k = 4
print(f"differentiating (x - p)**{k} and then setting x = p")
print()
print(f"{'j':>4} {'value at x = p':>18} {'case':>10}")
print("-" * 36)
for j in range(k + 3):
derivative = sp.diff((x - p) ** k, x, j)
value = sp.simplify(derivative.subs(x, p))
case = "j < k" if j < k else ("j = k" if j == k else "j > k")
print(f"{j:>4} {str(value):>18} {case:>10}")
differentiating (x - p)**4 and then setting x = p
j value at x = p case
------------------------------------
0 0 j < k
1 0 j < k
2 0 j < k
3 0 j < k
4 24 j = k
5 0 j > k
6 0 j > k
\(j=4\)인 줄만 \(0\)이 아니고 그 값이 \(24 = 4!\)이다.
세 경우가 그대로 보인다. \(j<k\)이면 \((x-p)^{k-j}\)가 남아 \(x=p\)에서 사라지고, \(j>k\)이면 상수를 미분해 사라지며, \(j=k\)이면 지수가 차례로 내려와 곱해진 \(k!\)이 남는다.
# Where the k! comes from: the exponents come down one at a time
k = 5
expression = (x - p) ** k
for j in range(1, k + 1):
expression = sp.diff(expression, x)
print(f"after {j} differentiations : {expression}")
print()
print(f"the accumulated factor is {factorial(k)} = {k}!")
after 1 differentiations : 5*(-p + x)**4
after 2 differentiations : 20*(-p + x)**3
after 3 differentiations : 60*(-p + x)**2
after 4 differentiations : -120*p + 120*x
after 5 differentiations : 120
the accumulated factor is 120 = 5!
2. 정의대로 다항식을 만든다#
정의 56를 그대로 옮긴다. sp.series를 쓰지 않는다.
# def-26-taylor, verbatim
def taylor_polynomial(expr, n, point=0):
return sum(sp.diff(expr, x, k).subs(x, point) / sp.factorial(k)
* (x - point) ** k
for k in range(n + 1))
네 예를 계산한다.
# The four examples of section 3, built from the definition
examples = [("exp", sp.exp(x), 5),
("sin", sp.sin(x), 7),
("cos", sp.cos(x), 6),
("1/(1-x)", 1 / (1 - x), 5)]
for name, expr, n in examples:
poly = sp.expand(taylor_polynomial(expr, n))
print(f"{name:>9} T_{n} = {poly}")
exp T_5 = x**5/120 + x**4/24 + x**3/6 + x**2/2 + x + 1
sin T_7 = -x**7/5040 + x**5/120 - x**3/6 + x
cos T_6 = -x**6/720 + x**4/24 - x**2/2 + 1
1/(1-x) T_5 = x**5 + x**4 + x**3 + x**2 + x + 1
계수만 따로 본다. 서술 파트 3절이 각 함수마다 규칙을 찾은 것이 이 표이다.
# The coefficients, side by side
print(f"{'k':>3} {'exp':>10} {'sin':>10} {'cos':>10} {'1/(1-x)':>10}")
print("-" * 48)
for k in range(7):
row = [sp.diff(expr, x, k).subs(x, 0) / sp.factorial(k)
for _, expr, _ in examples]
print(f"{k:>3} " + " ".join(f"{str(v):>10}" for v in row))
k exp sin cos 1/(1-x)
------------------------------------------------
0 1 0 1 1
1 1 1 0 1
2 1/2 0 -1/2 1
3 1/6 -1/6 0 1
4 1/24 0 1/24 1
5 1/120 1/120 0 1
6 1/720 0 -1/720 1
네 열이 각각 다른 규칙을 보인다.
\(\exp\)는 모든 계수가 \(\frac{1}{k!}\)이다. 도함수가 자기 자신이기 때문이다.
\(\sin\)과 \(\cos\)은 한 칸씩 건너뛰며 부호가 바뀐다. 네 번 미분하면 제자리로 돌아오기 때문이고, 20회 정리 65의 홀·짝 대칭이 그대로 나타난다.
\(\frac{1}{1-x}\)는 계수가 모두 \(1\)이다. \(f^{(k)}(0)=k!\)이라 \(k!\)이 약분된다.
라이브러리와 맞춰 본다.
# Only now do we use the library, and only as a check
for name, expr, n in examples:
ours = sp.expand(taylor_polynomial(expr, n))
theirs = sp.expand(sp.series(expr, x, 0, n + 1).removeO())
print(f"{name:>9} : agree = {sp.simplify(ours - theirs) == 0}")
exp : agree = True
sin : agree = True
cos : agree = True
1/(1-x) : agree = True
3. 나머지항을 손으로 아는 경우#
정의 57가 \(R_n(x)=f(x)-T_n(x)\)라 했고, 서술 파트는 그 크기를 말하지 않았다.
그런데 \(\frac{1}{1-x}\)에서만 손으로 적을 수 있었다. (26.7)의 닫힌 형태가 있기 때문이다.
# (26.7): the geometric sum has a closed form
for n in (3, 5, 8):
left = sum(x ** k for k in range(n + 1))
right = (1 - x ** (n + 1)) / (1 - x)
print(f"n = {n} : closed form matches = {sp.simplify(left - right) == 0}")
n = 3 : closed form matches = True
n = 5 : closed form matches = True
n = 8 : closed form matches = True
# And therefore the remainder is known exactly for this one function
n = 4
remainder = sp.simplify(1 / (1 - x) - taylor_polynomial(1 / (1 - x), n))
print("R_4(x) =", sp.simplify(remainder))
print("matches x**(n+1)/(1-x) :",
sp.simplify(remainder - x ** (n + 1) / (1 - x)) == 0)
R_4(x) = -x**5/(x - 1)
matches x**(n+1)/(1-x) : True
다른 세 함수에는 이런 형태가 없다. 확인해 둔다.
# For exp there is no closed form for the remainder in elementary terms
n = 4
for name, expr, _ in examples[:3]:
rest = sp.simplify(expr - taylor_polynomial(expr, n))
print(f"{name:>5} R_{n}(x) = {rest}")
exp R_4(x) = -x**4/24 - x**3/6 - x**2/2 - x + exp(x) - 1
sin R_4(x) = x**3/6 - x + sin(x)
cos R_4(x) = -x**4/24 + x**2/2 + cos(x) - 1
세 줄 모두 원래 함수에서 다항식을 뺀 것을 다시 적었을 뿐이다. 새 정보가 없다.
\(\frac{1}{1-x}\)가 특별했던 것이지 일반적인 것이 아니다.
4. 관찰한다#
서술 파트 5절이 두 줄을 관찰했다. 각각 수치로 잰다.
# Observation 1: at a fixed x, raising n improves the approximation
target = np.exp(0.1)
print("approximating e**0.1")
print(f"{'n':>3} {'T_n(0.1)':>18} {'error':>14} {'ratio':>10}")
print("-" * 50)
partial = 0.0
previous = None
for n in range(6):
partial += 0.1 ** n / factorial(n)
error = abs(target - partial)
ratio = "" if previous is None else f"{previous / error:.1f}"
print(f"{n:>3} {partial:>18.12f} {error:>14.3e} {ratio:>10}")
previous = error
print(f"\n e**0.1 = {target:.12f}")
approximating e**0.1
n T_n(0.1) error ratio
--------------------------------------------------
0 1.000000000000 1.052e-01
1 1.100000000000 5.171e-03 20.3
2 1.105000000000 1.709e-04 30.3
3 1.105166666667 4.251e-06 40.2
4 1.105170833333 8.474e-08 50.2
5 1.105170916667 1.409e-09 60.1
e**0.1 = 1.105170918076
마지막 열이 오차가 줄어드는 비율이다. \(20\), \(30\), \(40\)으로 커진다. 서술 파트가 적은 대로 다음 항의 분모에 \((n+1)!\)이 있기 때문이다.
# Observation 2: at a fixed n, moving away from p makes it worse
n = 3
print(f"the degree {n} polynomial of exp, at several points")
print(f"{'x':>6} {'T_3(x)':>16} {'exp(x)':>16} {'error':>14}")
print("-" * 56)
for v in (0.1, 0.5, 1.0, 2.0, 4.0):
approx = sum(v ** k / factorial(k) for k in range(n + 1))
print(f"{v:>6} {approx:>16.9f} {np.exp(v):>16.9f} "
f"{abs(np.exp(v) - approx):>14.6f}")
the degree 3 polynomial of exp, at several points
x T_3(x) exp(x) error
--------------------------------------------------------
0.1 1.105166667 1.105170918 0.000004
0.5 1.645833333 1.648721271 0.002888
1.0 2.666666667 2.718281828 0.051615
2.0 6.333333333 7.389056099 1.055723
4.0 23.666666667 54.598150033 30.931483
\(x=4\)에서 오차가 \(20\)을 넘는다. 같은 다항식인데 \(x\)에 따라 쓸모가 완전히 달라진다.
코드도 관찰만 한다
왼쪽 그림에서 차수가 높아질수록 겹치는 구간이 넓어지고, 오른쪽 그림에서 오차가 곧게 떨어진다.
둘 다 증명이 아니다. 왼쪽은 \([-8,8]\)의 \(600\)개 점에서 그린 것이고 오른쪽은 \(n\)을 \(8\)까지만 본 것이다.
그리고 이 두 그림은 \(\sin\)과 \(\exp\)에 대한 것이다. 다른 함수에서도 그런지는 이 노트북이 말할 수 없다. 서술 파트가 흔한 오해 상자에서 “\(n\)을 키우면 언제나 좋아지는가 — 이 회차의 도구로는 알 수 없다”고 적은 것이 그것이다.
\(R_n\)의 크기를 통제하는 것이 다음 회차의 일이다.
마치며#
정의 54을 재귀 함수로 옮겼고 sympy 와 값이 맞았다.
\((x-p)^4\)를 \(j\)번 미분하고 \(x=p\)를 넣었더니 \(j=4\)인 줄만 살아남았고 그 값이 \(24=4!\)이었다. 정리 82 증명의 세 경우가 표 하나에 그대로 나온다. 미분을 한 번씩 따라가 보면 지수가 \(5,4,3,2,1\)로 내려오며 곱해지는 것도 보인다.
sp.series를 쓰지 않고 정의 56만으로 네 다항식을 만들었고, 라이브러리와 맞았다. 계수 표에서 네 함수가 각각 다른 규칙을 보였다.
\(\frac{1}{1-x}\)에서만 \(R_n(x)=\frac{x^{n+1}}{1-x}\)을 얻었다. 나머지 셋은 정의를 다시 적은 것에 지나지 않았다.
오차를 두 방향으로 쟀다. \(n\)을 키우면 \(20\), \(30\), \(40\)배씩 줄었고, \(x\)를 \(4\)까지 밀면 오차가 \(20\)을 넘었다. 둘 다 관찰이며 증명이 아니다.
더 해 볼 것#
1절의 \(k\)를 \(6\)으로 바꾸고 \(j\)를 \(0\)부터 \(8\)까지 돌려 보자. \(0\)이 아닌 줄이 몇 개인가.
2절의
taylor_polynomial에 \(p=1\)을 넣어 \(\log x\)의 테일러 다항식을 구해 보자. 계수에 어떤 규칙이 있는가. \(p=0\)에서는 왜 안 되는가.4절의 둘째 표를 \(n=6\)으로 다시 뽑아 보자. \(x=4\)에서 오차가 얼마나 줄어드는가. \(x=8\)에서는 어떠한가.
\(f(x)=\frac{1}{1+x^2}\)의 매클로린 다항식을 \(n=6\)까지 구해 보자. 계수에 규칙이 보이는가. \(x=2\)에서 \(T_6\)과 \(f\)를 비교하면 무슨 일이 일어나는가.
짝이 되는 서술