Numerical Differentiation and Integration Numerical Differentiation f ′ ( x ) = f ( x + h ) − f ( x ) h − h 2 f ′ ′ ( c )
f^{\prime}(x) = \frac{f(x+h) - f(x)}{h} - \frac{h}{2}f^{\prime\prime}(c)
f ′ ( x ) = h f ( x + h ) − f ( x ) − 2 h f ′′ ( c ) treating the last term h 2 f ′ ′ ( c ) \frac{h}{2}f^{\prime\prime}(c) 2 h f ′′ ( c ) x x x x + h x+h x + h
if error is O ( h n ) O(h^n) O ( h n ) n n n
A second-order formula can be developed if f f f
according to Taylor's Theorem
f ( x + h ) = f ( x ) + h f ′ ( x ) + h 2 2 f ′ ′ ( x ) + h 3 6 f ′ ′ ′ ( c 1 ) f ( x − h ) = f ( x ) − h f ′ ( x ) + h 2 2 f ′ ′ ( x ) − h 3 6 f ′ ′ ′ ( c 1 ) f ′ ( x ) = f ( x + h ) − f ( x − h ) 2 h − h 2 12 f ′ ′ ′ ( c 1 ) − h 2 12 f ′ ′ ′ ( c 2 )
f(x+h) = f(x) + hf'(x) + \frac{h^2}{2}f''(x) + \frac{h^3}{6}f'''(c_1) \\\\
f(x-h) = f(x) - hf'(x) + \frac{h^2}{2}f''(x) - \frac{h^3}{6}f'''(c_1) \\\\
f'(x) = \frac{f(x+h) - f(x-h)}{2h} - \frac{h^2}{12}f'''(c_1) - \frac{h^2}{12}f'''(c_2)
f ( x + h ) = f ( x ) + h f ′ ( x ) + 2 h 2 f ′′ ( x ) + 6 h 3 f ′′′ ( c 1 ) f ( x − h ) = f ( x ) − h f ′ ( x ) + 2 h 2 f ′′ ( x ) − 6 h 3 f ′′′ ( c 1 ) f ′ ( x ) = 2 h f ( x + h ) − f ( x − h ) − 12 h 2 f ′′′ ( c 1 ) − 12 h 2 f ′′′ ( c 2 ) where x − h < c 2 < x < c 1 < x + h x-h < c_2 < x < c_1 < x+h x − h < c 2 < x < c 1 < x + h
we can use Generalized Intermediate Value Theorem to replace − h 2 12 f ′ ′ ′ ( c 1 ) − h 2 12 f ′ ′ ′ ( c 2 ) -\frac{h^2}{12}f'''(c_1) - \frac{h^2}{12}f'''(c_2) − 12 h 2 f ′′′ ( c 1 ) − 12 h 2 f ′′′ ( c 2 ) − h 2 6 f ′ ′ ′ ( c ) -\frac{h^2}{6}f'''(c) − 6 h 2 f ′′′ ( c )
a 1 f ( x i ) + ⋯ + a n f ( x i ) ≤ a 1 f ( x 1 ) + ⋯ + a n f ( x n ) ≤ a 1 f ( x j ) + ⋯ + a n f ( x j ) f ( x i ) ≤ a 1 f ( x 1 ) + ⋯ + a n f ( x n ) a 1 + ⋯ + a n ≤ f ( x j )
\begin{align}
&
a_1f(x_i) + \cdots + a_nf(x_i) \le
a_1f(x_1) + \cdots + a_nf(x_n) \le
a_1f(x_j) + \cdots + a_nf(x_j)
\\\\
&
f(x_i) \le \frac{a_1f(x_1) + \cdots + a_nf(x_n)} {a_1 + \cdots + a_n} \le f(x_j)
\end{align}
a 1 f ( x i ) + ⋯ + a n f ( x i ) ≤ a 1 f ( x 1 ) + ⋯ + a n f ( x n ) ≤ a 1 f ( x j ) + ⋯ + a n f ( x j ) f ( x i ) ≤ a 1 + ⋯ + a n a 1 f ( x 1 ) + ⋯ + a n f ( x n ) ≤ f ( x j ) by the Intermediate Value Theorem, there is a number c c c x i x_i x i x j x_j x j
f ( c ) = a 1 f ( x 1 ) + ⋯ + a n f ( x n ) a 1 + ⋯ + a n
f(c) = \frac{a_1f(x_1) + \cdots + a_nf(x_n)} {a_1 + \cdots + a_n}
f ( c ) = a 1 + ⋯ + a n a 1 f ( x 1 ) + ⋯ + a n f ( x n ) f ′ ( x ) = f ( x + h ) − f ( x − h ) 2 h − h 2 6 f ′ ′ ′ ( c )
f'(x) = \frac{f(x+h) - f(x-h)}{2h} - \frac{h^2}{6}f'''(c)
f ′ ( x ) = 2 h f ( x + h ) − f ( x − h ) − 6 h 2 f ′′′ ( c ) f ( x + h ) = f ( x ) + h f ′ ( x ) + h 2 2 f ′ ′ ( x ) + h 3 6 f ′ ′ ′ ( x ) + h 4 24 f ( i v ) ( c 1 ) f ( x − h ) = f ( x ) − h f ′ ( x ) + h 2 2 f ′ ′ ( x ) − h 3 6 f ′ ′ ′ ( x ) + h 4 24 f ( i v ) ( c 1 ) f ′ ′ ( x ) = f ( x − h ) − 2 f ( x ) + f ( x + h ) h 2 − h 2 12 f ( i v ) ( c )
\begin{align}
&
f(x+h) = f(x) + hf'(x) +\frac{h^2}{2} f''(x) + \frac{h^3}{6}f'''(x) + \frac{h^4}{24}f^{(iv)}(c_1) \\\\
&
f(x-h) = f(x) - hf'(x) +\frac{h^2}{2} f''(x) - \frac{h^3}{6}f'''(x) + \frac{h^4}{24}f^{(iv)}(c_1) \\\\
&
f''(x) = \frac{f(x-h) - 2f(x) + f(x+h)}{h^2} - \frac{h^2}{12}f^{(iv)}(c)
\end{align}
f ( x + h ) = f ( x ) + h f ′ ( x ) + 2 h 2 f ′′ ( x ) + 6 h 3 f ′′′ ( x ) + 24 h 4 f ( i v ) ( c 1 ) f ( x − h ) = f ( x ) − h f ′ ( x ) + 2 h 2 f ′′ ( x ) − 6 h 3 f ′′′ ( x ) + 24 h 4 f ( i v ) ( c 1 ) f ′′ ( x ) = h 2 f ( x − h ) − 2 f ( x ) + f ( x + h ) − 12 h 2 f ( i v ) ( c ) where c between x − h x-h x − h x + h x+h x + h
Rounding error denote the floating point version of the input f ( x + h ) f(x+h) f ( x + h ) f ^ ( x + h ) \hat{f}(x+h) f ^ ( x + h )
f ^ ( x + h ) = f ( x + h ) + ϵ 1 f ^ ( x − h ) = f ( x − h ) + ϵ 2
\hat{f}(x+h) = f(x+h) + \epsilon_1 \\
\hat{f}(x-h) = f(x-h) + \epsilon_2 \\
f ^ ( x + h ) = f ( x + h ) + ϵ 1 f ^ ( x − h ) = f ( x − h ) + ϵ 2 where ∣ ϵ 1 ∣ , ∣ ϵ 2 ∣ ≈ ϵ m a c h |\epsilon_1|,|\epsilon_2|\approx \epsilon_{mach} ∣ ϵ 1 ∣ , ∣ ϵ 2 ∣ ≈ ϵ ma c h
f ′ ( x ) c o r r e c t − f ′ ( x ) m a c h i n e = f ′ ( x ) − f ^ ( x + h ) − f ^ x − h 2 h = f ′ ( x ) − f ( x + h ) + ϵ 1 − ( f ( x − h ) + ϵ 2 ) 2 h = ( f ′ ( x ) − f ( x + h ) − f ( x − h ) 2 h ) + ϵ 2 − ϵ 1 2 h = ( f ′ ( x ) c o r r e c t − f ′ ( x ) f o r m u l a ) + e r r o r r o u n d i n g
\begin{align}
f'(x)_{correct} - f'(x)_{machine} &=
f'(x) - \frac{\hat{f}(x+h) - \hat{f}{x-h}}{2h} \\
&=f'(x) - \frac{f(x+h) + \epsilon_1 - (f(x-h) + \epsilon_2)}{2h} \\
&=(f'(x) - \frac{f(x+h) - f(x-h)}{2h}) + \frac{\epsilon_2 - \epsilon_1}{2h} \\
&=(f'(x)_{correct} - f'(x)_{formula}) + error_{rounding}
\end{align}
f ′ ( x ) correc t − f ′ ( x ) ma c hin e = f ′ ( x ) − 2 h f ^ ( x + h ) − f ^ x − h = f ′ ( x ) − 2 h f ( x + h ) + ϵ 1 − ( f ( x − h ) + ϵ 2 ) = ( f ′ ( x ) − 2 h f ( x + h ) − f ( x − h ) ) + 2 h ϵ 2 − ϵ 1 = ( f ′ ( x ) correc t − f ′ ( x ) f or m u l a ) + erro r ro u n d in g ∣ ϵ 2 − ϵ 1 2 h ∣ ≤ 2 ϵ m a c h 2 h = ϵ m a c h h
\left|\frac{\epsilon_2 - \epsilon_1}{2h}\right| \le \frac{2\epsilon_{mach}}{2h} = \frac{\epsilon_{mach}}{h}
2 h ϵ 2 − ϵ 1 ≤ 2 h 2 ϵ ma c h = h ϵ ma c h E ( h ) ≡ h 2 6 f ′ ′ ′ ( c ) + ϵ m a c h h
E(h) \equiv \frac{h^2}{6}f'''(c) + \frac{\epsilon_{mach}}{h}
E ( h ) ≡ 6 h 2 f ′′′ ( c ) + h ϵ ma c h where x − h < c < x + h x-h < c < x+h x − h < c < x + h
the minimum of E ( h ) E(h) E ( h ) 0 = E ′ ( h ) 0 = E'(h) 0 = E ′ ( h )
Assume that we are presented with an order n n n F ( h ) F(h) F ( h )
Q ≈ F ( h ) + K h n
Q\approx F(h) + Kh^n
Q ≈ F ( h ) + K h n where K K K h h h
f ′ ( x ) = f ( x + h ) − f ( x − h ) 2 h − f ′ ′ ′ ( c h ) 6 h 2
f'(x) = \frac{f(x+h) - f(x-h)}{2h} - \frac{f'''(c_h)}{6}h^2
f ′ ( x ) = 2 h f ( x + h ) − f ( x − h ) − 6 f ′′′ ( c h ) h 2 the point c h c_h c h x x x x + h x+h x + h
F ( h ) + K ( h ) n − F ( h / 2 ) ≈ 1 2 n ( F ( h ) + K h n − F ( h ) ) Q − F ( h / 2 ) ≈ 1 2 n ( Q − F ( h ) )
\begin{align}
F(h) + K(h)^n - F(h/2) &\approx \frac{1}{2^n}(F(h) + Kh^n - F(h))\\
Q - F(h/2) &\approx \frac{1}{2^n}(Q - F(h))\\
\end{align}
F ( h ) + K ( h ) n − F ( h /2 ) Q − F ( h /2 ) ≈ 2 n 1 ( F ( h ) + K h n − F ( h )) ≈ 2 n 1 ( Q − F ( h )) Q ≈ 2 n F ( h / 2 ) − F ( h ) 2 n − 1
Q \approx \frac{ 2^nF(h/2) - F(h)}{2^n-1}
Q ≈ 2 n − 1 2 n F ( h /2 ) − F ( h ) use this extrapolation formula (sometimes called Richardson extrapolation ) to leverage an order n n n
(a)( 0 , 0 ) (0,0) ( 0 , 0 ) ( h , 1 ) (h,1) ( h , 1 )
the area is
∫ 0 h x h d x = h / 2
\int_0^h\frac{x}{h}dx = h/2
∫ 0 h h x d x = h /2 (b)P ( x ) P(x) P ( x ) ( − h , 0 ) , ( 0 , 1 ) (-h,0),(0,1) ( − h , 0 ) , ( 0 , 1 ) ( h , 0 ) (h,0) ( h , 0 )
∫ − h h x 2 − h 2 h 2 d x = 4 3 h
\int_{-h}^{h} \frac{x^2-h^2}{h^2}dx = \frac{4}{3}h
∫ − h h h 2 x 2 − h 2 d x = 3 4 h (c)( − h , 1 ) , ( 0 , 0 ) (-h,1),(0,0) ( − h , 1 ) , ( 0 , 0 ) ( h , 0 ) (h,0) ( h , 0 )
∫ − h 0 x 2 − x h 2 h 2 + ∫ 0 h x 2 − x h 2 h 2 = 1 3 h
\int_{-h}^{0}\frac{x^2-xh}{2h^2} +
\int_{0}^{h}\frac{x^2-xh}{2h^2} = \frac{1}{3}h
∫ − h 0 2 h 2 x 2 − x h + ∫ 0 h 2 h 2 x 2 − x h = 3 1 h Trapezoid Rule ∫ x 0 x 1 f ( x ) d x = h 2 ( y 0 + y 1 ) − h 3 12 f ′ ′ ( c )
\int_{x_0}^{x_1} f(x)dx = \frac{h}{2}(y_0+y_1) - \frac{h^3}{12}f''(c)
∫ x 0 x 1 f ( x ) d x = 2 h ( y 0 + y 1 ) − 12 h 3 f ′′ ( c ) where h = x 1 − x 0 h = x_1 - x_0 h = x 1 − x 0 c c c x 0 x_0 x 0 x 1 x_1 x 1
details
consider the degree 1 interpolating polynomial P 1 ( x ) P_1(x) P 1 ( x ) ( x 0 , y 0 ) (x_0,y_0) ( x 0 , y 0 ) ( x 1 , y 1 ) (x_1,y_1) ( x 1 , y 1 )
f ( x ) = y 0 x − x 1 x 0 − x 1 + y 1 x − x 0 x 1 − x 0 + ( x − x 0 ) ( x − x 1 ) 2 ! f ′ ′ ( c x ) = P ( x ) + E ( x ) ∫ x 0 x 1 f ( x ) d x = ∫ x 0 x 1 P ( x ) d x + ∫ x 0 x 1 E ( x ) d x ∫ x 0 x 1 P ( x ) d x = y 0 ∫ x 0 x 1 x − x 1 x 0 − x 1 d x + y 1 ∫ x 0 x 1 x − x 0 x 1 − x 0 d x = y 0 x 1 − x 0 2 + y 1 x 1 − x 0 2 = ( x 1 − x 0 ) y 0 + y 1 2 ∫ x 0 x 1 E ( x ) d x = − ( x 1 − x 0 ) 3 12 f ′ ′ ( c )
\begin{align}
&
f(x) =
y_0 \frac{x-x_1}{x_0 -x_1} +
y_1 \frac{x-x_0}{x_1 -x_0} +
\frac{(x-x_0)(x-x_1)}{2!}f''(c_x) = P(x) + E(x) \\\\
&
\int_{x_0}^{x_1} f(x)dx = \int_{x_0}^{x_1}P(x) dx + \int_{x_0}^{x_1} E(x)dx \\\\
&
\begin{align}
\int_{x_0}^{x_1}P(x)dx &= y_0\int_{x_0}^{x_1}\frac{x-x_1}{x_0-x_1}dx + y_1\int_{x_0}^{x_1}\frac{x-x_0}{x_1-x_0}dx\\
&=y_0\frac{x_1 - x_0}{2} + y_1\frac{x_1 - x_0}{2} = (x_1 - x_0)\frac{y_0+y_1}{2}
\end{align}
\\\\
&
\int_{x_0}^{x_1} E(x)dx = -\frac{(x_1-x_0)^3}{12} f''(c)
\end{align}
f ( x ) = y 0 x 0 − x 1 x − x 1 + y 1 x 1 − x 0 x − x 0 + 2 ! ( x − x 0 ) ( x − x 1 ) f ′′ ( c x ) = P ( x ) + E ( x ) ∫ x 0 x 1 f ( x ) d x = ∫ x 0 x 1 P ( x ) d x + ∫ x 0 x 1 E ( x ) d x ∫ x 0 x 1 P ( x ) d x = y 0 ∫ x 0 x 1 x 0 − x 1 x − x 1 d x + y 1 ∫ x 0 x 1 x 1 − x 0 x − x 0 d x = y 0 2 x 1 − x 0 + y 1 2 x 1 − x 0 = ( x 1 − x 0 ) 2 y 0 + y 1 ∫ x 0 x 1 E ( x ) d x = − 12 ( x 1 − x 0 ) 3 f ′′ ( c ) Simpson's Rule ∫ x 0 x 2 f ( x ) d x = h 3 ( y 0 + 4 y 1 + y 2 ) − h 5 90 f ( i v ) ( c )
\int_{x_0} ^{x_2} f(x) dx = \frac{h}{3}(y_0+4y_1+y_2) - \frac{h^5}{90}f^{(iv)}(c)
∫ x 0 x 2 f ( x ) d x = 3 h ( y 0 + 4 y 1 + y 2 ) − 90 h 5 f ( i v ) ( c ) where h = x 2 − x 1 = x 1 − x 0 h = x_2 - x_1 = x_1 - x_0 h = x 2 − x 1 = x 1 − x 0 c c c x 0 x_0 x 0 x 2 x_2 x 2
details
simpson's rule is similar to trapezoid rule , except that the degree 1 interpolant is replaced by a parabola
Simpson's 3/8 Rule ∫ x 0 x 3 f ( x ) d x ≈ 3 h 8 ( y 0 + 3 y 1 + 3 y 2 + y 3 )
\int_{x_0}^{x_3} f(x)dx \approx \frac{3h}{8}(y_0 + 3y_1 + 3y_2 + y_3)
∫ x 0 x 3 f ( x ) d x ≈ 8 3 h ( y 0 + 3 y 1 + 3 y 2 + y 3 ) dividing the interval up into several subintervals, applying the rule separately on each one, and then totaling up . The strategy is called composite numerical integration
Composite Trapezoid Rule ∫ a b f ( x ) d x = h 2 [ f ( a ) + f ( b ) + 2 ∑ i = 1 m − 1 f ( x i ) ] − ∑ i = 0 m − 1 h 3 12 f ′ ′ ( c i )
\int_a^b f(x)dx = \frac{h}{2}\left[f(a) + f(b) + 2\sum_{i=1}^{m-1}f(x_i)\right] -\sum_{i=0}^{m-1}\frac{h^3}{12}f''(c_i)
∫ a b f ( x ) d x = 2 h [ f ( a ) + f ( b ) + 2 i = 1 ∑ m − 1 f ( x i ) ] − i = 0 ∑ m − 1 12 h 3 f ′′ ( c i ) the error term can be written
h 3 12 ∑ i = 1 m − 1 f ′ ′ ( c i ) = h 3 12 m f ′ ′ ( c )
\frac{h^3}{12}\sum_{i=1}^{m-1}f''(c_i) = \frac{h^3}{12}mf''(c)
12 h 3 i = 1 ∑ m − 1 f ′′ ( c i ) = 12 h 3 m f ′′ ( c ) where m h = ( b − a ) mh = (b-a) mh = ( b − a ) a < c < b a<c<b a < c < b f ′ ′ f'' f ′′ [ a , b ] [a,b] [ a , b ]
∫ a b f ( x ) d x = h 2 [ y 0 + y m + 2 ∑ i = 1 m − 1 y i ] − ( b − a ) h 2 12 f ′ ′ ( c )
\int_a^b f(x)dx = \frac{h}{2}\left[y_0 + y_m + 2\sum_{i=1}^{m-1}y_i\right] -\frac{(b-a)h^2}{12}f''(c)
∫ a b f ( x ) d x = 2 h [ y 0 + y m + 2 i = 1 ∑ m − 1 y i ] − 12 ( b − a ) h 2 f ′′ ( c ) Composite Simpson's Rule ∫ a b f ( x ) d x = h 3 [ y 0 + y 2 m + 4 ∑ i = 1 m y 2 i − 1 + 2 ∑ i = 1 m − 1 y 2 i ] − ( b − a ) h 4 180 f ( i v ) ( c )
\int_a^b f(x)dx = \frac{h}{3}\left[y_0 + y_{2m} + 4\sum_{i=1}^{m}y_{2i-1}+2\sum_{i=1}^{m-1}y_{2i}\right] - \frac{(b-a)h^4}{180}f^{(iv)}(c)
∫ a b f ( x ) d x = 3 h [ y 0 + y 2 m + 4 i = 1 ∑ m y 2 i − 1 + 2 i = 1 ∑ m − 1 y 2 i ] − 180 ( b − a ) h 4 f ( i v ) ( c ) where c c c a a a b b b
∫ x 2 i x 2 i + 2 f ( x ) d x = h 3 [ f ( x 2 i ) + 4 f ( x 2 i + 1 ) + f ( x 2 i + 2 ) ] − h 5 90 f ( i v ) ( c i )
\int_{x_{2i}}^{x_{2i+2}}f(x)dx = \frac{h}{3}[f(x_{2i}) + 4f(x_{2i+1}) + f(x_{2i+2})] - \frac{h^5}{90}f^{(iv)}(c_i)
∫ x 2 i x 2 i + 2 f ( x ) d x = 3 h [ f ( x 2 i ) + 4 f ( x 2 i + 1 ) + f ( x 2 i + 2 )] − 90 h 5 f ( i v ) ( c i ) 
