수치적분

사다리꼴 공식과 Lagrange 보간다항식을 이용한 수치적분

1. 학습 목표

이번 강의가 끝나면 다음을 할 수 있어야 한다.

  • 정적분을 함수 아래 넓이로 설명할 수 있다.
  • 사다리꼴 하나의 넓이를 계산할 수 있다.
  • 적분 구간을 여러 부분으로 나누어 사다리꼴 넓이를 더할 수 있다.
  • 간단한 수치적분을 Python으로 계산할 수 있다.
  • 주어진 점을 지나는 Lagrange 보간다항식을 만들고 적분할 수 있다.

2. 수치적분이란

정적분

\[I=\int_a^b f(x)\,dx\]

은 $x=a$부터 $x=b$까지 함수와 $x$축 사이의 부호 있는 넓이를 나타낸다.

간단한 함수는 적분 공식을 사용하여 정확한 값을 구할 수 있다. 하지만 함수가 복잡하거나 실험 결과가 표로만 주어지면 적분 공식을 사용하기 어렵다.

이때 함수 아래 영역을 단순한 도형으로 나누어 넓이를 더한다. 이를 수치적분(numerical integration)이라고 한다.

먼저 사다리꼴 공식을 익힌 뒤, 점들을 지나는 Lagrange 보간다항식을 적분하는 방법으로 확장한다.

3. 사다리꼴 공식

3.1. 한 구간에서 계산하기

두 점 $(a,f(a))$와 $(b,f(b))$를 직선으로 연결하면 함수 아래 영역을 사다리꼴로 근사할 수 있다.

사다리꼴의 평행한 두 변의 길이는 $f(a)$와 $f(b)$이고, 폭은 $b-a$이다. 따라서 넓이는

\[\boxed{ \int_a^b f(x)\,dx \approx \frac{b-a}{2}\left[f(a)+f(b)\right] }\]

이다.

즉, 다음 두 값을 곱하면 된다.

\[\text{사다리꼴 넓이} = \text{구간의 폭} \times \text{양 끝 함수값의 평균}\]

3.2. 여러 구간으로 나누기

넓은 구간을 사다리꼴 하나로 나타내면 곡선과 직선의 차이가 클 수 있다. 구간을 여러 개의 작은 구간으로 나누면 곡선의 모양을 더 잘 따라갈 수 있다.

구간 $[a,b]$를 $n$개로 똑같이 나누면 각 구간의 폭은

\[h=\frac{b-a}{n}\]

이다. 나눈 점을

\[x_0=a,\quad x_1=a+h,\quad \ldots,\quad x_n=b\]

라고 하자. 각 작은 사다리꼴의 넓이를 모두 더하면

\[\boxed{ I\approx \frac{h}{2} \sum_{i=0}^{n-1} \left[f(x_i)+f(x_{i+1})\right] }\]

이다.

복잡해 보이지만 실제 계산은 다음 세 단계뿐이다.

  1. 구간의 폭 $h$를 구한다.
  2. 서로 이웃한 두 점의 함수값으로 사다리꼴 넓이를 구한다.
  3. 모든 사다리꼴의 넓이를 더한다.

4. 손으로 계산하는 예제

다음 적분을 사다리꼴 공식으로 근사해 보자.

\[I=\int_0^2 x^2\,dx\]

구간 $[0,2]$를 두 부분으로 나누면

\[x_0=0,\qquad x_1=1,\qquad x_2=2\]

이고, 각 구간의 폭은 $h=1$이다. 함수값은

\[f(0)=0,\qquad f(1)=1,\qquad f(2)=4\]

이다.

첫 번째 사다리꼴의 넓이는

\[A_1 =\frac{1}{2}[f(0)+f(1)] =\frac{1}{2}(0+1) =0.5\]

이다. 두 번째 사다리꼴의 넓이는

\[A_2 =\frac{1}{2}[f(1)+f(2)] =\frac{1}{2}(1+4) =2.5\]

이다. 따라서 수치적분 결과는

\[I\approx A_1+A_2=3.0\]

이다.

정확한 값은

\[\int_0^2x^2\,dx=\frac{8}{3}\approx2.6667\]

이므로 오차는

\[|3.0-2.6667|\approx0.3333\]

이다.

5. Python으로 계산하기

먼저 손으로 계산한 과정을 그대로 Python으로 작성해 보자.

def f(x):
    return x**2

x = [0.0, 1.0, 2.0]
area = 0.0

for i in range(len(x) - 1): ## 0 ~ (n-1)
    width = x[i + 1] - x[i]
    average_height = (f(x[i]) + f(x[i + 1])) / 2.0
    area = area + width * average_height

print(area)

출력값은 $3.0$이다.

같은 계산을 여러 문제에 사용할 수 있도록 함수로 만들 수 있다.

import numpy as np

def trapezoidal(f, a, b, n):
    x = np.linspace(a, b, n + 1)
    area = 0.0

    for i in range(n):
        width = x[i + 1] - x[i]
        average_height = (f(x[i]) + f(x[i + 1])) / 2.0
        area = area + width * average_height

    return area

def func1(x):
    return x**2

def func2(x):
    return x**3

result1 = trapezoidal(func1, 0.0, 2.0, n=2)
print(result1)
result2 = trapezoidal(func2, 0.0, 2.0, n=2)
print(result1)

6. 구간을 더 잘게 나누면

같은 적분을 구간 수 $n$만 바꾸어 계산해 보자.

exact = 8.0 / 3.0

for n in [1, 2, 4, 8, 16]:
    result = trapezoidal(func1, 0.0, 2.0, n)
    error = abs(exact - result)

    print(f"n={n:2d}, result={result:.6f}, error={error:.6f}")

예상되는 결과는 다음과 같다.

n= 1, result=4.000000, error=1.333333
n= 2, result=3.000000, error=0.333333
n= 4, result=2.750000, error=0.083333
n= 8, result=2.687500, error=0.020833
n=16, result=2.671875, error=0.005208

구간을 더 잘게 나눌수록 사다리꼴이 곡선의 모양을 더 잘 따라가므로 결과가 정확한 값에 가까워진다.

그러나 구간 수를 무조건 크게 하는 것이 항상 좋은 것은 아니다. 계산 시간이 늘어나기 때문이다. 실제 문제에서는 구간 수를 조금씩 늘리면서 결과가 충분히 변하지 않는지 확인한다.

7. 재료공학 예제

재료의 비열 $C_p$는 온도에 따라 달라질 수 있다. 온도와 비열이 다음 표와 같이 측정되었다고 하자.

\[\begin{array}{c|cccc} T\ (\mathrm K) & 300 & 400 & 500 & 600\\ \hline C_p\ (\mathrm{J\,kg^{-1}K^{-1}}) & 500 & 520 & 550 & 590 \end{array}\]

재료 1 kg을 300 K에서 600 K까지 가열하는 데 필요한 열량은

\[Q=\int_{300}^{600}C_p(T)\,dT\]

로 계산한다. 각 온도 구간에서 사다리꼴 넓이를 구하면

\[\begin{aligned} Q_1 &=100\times\frac{500+520}{2}=51000,\\ Q_2 &=100\times\frac{520+550}{2}=53500,\\ Q_3 &=100\times\frac{550+590}{2}=57000 \end{aligned}\]

이다. 따라서

\[Q\approx Q_1+Q_2+Q_3 =161500\ \mathrm{J/kg} =161.5\ \mathrm{kJ/kg}\]

이다.

Python에서는 측정값을 배열로 저장하여 같은 계산을 할 수 있다.

temperature = np.array([300.0, 400.0, 500.0, 600.0])
heat_capacity = np.array([500.0, 520.0, 550.0, 590.0])

heat = 0.0

for i in range(len(temperature) - 1):
    width = temperature[i + 1] - temperature[i]
    average_cp = (heat_capacity[i] + heat_capacity[i + 1]) / 2.0
    heat = heat + width * average_cp

print(f"Q = {heat / 1000.0:.1f} kJ/kg")

8. Lagrange 보간다항식의 적분

사다리꼴 공식은 두 점 사이의 함수를 직선으로 근사한 뒤 적분한 결과이다. 세 점을 이용하면 곡선의 굽은 모양을 이차다항식(2nd order polynomial)으로 표현할 수 있다. 이처럼 주어진 점들을 정확히 지나는 다항식을 만드는 것을 보간 (interpolation)이라고 한다. 보간한 다항식을 적분하여 원래 함수의 적분을 근사한다.

8.1. 보간다항식 만들기

서로 다른 위치 $x_0,x_1,\ldots,x_{\textcolor{green}{n}}$에서 함수값 $y_i=f(x_i)$를 알고 있다고 하자. Lagrange 보간다항식은

\[p_{\textcolor{green}{n}}(x)=\sum_{i=0}^{\textcolor{green}{n}}y_iL_i(x),\qquad L_{\textcolor{red}{i}}(x)=\prod_{\substack{\textcolor{blue}{j}=0\\ \textcolor{blue}{j}\ne \textcolor{red}{i}}}^{\textcolor{green}{n}} \frac{x-x_{\textcolor{blue}{j}}}{x_{\textcolor{red}{i}}-x_{\textcolor{blue}{j}}}\]

이다. 곱 기호 $\prod_{\substack{i=0\j\ne i}}^n$는 $j=i$를 제외한 모든 항($i=0,…,n)$을 곱하라는 뜻이다. 각 기저다항식 $L_i$는 자기 점 $x_i$에서는 1이고 다른 점에서는 0이다. 따라서 $p_n(x_i)=y_i$가 되어 모든 데이터를 정확히 지난다.

예를 들어 세 점에서는

\[\begin{aligned} L_0(x)&=\cancel{\bigg(\frac{x-x_0}{x_0-x_0}\bigg)} \bigg(\frac{x-x_1}{x_0-x_1}\bigg)\bigg(\frac{x-x_2}{x_0-x_2}\bigg)=\frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)},\\ L_1(x)&=\bigg(\frac{x-x_0}{x_0-x_0}\bigg) \cancel{\bigg(\frac{x-x_1}{x_0-x_1}\bigg)}\bigg(\frac{x-x_2}{x_0-x_2}\bigg)=\frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)},\\ L_2(x)&=\bigg(\frac{x-x_0}{x_0-x_0}\bigg) \bigg(\frac{x-x_1}{x_0-x_1}\bigg)\cancel{\bigg(\frac{x-x_2}{x_0-x_2}\bigg)}=\frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)}. \end{aligned}\]

이를 적분하면 함수값의 가중합(weighted sum)으로 계산할 수 있다.

\[\boxed{\int_a^b f(x)\,dx\approx\int_a^b p_n(x)\,dx =\int_a^b \sum_{i=0}^n y_iL_i(x)\,dx =\sum_{i=0}^n\int_a^bL_i(x)y_i\,dx =\sum_{i=0}^{n}w_i y_i},\qquad w_i=\int_a^b L_i(x)\,dx\]

가중치 $w_i$는 점의 위치($x_0,x_1,…,x_{n}$)와 적분 구간($a,b$)으로 결정된다. 점의 간격이 같지 않아도 보간은 가능하지만, 위치가 바뀌면 가중치도 다시 계산해야 한다.

8.2. 두 점을 사용하면 사다리꼴 공식

양 끝점 $x_0=a$, $x_1=b$를 사용하면

\[p_1(x)=\sum_{i=0}^1y_iL_i(x)=y_0L_0(x)+y_1L_1(x)=y_0\frac{x-x_1}{x_0-x_1}+y_1\frac{x-x_0}{x_1-x_0}=y_0\frac{b-x}{b-a}+y_1\frac{x-a}{b-a}\]

이다. 두 기저함수를 각각 적분하면 모두 $(b-a)/2$가 된다, 즉 \(w_0=w_1=\frac{b-a}{2}\)

따라서

\[\int_a^b p_1(x)\,dx=\sum_{i=0}^nw_iy_i=w_0y_0+w_1y_1=\frac{b-a}{2}(y_0+y_1)\]

이다. 앞에서 배운 사다리꼴 공식이 바로 일차 보간다항식의 적분이다.

8.3. 세 점을 사용한 손계산

4절과 동일하게 $(x_0,y_0)=(0,0)$, $(x_1,y_1)=(1,1)$, $(x_2,y_2)=(2,4)$를 사용하자. 기저다항식은

\[L_0=\frac{(x-1)(x-2)}{2},\qquad L_1=2x-x^2,\qquad L_2=\frac{x(x-1)}{2}\]

이다. 각 기저를 $0$부터 $2$까지 적분하면

\[\begin{aligned} w_0&=\left[\frac{x^3}{6}-\frac{3x^2}{4}+x\right]_0^2=\frac13,\\ w_1&=\left[x^2-\frac{x^3}{3}\right]_0^2=\frac43,\\ w_2&=\left[\frac{x^3}{6}-\frac{x^2}{4}\right]_0^2=\frac13. \end{aligned}\]

따라서 적분값은

\[I\approx\frac13(0)+\frac43(1)+\frac13(4)=\frac83\]

이다. 이 예제에서는 $p_2(x)=0L_0+L_1+4L_2=x^2$이므로 정확한 적분값과 같다. 같은 세 점으로 사다리꼴 두 개를 만들면 3이었지만, 이차 보간은 곡률을 반영한다. 일반 함수에서는 보간다항식과 원래 함수 사이의 차이 때문에 적분 오차가 남는다.

8.4. Simpson 공식과 Python 계산

등간격 세 점 $x_0=a$, $x_1=a+h$, $x_2=a+2h=b$로 위 계산(즉 이차 다항식으로 보간)을 반복하면 가중치가 $h/3$, $4h/3$, $h/3$이 된다. 이것이 Simpson의 1/3 공식이다.

\[\boxed{\int_a^b f(x)\,dx\approx \frac{h}{3}(y_0+4y_1+y_2)},\qquad h=\frac{b-a}{2}\]

여기서 $h$는 전체 폭이 아니라 이웃한 두 점 사이의 거리이다. 다음 코드는 등간격 세 점의 이차 Lagrange 보간다항식을 적분한 공식을 구현한다.

def integrate_three_points(a, b, y0, y1, y2):
    # y0, y1, y2는 a, (a+b)/2, b에서의 함수값이다.
    if b <= a:
        raise ValueError("이 예제에서는 a < b인 구간을 사용합니다.")
    h = (b - a) / 2.0
    return h * (y0 + 4.0*y1 + y2) / 3.0


result = integrate_three_points(0.0, 2.0, 0.0, 1.0, 4.0)
print(f"적분값 = {result:.6f}")  # 2.666667

재료공학 예제로 7절의 첫 세 측정점만 사용해 300 K에서 500 K까지의 단위 질량당 열량을 계산해 보자. 다시 테이블을 옮기자면 아래와 같다.

\[\begin{array}{c|cccc} T\ (\mathrm K) & 300 & 400 & 500 & 600\\ \hline C_p\ (\mathrm{J\,kg^{-1}K^{-1}}) & 500 & 520 & 550 & 590 \end{array}\] \[\int_{300}^{500}C_p(T)\,dT \approx\frac{100}{3}(500+4\times520+550) =104333.3\ \mathrm{J/kg} \approx 104.3\ \mathrm{kJ/kg}\] 
heat_per_mass = integrate_three_points(
    300.0, 500.0, 500.0, 520.0, 550.0
)
print(f"단위 질량당 열량 = {heat_per_mass / 1000.0:.3f} kJ/kg")

같은 구간의 사다리꼴 결과는 $104.5\ \mathrm{kJ/kg}$이다. 측정점 사이의 실제 비열 함수를 모르므로 두 결과의 차이만으로 정확한 오차를 알 수는 없으나, 사다리꼴 결과와 보간다항식 적분 결과가 매우 유사한 것을 알 수 있다.

점의 수를 계속 늘려 높은 차수의 다항식 하나로 보간하면 진동이 커질 수 있다. 실제로는 작은 구간마다 낮은 차수의 보간다항식을 적분하는 방법을 많이 사용한다. 위 Simpson 공식을 여러 구간에 반복 적용하려면 등간격의 작은 구간을 두 개씩 묶어야 하므로 전체 작은 구간 수는 짝수여야 한다.

9. 정리

  • 수치적분은 복잡한 함수 아래의 넓이를 단순한 도형의 넓이로 근사하는 방법이다.
  • 사다리꼴 넓이는 구간의 폭과 양 끝 함수값의 평균을 곱하여 구한다.
  • 여러 개의 작은 사다리꼴을 사용하면 일반적으로 더 정확한 결과를 얻는다.
  • 실험 데이터가 표로만 주어진 경우에도 수치적분을 사용할 수 있다.
  • 구간 수를 늘리면서 계산 결과가 일정한 값에 가까워지는지 확인해야 한다.
  • Lagrange 보간다항식의 적분은 주어진 함수값에 가중치를 곱해 더하는 계산이다.
  • 두 점의 일차 보간에서 사다리꼴 공식, 등간격 세 점의 이차 보간에서 Simpson 공식을 얻는다.

10. 연습 문제

문제 1. 사다리꼴 넓이

폭이 2이고 양 끝의 높이가 각각 3과 5인 사다리꼴의 넓이를 구하라.

문제 2. 한 구간의 적분

사다리꼴 공식을 한 번 사용하여 다음 적분을 근사하라.

\[\int_0^2 x\,dx\]

문제 3. 구간의 폭

구간 $[0,4]$를 같은 크기의 4개 구간으로 나눌 때 $h$와 나눈 점들을 구하라.

문제 4. 두 사다리꼴의 합

다음 표를 이용하여 $x=0$부터 $x=2$까지의 넓이를 구하라.

\[\begin{array}{c|ccc} x & 0 & 1 & 2\\ \hline f(x) & 1 & 2 & 3 \end{array}\]

문제 5. 구간 수와 정확도

일반적으로 적분 구간을 더 잘게 나누면 계산 결과가 더 정확해지는 이유를 한 문장으로 설명하라.

문제 6. Python 실습

위의 trapezoidal 함수를 사용하여

\[\int_0^1 x^2\,dx\]

를 $n=1$, $n=2$, $n=4$로 계산하라. 정확한 값 $1/3$과 비교하라.

문제 7. 열량 계산

비열이 300 K에서 $400\ \mathrm{J\,kg^{-1}K^{-1}}$이고, 400 K에서 $440\ \mathrm{J\,kg^{-1}K^{-1}}$이다. 사다리꼴 공식으로 1 kg의 재료를 300 K에서 400 K까지 가열하는 데 필요한 열량을 구하라.

문제 8. 기저다항식 확인

8절의 $L_0(x)=(x-1)(x-2)/2$에 $x=0,1,2$를 각각 대입하라.

문제 9. 세 점의 가중합

$x=0,1,2$에서 함수값이 각각 $1,2,5$이다. 가중치 $1/3,4/3,1/3$을 사용하여 $0$부터 $2$까지의 적분을 근사하라.

문제 10. 간격과 적분값

$x=0,2,4$에서 함수값이 각각 $0,4,16$이다. 이웃한 점의 간격 $h$를 구하고 Simpson 공식으로 적분하라. 8절의 Python 함수로도 확인하라.

문제 11. 공식의 적용 조건

$x=0,1,3$에서 얻은 함수값에 $h(y_0+4y_1+y_2)/3$을 바로 적용해도 되는가?