함수의 정확한 형태가 주어지지 않고, 몇몇 점들의 함수값들만 알고 있다면 미분을 어떻게 근사할 수 있을까?
이런 문제에 접근하는 근사 방법들을 일컬어 수치 미분Numerical Differentiation이라고 부릅니다.
가장 쉬운 근사 가장 간단한 근사로 미분의 정의를 활용해볼 수 있습니다. 함수 f f f x 0 x_0 x 0
f ′ ( x 0 ) = lim h → 0 f ( x 0 + h ) − f ( x 0 ) h
f'(x_0) = \lim _{h \to 0 } \frac {f(x_0 + h) - f(x_0)} h
f ′ ( x 0 ) = h → 0 lim h f ( x 0 + h ) − f ( x 0 )
따라서 $ \dfrac {f(x_0 + h) - f(x_0)} h를 작은 $h f ′ ( x 0 ) f'(x_0) f ′ ( x 0 )
수치적 안정성의 문제 하지만 컴퓨터한테 위 근사를 시키면 수를 표현할 때 사용하는 비트 수에 한계가 있어 반올림을 하게 되는데 이 때 발생하는 오차를 반올림 오차round-off-error라고 합니다.
f ( x 0 + h ) − f ( x 0 ) = δ + c p ( f ( x 0 + h ) − f ( x 0 ) )
f(x_0 + h) - f(x_0) =\delta+ cp(f(x_0+h) - f(x_0))
f ( x 0 + h ) − f ( x 0 ) = δ + c p ( f ( x 0 + h ) − f ( x 0 ) )
계산할 때 분자 부분의 반올림 오차를 위와 같이 δ \delta δ c p ( e x p r ) cp(expr) c p ( e x p r )
이제 양변을 h로 나누면
f ( x 0 + h ) − f ( x 0 ) h = δ h + c p ( f ( x 0 + h ) − f ( x 0 ) ) h
\frac{f(x_0 + h) - f(x_0)} h =\frac \delta h+ \frac {cp(f(x_0+h) - f(x_0))} h
h f ( x 0 + h ) − f ( x 0 ) = h δ + h c p ( f ( x 0 + h ) − f ( x 0 ) )
이기 때문에 전체 오차는 δ / h \delta/h δ / h h h h δ / h \delta / h δ / h
넘어가기 전에 라그랑주 보간법Lagrange interpolation은 함수 f를 n개의 점을 지나는 다항식으로 근사할 수 있습니다. 구간 I 위에서 정의된 f ∈ C n + 1 ( I ) f \in C^{n+1} (I) f ∈ C n + 1 ( I ) x 0 , x 1 , ⋯ , x n x_0, x_1, \cdots, x_n x 0 , x 1 , ⋯ , x n f ( x 0 ) , f ( x 1 ) , ⋯ , f ( x n ) f(x_0), f(x_1), \cdots, f(x_n) f ( x 0 ) , f ( x 1 ) , ⋯ , f ( x n )
f ( x ) = ∑ k = 0 n f ( x k ) L k ( x ) + f ( n + 1 ) ( ξ ( x ) ) ( n + 1 ) ! ( x − x 0 ) ⋯ ( x − x n ) where L k ( x ) = ∏ i = 0 i ≠ k n x − x i x k − x i
f(x) = \sum _{k=0} ^n f(x_k) L_k(x) + \frac {f^{(n+1)}(\xi(x))}{(n+1)!} (x-x_0)\cdots(x-x_n)\\
\text{where } L_k(x) = \prod _{i=0\\i \ne k} ^{n} \frac {x-x_i}{x_k-x_i}
f ( x ) = k = 0 ∑ n f ( x k ) L k ( x ) + ( n + 1 ) ! f ( n + 1 ) ( ξ ( x ) ) ( x − x 0 ) ⋯ ( x − x n ) where L k ( x ) = i = 0 i = k ∏ n x k − x i x − x i
라는 정리가 알려져 있습니다. 이 때, ξ ( x ) \xi(x) ξ ( x ) I I I
f ∈ C 2 [ a , b ] f \in C^2 [a,b] f ∈ C 2 [ a , b ] x 0 , x 1 = x 0 + h x_0, x_1=x_0 + h x 0 , x 1 = x 0 + h
f ( x ) = f ( x 0 ) x − x 0 − h − h + f ( x 1 ) x − x 0 h + f ′ ′ ( ξ ( x ) ) 2 ( x − x 0 ) ( x − x 0 − h )
f(x) = f(x_0) \frac {x-x_0-h}{-h} + f(x_1) \frac {x-x_0} h + \frac {f''(\xi(x))} 2 (x-x_0) (x-x_0 - h)
f ( x ) = f ( x 0 ) − h x − x 0 − h + f ( x 1 ) h x − x 0 + 2 f ′ ′ ( ξ ( x ) ) ( x − x 0 ) ( x − x 0 − h )
여기서 한 번 더 미분하면
f ′ ( x ) = − f ( x 0 ) h + f ( x 1 ) h + ( 2 ( x − x 0 ) − h ) 2 f ′ ′ ( ξ ( x ) ) + ( x − x 0 ) ( x − x 0 − h ) 2 d d x f ′ ′ ( ξ ( x ) )
f'(x) = - \frac {f(x_0)} h + \frac {f(x_1)} h + \frac {(2(x-x_0) - h)} 2 f'' (\xi(x)) + \frac { (x-x_0)(x-x_0-h)} 2 \frac {d}{dx} f''(\xi(x))
f ′ ( x ) = − h f ( x 0 ) + h f ( x 1 ) + 2 ( 2 ( x − x 0 ) − h ) f ′ ′ ( ξ ( x ) ) + 2 ( x − x 0 ) ( x − x 0 − h ) d x d f ′ ′ ( ξ ( x ) )
가 되고, x = x 0 x= x_0 x = x 0
f ′ ( x 0 ) = f ( x 0 + h ) − f ( x 0 ) h − h 2 f ′ ′ ( ξ ( x 0 ) )
f'(x_0) = \frac { f(x_0 + h) - f(x_0)} h - \frac h 2 f''(\xi(x_0))
f ′ ( x 0 ) = h f ( x 0 + h ) − f ( x 0 ) − 2 h f ′ ′ ( ξ ( x 0 ) )
를 얻습니다. 이는 처음에 미분 정의로 근사했던 게 점 두개로 라그랑주 보간법을 사용한 결과와 상통합니다. 심지어
∣ f ′ ( x 0 ) − f ( x 0 + h ) − f ( x 0 ) h ∣ = ∣ h 2 f ′ ′ ( ξ ( x 0 ) ) ∣ = ∣ h ∣ 2 max a ≤ x ≤ b ∣ f ′ ′ ( x ) ∣
\left| f'(x_0) - \frac { f(x_0 + h) - f(x_0)} h \right| = \left|\frac h 2 f''(\xi(x_0))\right| = \frac {|h|} 2 \max _{a \le x \le b} |f''(x)|
∣ ∣ ∣ ∣ ∣ f ′ ( x 0 ) − h f ( x 0 + h ) − f ( x 0 ) ∣ ∣ ∣ ∣ ∣ = ∣ ∣ ∣ ∣ ∣ 2 h f ′ ′ ( ξ ( x 0 ) ) ∣ ∣ ∣ ∣ ∣ = 2 ∣ h ∣ a ≤ x ≤ b max ∣ f ′ ′ ( x ) ∣
라는 error bound까지 얻을 수 있네요!
이 근사법에서
h > 0이면 forward-difference formula h < 0이면 backward-difference formula 라고 부릅니다.
라그랑주 보간법은 여러 개의 점을 통한 근사 방법을 제공하기 때문에 Two-Point Formulas를 넘어서 Three-Point, Four-Point, Five-Point Formulas, …로 자연스럽게 확장이 가능합니다. 점 x 0 , ⋯ , x n x_0, \cdots, x_n x 0 , ⋯ , x n x = x j x=x_j x = x j
f ′ ( x j ) = ∑ k = 0 n f ( x k ) L k ′ ( x j ) + f ( n + 1 ) ( ξ ( x j ) ) ( n + 1 ) ! ∏ k = 0 k ≠ j n ( x j − x k )
f'(x_j) = \sum_{k=0} ^n f(x_k) L_k ' (x_j) + \frac {f^{(n+1)} (\xi(x_j))}{(n+1)!} \prod _{k=0\\k\ne j}^n (x_j-x_k)
f ′ ( x j ) = k = 0 ∑ n f ( x k ) L k ′ ( x j ) + ( n + 1 ) ! f ( n + 1 ) ( ξ ( x j ) ) k = 0 k = j ∏ n ( x j − x k )
이고, x 0 , ⋯ , x n x_0, \cdots, x_n x 0 , ⋯ , x n ( x j − x k ) (x_j - x_k) ( x j − x k ) O ( h n ) O(h^n) O ( h n )
여기서 주로 사용하는 버전은 three-point와 five-point formula입니다. 여기서는 간단하게 three-point formula만 소개하려 합니다.
x 0 , x 1 = x 0 + h , x 1 = x 0 + 2 h x_0, x_1= x_0 + h, x_1 = x_0 + 2h x 0 , x 1 = x 0 + h , x 1 = x 0 + 2 h
\begin{aligned}
f'(x_j)
&= f(x_0) \frac {2x_j - 2x_0 - 3h} {2h^2}
- f(x_0 + h) \frac {2x_j - 2x_0 - 2h}{h^2}
\\&+ f(x_0 + 2h) \frac {2x_j - x_0 - h}{2h^2}
+ \frac 1 6 f^{(3)} (\xi_j) \prod _{k=0\\k\ne j } ^2 (x_j - x_k)
\end{aligned}
f(x_0 + h) \frac {2x_j - 2x_0 - 2h}{h^2}
\&+ f(x_0 + 2h) \frac {2x_j - x_0 - h}{2h^2} \frac 1 6 f^{(3)} (\xi_j) \prod _{k=0\k\ne j } ^2 (x_j - x_k)
\end{aligned}
$$ 이고, 결과적으로 다음의 결과를 얻습니다.
f ′ ( x 0 ) = 1 2 h [ − 3 f ( x 0 ) + 4 ( x 0 + h ) − f ( x 0 + 2 h ) ] + h 2 3 f ( 3 ) ( ξ 0 ) \displaystyle f'(x_0) = \frac {1} {2h} [-3f(x_0) + 4(x_0+ h) - f(x_0 + 2h)] + \frac {h^2} 3 f^{(3)}(\xi_0) f ′ ( x 0 ) = 2 h 1 [ − 3 f ( x 0 ) + 4 ( x 0 + h ) − f ( x 0 + 2 h ) ] + 3 h 2 f ( 3 ) ( ξ 0 ) f ′ ( x 1 ) = 1 2 h [ − f ( x 1 − h ) + f ( x 1 + h ) ] − h 2 6 f ( 3 ) ( ξ 1 ) \displaystyle f'(x_1) = \frac 1 {2h} [-f(x_1-h) + f(x_1 + h)] - \frac {h^2} 6 f^{(3)}(\xi_1) f ′ ( x 1 ) = 2 h 1 [ − f ( x 1 − h ) + f ( x 1 + h ) ] − 6 h 2 f ( 3 ) ( ξ 1 ) f ′ ( x 2 ) = 1 2 h [ f ( x 2 − 2 h ) − 4 f ( x 2 − h ) + 3 f ( x 2 ) ] + h 2 3 f ( 3 ) ( ξ 1 ) \displaystyle f'(x_2) = \frac 1 {2h} [f(x_2 - 2h) -4f(x_2 - h) + 3f(x_2)] + \frac {h^2} 3 f^{(3)}(\xi_1) f ′ ( x 2 ) = 2 h 1 [ f ( x 2 − 2 h ) − 4 f ( x 2 − h ) + 3 f ( x 2 ) ] + 3 h 2 f ( 3 ) ( ξ 1 ) 자세히 보면 x 1 x_1 x 1 f ( x 1 + h ) f(x_1 + h) f ( x 1 + h )
수치적 안정성 three-point midpoint formula는 다음과 같습니다.
f ′ ( x 0 ) = 1 2 h [ f ( x 0 + h ) − f ( x 0 − h ) ] − h 2 6 f ( 3 ) ( ξ ) where h > 0
f'(x_0) = \frac 1 {2h} [f(x_0 + h)-f(x_0-h)] - \frac {h^2} 6 f^{(3)}(\xi) \quad \text{ where } h > 0
f ′ ( x 0 ) = 2 h 1 [ f ( x 0 + h ) − f ( x 0 − h ) ] − 6 h 2 f ( 3 ) ( ξ ) where h > 0
이 때 각 함수값의 round-off error를 ϵ \epsilon ϵ
ϵ = max { ∣ f ( x 0 + h ) − c p ( f ( x 0 + h ) ) ∣ , ∣ f ( x 0 − h ) − c p ( f ( x 0 − h ) ) ∣ }
\epsilon = \max \{|f(x_0 + h) - cp(f(x_0 + h))|, |f(x_0-h) - cp(f(x_0 - h))| \}
ϵ = max { ∣ f ( x 0 + h ) − c p ( f ( x 0 + h ) ) ∣ , ∣ f ( x 0 − h ) − c p ( f ( x 0 − h ) ) ∣ }
그러면
e ( h ) = ∣ f ′ ( x 0 ) − 1 2 h [ c p ( f ( x 0 + h ) ) − c p ( f ( x 0 − h ) ) ] ∣ = ∣ 1 2 h [ f ( x 0 + h ) − f ( x 0 − h ) ] − h 2 6 f ( 3 ) ( ξ ) − 1 2 h [ c p ( f ( x 0 + h ) ) − c p ( f ( x 0 − h ) ) ] ∣ ≤ h 2 6 M + ϵ h where M = max x ∈ I ∣ f ( 3 ) ( x ) ∣
\begin{aligned}
e(h) &= \left|f'(x_0) - \frac {1}{2h} [cp(f(x_0 + h)) - cp(f(x_0 - h))]\right|
\\ &= \left | \frac 1 {2h} [f(x_0 + h)-f(x_0-h)] - \frac {h^2} 6 f^{(3)}(\xi) - \frac {1}{2h} [cp(f(x_0 + h)) - cp(f(x_0 - h))] \right|
\\ &\le \frac {h^2} 6 M + \frac \epsilon h\quad \text{ where } M = \max_{x \in I} {|f^{(3)}(x)|}
\end{aligned}
e ( h ) = ∣ ∣ ∣ ∣ ∣ f ′ ( x 0 ) − 2 h 1 [ c p ( f ( x 0 + h ) ) − c p ( f ( x 0 − h ) ) ] ∣ ∣ ∣ ∣ ∣ = ∣ ∣ ∣ ∣ ∣ 2 h 1 [ f ( x 0 + h ) − f ( x 0 − h ) ] − 6 h 2 f ( 3 ) ( ξ ) − 2 h 1 [ c p ( f ( x 0 + h ) ) − c p ( f ( x 0 − h ) ) ] ∣ ∣ ∣ ∣ ∣ ≤ 6 h 2 M + h ϵ where M = x ∈ I max ∣ f ( 3 ) ( x ) ∣
전항 h 2 M / 6 h^2M / 6 h 2 M / 6 ϵ / h \epsilon / h ϵ / h
truncation error과 round-off error 사이의 tradeoff
이를 최소화하는 선택은 e ( h ) e(h) e ( h ) h = 3 ϵ / M 3 h = \sqrt[3]{3\epsilon / M} h = 3 3 ϵ / M M M M h h h
마지막으로 이계 미분을 근사하는 기법에 대해서 소개하고 마치려 합니다. 3차 테일러 근사를 이용하면
f ( x 0 + h ) = f ( x 0 ) + f ′ ( x 0 ) h + 1 2 f ′ ′ ( x 0 ) h 2 + 1 6 f ′ ′ ′ ( x 0 ) h 3 + 1 24 f ( 4 ) ( ξ 1 ) h 4 f ( x 0 − h ) = f ( x 0 ) − f ′ ( x 0 ) h + 1 2 f ′ ′ ( x 0 ) h 2 − 1 6 f ′ ′ ′ ( x 0 ) h 3 + 1 24 f ( 4 ) ( ξ − 1 ) h 4 where x 0 − h < ξ − 1 < x 0 < ξ < x 0 + h
f(x_0 + h) = f(x_0) + f'(x_0) h + \frac 1 2 f''(x_0) h^2 + \frac 1 6 f'''(x_0) h^3 + \frac 1 {24} f^{(4)}(\xi_1)h^4\\
f(x_0 - h) = f(x_0) - f'(x_0) h + \frac 1 2 f''(x_0) h^2 - \frac 1 6 f'''(x_0) h^3 + \frac 1 {24} f^{(4)}(\xi_{-1})h^4\\
\text{where } x_0 - h < \xi_{-1} < x_0 < \xi < x_0 + h
f ( x 0 + h ) = f ( x 0 ) + f ′ ( x 0 ) h + 2 1 f ′ ′ ( x 0 ) h 2 + 6 1 f ′ ′ ′ ( x 0 ) h 3 + 2 4 1 f ( 4 ) ( ξ 1 ) h 4 f ( x 0 − h ) = f ( x 0 ) − f ′ ( x 0 ) h + 2 1 f ′ ′ ( x 0 ) h 2 − 6 1 f ′ ′ ′ ( x 0 ) h 3 + 2 4 1 f ( 4 ) ( ξ − 1 ) h 4 where x 0 − h < ξ − 1 < x 0 < ξ < x 0 + h
이고, 이 둘을 더하면 f ′ ( x 0 ) f'(x_0) f ′ ( x 0 )
f ′ ′ ( x 0 ) = 1 h 2 [ f ( x 0 − h ) − 2 f ( x 0 ) + f ( x 0 + h ) ] − h 2 24 [ f ( 4 ) ( ξ 1 ) + f ( 4 ) ( ξ − 1 ) ]
f''(x_0) = \frac 1 {h^2} [f(x_0-h) - 2f(x_0) + f(x_0 + h)] - \frac {h^2} {24} [f^{(4)}(\xi_1) + f^{(4)}(\xi_{-1})]
f ′ ′ ( x 0 ) = h 2 1 [ f ( x 0 − h ) − 2 f ( x 0 ) + f ( x 0 + h ) ] − 2 4 h 2 [ f ( 4 ) ( ξ 1 ) + f ( 4 ) ( ξ − 1 ) ]
인데, intermediate value theorem에 의하여
f ′ ′ ( x 0 ) = 1 h 2 [ f ( x 0 − h ) − 2 f ( x 0 ) + f ( x 0 + h ) ] − h 2 12 f ( 4 ) ( ξ )
f''(x_0) = \frac 1 {h^2} [f(x_0-h) - 2f(x_0) + f(x_0 + h)] - \frac {h^2}{12} f^{(4)}(\xi)
f ′ ′ ( x 0 ) = h 2 1 [ f ( x 0 − h ) − 2 f ( x 0 ) + f ( x 0 + h ) ] − 1 2 h 2 f ( 4 ) ( ξ )
가 도출됩니다. 이렇게 f ∈ C 4 [ x 0 − h , x 0 + h ] f \in C^4[x_0-h, x_0+h] f ∈ C 4 [ x 0 − h , x 0 + h ] O ( h 2 ) O(h^2) O ( h 2 ) f ′ ′ ( x 0 ) f \prime\prime (x_0) f ′ ′ ( x 0 )
References [1] KAIST 2021 Spring MAS365 Numerical Analysis