19회 실험 · 두 절차가 만난다#
이 실험
짝 서술 |
|
확인하는 것 |
두 정의를 각각 코드로 옮겨 공통 부분이 없음을 본다. 누적 함수를 상합·하합으로 만들고 그것을 차분몫으로 미분해 원래 함수가 나오는지 확인한다 |
서술 파트 1절이 두 정의를 나란히 놓고 “서로를 언급하지 않는다”고 하였다. 코드는 그 말을 검사할 수 있다.
1절에서 적분과 도함수를 각각 정의에 따라 구현한다. 하나는 분할을 만들고 계단함수의 넓이를 더해 상합과 하합을 좁히고, 다른 하나는 차분몫에 작은 수를 넣는다. 두 함수는 변수 하나도 공유하지 않는다.
그다음이 요점이다. 3절에서 적분으로 만든 함수를 미분한다. 상합과 하합으로 계산한 값들을 차분몫에 넣는 것이며, 두 절차 사이에 아무 연결이 없는데도 결과가 원래 함수와 맞는다. 정리 63이 주장하는 것이 그것이다.
4절이 반대 방향을 본다. 7회가 열 개의 식으로 얻은 값을 정리 64가 한 줄로 준다.
서술 파트의 내용 |
여기서 확인하는 방법 |
|---|---|
두 정의는 서로를 언급하지 않는다 |
1절의 두 함수가 코드를 공유하지 않는다 |
2절이 위끝과 아래끝을 바꾸어 부호를 본다 |
|
2절이 \(\lvert A(x)-A(p)\rvert\)와 \(M\lvert x-p\rvert\)를 견준다 |
|
3절이 \(A\)를 수치로 미분해 \(f\)와 맞춘다 |
|
불연속점에서 무너진다 |
3절이 도약함수의 좌우 차분몫을 따로 구한다 |
4절이 7회의 \(L_n\), \(U_n\)과 한 줄 계산을 나란히 놓는다 |
0. 준비#
이 회차에서 쓸 것만 불러온다. 상합과 하합을 다시 만드므로 그쪽 색 상수를 둔다.
1. 두 정의를 각각 코드로#
정의 15을 그대로 옮긴다. 분할을 만들고, 각 부분구간에서 최솟값과 최댓값으로 계단함수를 만들고, 그 넓이를 더한다.
# Definition 7: partitions, step functions below and above, and their areas
def partition(a, b, n):
return np.linspace(a, b, n + 1)
def lower_upper(f, a, b, n, samples=200):
points = partition(a, b, n)
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)))
정의 30도 그대로 옮긴다. 차분몫에 \(p\)에 가까운 값을 넣는다.
# Definition 16: the difference quotient, evaluated near p
def derivative_by_definition(f, p, step=1e-6):
return (f(p + step) - f(p - step)) / (2 * step)
두 함수를 나란히 놓고 보자. lower_upper에 있는 것은 partition, min, max, sum이다. derivative_by_definition에 있는 것은 뺄셈과 나눗셈뿐이다. 공유하는 것은 인자 이름 f와 p 같은 것뿐이며, 계산이 겹치는 부분이 하나도 없다.
# Both procedures applied to the same f, producing two unrelated numbers
square = lambda t: t ** 2
low, high = lower_upper(square, 0.0, 1.0, 2000)
print(f"lower sum on [0,1] : {low:.9f}")
print(f"upper sum on [0,1] : {high:.9f}")
print(f"they trap : {high - low:.2e} wide")
print()
print(f"difference quotient at p = 3 : "
f"{derivative_by_definition(square, 3.0):.9f}")
lower sum on [0,1] : 0.333083375
upper sum on [0,1] : 0.333583375
they trap : 5.00e-04 wide
difference quotient at p = 3 : 6.000000001
왼쪽 세 줄은 넓이를 좁힌 것이고 마지막 줄은 기울기를 잰 것이다. 두 결과 사이에 아무 관계도 보이지 않는다.
2. 누적 함수#
이제 위끝을 움직인다. 정의 42을 코드로 옮기되, 정의 41의 약속도 함께 넣는다.
# Definition 19.2 with the orientation convention (19.1) built in
def integral_by_definition(f, a, b, n=4000):
if a == b:
return 0.0
if b < a:
return -integral_by_definition(f, b, a, n)
low, high = lower_upper(f, a, b, n)
return 0.5 * (low + high)
def accumulation(f, a):
return lambda x: integral_by_definition(f, a, x)
A = accumulation(square, 0.0)
\(A\)가 무엇인지 값으로 확인한다. 서술 파트 6절이 \(\int_0^b x^2 = b^3/3\)이라 했으므로 \(A(x)\)는 \(x^3/3\)이어야 한다.
# A(x) built from step functions, against the closed form
print(f"{'x':>8} {'A(x)':>16} {'x**3 / 3':>16} {'difference':>14}")
print("-" * 58)
for x in (0.0, 0.5, 1.0, 1.5, 2.0):
print(f"{x:>8} {A(x):>16.9f} {x ** 3 / 3:>16.9f} "
f"{abs(A(x) - x ** 3 / 3):>14.2e}")
x A(x) x**3 / 3 difference
----------------------------------------------------------
0.0 0.000000000 0.000000000 0.00e+00
0.5 0.041666668 0.041666667 1.30e-09
1.0 0.333333344 0.333333333 1.04e-08
1.5 1.125000035 1.125000000 3.52e-08
2.0 2.666666750 2.666666667 8.33e-08
아직 \(A\)가 \(x^3/3\)이라는 것을 증명한 것은 아니다. 표는 다섯 점에서 값이 맞는다는 것을 보일 뿐이며, 각 값도 \(n=4000\)의 상합과 하합을 평균한 근삿값이다.
방향 규약을 확인한다.
# (19.1): the empty interval gives zero, and swapping the ends flips the sign
print("A(0) :", integral_by_definition(square, 0.0, 0.0))
print("integral from 2 to 1 :", f"{integral_by_definition(square, 2.0, 1.0):.9f}")
print("minus 1 to 2 :", f"{-integral_by_definition(square, 1.0, 2.0):.9f}")
A(0) : 0.0
integral from 2 to 1 : -2.333333344
minus 1 to 2 : -2.333333344
립시츠 부등식#
(19.3)은 \(|A(x)-A(p)| \le M|x-p|\)라고 했다. \([0,2]\)에서 \(|x^2| \le 4\)이므로 \(M=4\)이다.
# (19.3): A cannot be steeper than M
bound = 4.0
print(f"{'p':>6} {'x':>6} {'|A(x)-A(p)|':>16} {'M|x-p|':>12} {'holds':>7}")
print("-" * 52)
for p, x in ((0.5, 0.6), (1.0, 1.2), (1.9, 2.0), (0.0, 2.0)):
left = abs(A(x) - A(p))
right = bound * abs(x - p)
print(f"{p:>6} {x:>6} {left:>16.9f} {right:>12.4f} {str(left <= right):>7}")
p x |A(x)-A(p)| M|x-p| holds
----------------------------------------------------
0.5 0.6 0.030333334 0.4000 True
1.0 1.2 0.242666674 0.8000 True
1.9 2.0 0.380333345 0.4000 True
0.0 2.0 2.666666750 8.0000 True
네 줄 모두 성립한다. 마지막 줄은 구간 전체이고 여유가 가장 크다. 부등식이 팽팽하지 않은 것은 \(M\)이 최댓값이라 대부분의 자리에서 넉넉하기 때문이다.
3. 적분한 뒤 미분한다#
이 절이 이 노트북의 요점이다.
1절에서 만든 두 함수를 겹쳐 쓴다. A는 lower_upper로 만들어진 함수이고, 그것을 derivative_by_definition에 넣는다.
# Feed the integral-built function into the derivative-built procedure
print(f"{'p':>6} {'derivative of A':>18} {'f(p) = p**2':>14} {'difference':>14}")
print("-" * 56)
for p in (0.5, 1.0, 1.5, 1.9):
slope = derivative_by_definition(A, p, step=1e-4)
print(f"{p:>6} {slope:>18.9f} {p ** 2:>14.4f} {abs(slope - p ** 2):>14.2e}")
p derivative of A f(p) = p**2 difference
--------------------------------------------------------
0.5 0.250000011 0.2500 1.11e-08
1.0 1.000000035 1.0000 3.46e-08
1.5 2.250000074 2.2500 7.36e-08
1.9 3.610000116 3.6100 1.16e-07
네 줄 모두 맞는다. 왼쪽 열은 계단함수의 넓이를 더해 얻은 값들의 차분몫이고 오른쪽 열은 그냥 \(p^2\)이다. 두 절차 사이에 연결이 없는데 값이 같다.
정리 63이 주장하는 것이 이것이다.
수치가 증명은 아니다
차이 열이 \(10^{-7}\) 수준인 것은 두 가지 근사가 겹친 결과이다. \(A\)의 각 값이 \(n=4000\)의 상합과 하합을 평균한 것이고, 차분몫은 step=1e-4를 넣은 것이다.
16회 실험 2절에서 보았듯 step을 줄이면 어느 지점부터 오차가 커진다. 여기서는 \(A\)의 값 자체에도 오차가 있으므로 그 지점이 더 일찍 온다.
코드가 보이는 것은 정리가 참일 때 무엇이 보이는지이다. \(\epsilon\)과 \(\delta\)를 다루는 것은 서술 파트의 증명이다.
불연속점에서는 무너진다#
서술 파트의 흔한 오해 상자가 든 예를 만든다. \([-1,1]\)에서 \(t<0\)이면 \(0\)이고 \(t \ge 0\)이면 \(1\)인 함수이다.
이 함수는 계단함수이므로 누적 함수를 7회 (7.1)로 정확히 적을 수 있다. 근사할 필요가 없다.
# A step function and its accumulation, both written exactly
def jump(t):
return 0.0 if t < 0.0 else 1.0
def jump_accumulation(x):
# the area from -1 to x: zero while x < 0, then x
return 0.0 if x < 0.0 else x
print(f"{'x':>8} {'B(x)':>10}")
print("-" * 20)
for x in (-0.5, -0.1, 0.0, 0.1, 0.5):
print(f"{x:>8} {jump_accumulation(x):>10.4f}")
x B(x)
--------------------
-0.5 0.0000
-0.1 0.0000
0.0 0.0000
0.1 0.1000
0.5 0.5000
\(B\)는 \(0\)에서 연속이다. 정리 62가 보장한 대로이다. 그런데 미분 가능하지 않다.
# One-sided difference quotients at 0 disagree
for step in (1e-3, 1e-6):
right = (jump_accumulation(step) - jump_accumulation(0.0)) / step
left = (jump_accumulation(-step) - jump_accumulation(0.0)) / (-step)
print(f"step = {step:<8} right = {right:.6f} left = {left:.6f}")
step = 0.001 right = 1.000000 left = -0.000000
step = 1e-06 right = 1.000000 left = -0.000000
우도함수값이 \(1\)이고 좌도함수값이 \(0\)이다. \(f\)의 불연속이 \(A\)의 미분 불가능으로 옮겨 갔다.
\(f(0)=1\)이므로 정리 63이 성립했다면 \(B'(0)=1\)이어야 한다. 우도함수값만 그렇고 좌도함수값은 아니다. 정리가 \(p\)에서의 연속을 요구한 이유가 여기 있다.
오른쪽 그림에서 \(B\)의 꺾인 자리와 계단함수가 뛰는 자리가 같은 \(x\)이다.
4. 미분한 뒤 적분한다#
7회 6절이 \(\int_0^b x^2\)을 정의에서 계산할 때 쓴 식을 sympy로 다시 만든다.
# The lower and upper sums of lecture 7, section 6, rebuilt symbolically
b, n, k = sp.symbols("b n k", positive=True)
n = sp.Symbol("n", positive=True, integer=True)
k = sp.Symbol("k", integer=True)
L_n = sp.simplify(b ** 3 / n ** 3 * sp.summation(k ** 2, (k, 0, n - 1)))
U_n = sp.simplify(b ** 3 / n ** 3 * sp.summation(k ** 2, (k, 1, n)))
print("L_n =", sp.factor(L_n))
print("U_n =", sp.factor(U_n))
print("U_n - L_n =", sp.simplify(U_n - L_n))
L_n = b**3*(n - 1)*(2*n - 1)/(6*n**2)
U_n = b**3*(n + 1)*(2*n + 1)/(6*n**2)
U_n - L_n = b**3/n
\(U_n - L_n = b^3/n\)이다. 7회가 이 값이 \(0\)으로 갈 수 있음을 5회 따름정리 4로 보였다.
\(b=1\), \(n=10\)을 넣으면 7회가 적은 값이 나온다.
# The two numbers lecture 7 reports at b = 1, n = 10
print("L_10 =", L_n.subs({b: 1, n: 10}), "=", float(L_n.subs({b: 1, n: 10})))
print("U_10 =", U_n.subs({b: 1, n: 10}), "=", float(U_n.subs({b: 1, n: 10})))
print("the two trap 1/3 =", float(sp.Rational(1, 3)))
L_10 = 57/200 = 0.285
U_10 = 77/200 = 0.385
the two trap 1/3 = 0.3333333333333333
이제 정리 64로 같은 값을 구한다.
# The same value from a primitive: find F with F' = f, then subtract
t = sp.Symbol("t", real=True)
F = t ** 3 / 3
print("F =", F)
print("F' =", sp.diff(F, t), " equals f :", sp.diff(F, t) == t ** 2)
print("F(b) - F(0) =", sp.simplify(F.subs(t, b) - F.subs(t, 0)))
F = t**3/3
F' = t**2 equals f : True
F(b) - F(0) = b**3/3
세 줄이다. 도함수가 \(x^2\)인 함수를 하나 적고, 두 끝에서 값을 빼면 끝난다.
다른 함수로도 해 본다. 7회는 \(x^2\) 하나에 한 절을 썼는데, 정리 64가 있으면 여러 개를 한꺼번에 다룰 수 있다.
# Several integrals at once, each needing only a primitive
print(f"{'f':>14} {'primitive F':>16} {'integral 0..2':>16} {'by definition':>16}")
print("-" * 68)
for expr, primitive in ((t ** 2, t ** 3 / 3),
(t ** 3, t ** 4 / 4),
(3 * t ** 2 + 2 * t, t ** 3 + t ** 2),
(t ** sp.Rational(1, 2), sp.Rational(2, 3) * t ** sp.Rational(3, 2))):
exact = sp.simplify(primitive.subs(t, 2) - primitive.subs(t, 0))
numeric = integral_by_definition(sp.lambdify(t, expr, "numpy"), 0.0, 2.0)
print(f"{str(expr):>14} {str(primitive):>16} {str(exact):>16} {numeric:>16.6f}")
f primitive F integral 0..2 by definition
--------------------------------------------------------------------
t**2 t**3/3 8/3 2.666667
t**3 t**4/4 4 4.000000
3*t**2 + 2*t t**3 + t**2 12 12.000000
sqrt(t) 2*t**(3/2)/3 4*sqrt(2)/3 1.885616
네 번째 줄이 17회에서 만든 거듭제곱근이다. 정리 57가 준 공식 덕분에 원시함수를 적을 수 있다.
마지막 열의 수치 적분은 정의로 계산한 것이고, 셋째 열은 원시함수로 계산한 것이다. 두 열이 맞는다.
원시함수를 못 찾으면#
서술 파트 6절 끝이 \(1/x\)를 들었다. 원시함수를 적을 수 없으니 셋째 열을 채울 수 없다.
# 1/x has no primitive among the functions this book has built so far
reciprocal = lambda v: 1.0 / v
print(f"{'x':>6} {'integral from 1 to x':>22}")
print("-" * 30)
for x in (1.0, 2.0, 2.718281828, 4.0, 8.0):
print(f"{x:>6} {integral_by_definition(reciprocal, 1.0, x):>22.9f}")
x integral from 1 to x
------------------------------
1.0 0.000000000
2.0 0.693147184
2.718281828 1.000000013
4.0 1.386294405
8.0 2.079441793
넷째 줄의 값이 \(1\)에 가깝다. 그 자리의 \(x\)가 무엇인지, 그리고 이 표의 값들이 무슨 함수인지는 다음 회차의 일이다.
마치며#
정의 15과 정의 30를 각각 코드로 옮겼다. 한쪽은 분할과 최솟값·최댓값과 합이고 다른 쪽은 뺄셈과 나눗셈이다. 두 함수는 계산을 하나도 공유하지 않는다.
그런데 3절에서 앞의 것으로 만든 \(A\)를 뒤의 것에 넣었더니 네 점 모두에서 \(p^2\)이 나왔다. 정리 63이 주장하는 것이 그것이다.
\(f\)가 도약하는 자리에서는 무너졌다. \(B\)는 연속인데 \(0\)에서 좌우 차분몫이 \(0\)과 \(1\)로 갈렸다. 정리가 그 점에서의 연속을 요구한 이유가 눈에 보인다.
4절에서 7회의 \(L_n\)과 \(U_n\)을 sympy로 다시 만들어 \(U_n-L_n=b^3/n\)을 확인했고, \(b=1\), \(n=10\)에서 \(0.285\)와 \(0.385\)를 얻었다. 같은 값을 원시함수로는 세 줄에 구했다.
그리고 \(1/x\)에서는 셋째 열을 채울 수 없었다. 표의 값들은 수치로만 있고 이름이 없다.
서술 파트의 내용 |
이 노트북의 코드 |
|---|---|
두 정의가 서로를 언급하지 않는다 |
1절 |
2절 |
|
2절 |
|
3절 |
|
\(p\)에서의 연속이 필요하다 |
3절 좌우 차분몫이 \(0\)과 \(1\) |
3절 오른쪽 그림에서 \(B\)가 이어져 있는 것 |
|
7회 6절의 열 개의 식 |
4절 |
4절 |
더 해 볼 것#
2절의
integral_by_definition에서n을 \(4000\)에서 \(100\)으로 줄이고 3절의 표를 다시 뽑아 보자. 차이 열이 얼마나 커지는가.step을 함께 줄이면 나아지는가 나빠지는가.3절의 도약함수를 \(t<0\)에서 \(-1\)이고 \(t\ge0\)에서 \(1\)인 함수로 바꾸어 보자. \(B\)의 좌우 차분몫이 각각 얼마가 되는가. \(f(0)=1\)과 어느 쪽이 맞는가.
4절의 표에 \(f(t)=t^{-2}\)를 넣어 보자. 원시함수를 17회 공식으로 적을 수 있는가. 구간을 \([1,2]\)로 바꾸어야 하는 이유는 무엇인가.
4절 마지막 표에서 값이 정확히 \(1\)이 되는 \(x\)를 이분법으로 찾아보자. 그 수가 무엇인지 짐작이 가는가.
짝이 되는 서술