Python에서 미분(차분)을 한다는 것. 데이터, 혹은 함수

어떤 형태든 센서 신호를 만지작 거리고, 모터를 구동하고, 뭐 그러다보면 미분(혹은 차분)을 수행해야할 경우가 생깁니다. 오늘은 파이썬으로 미분(차분)하는 일에 대해 이야기를 해보려고 합니다.^^

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

t = np.arange(0, 2*np.pi, 0.1)
y1 = np.sin(t)

plt.figure(figsize=(12,6))
plt.plot(t, y1);

일단, 시험 신호를 하나 만들어 두겠습니다. 시간축 t는 0부터 2pi까지 0.1 간격을 가지도록 했습니다.

이렇게 생긴거죠^^ 먼저, 위 시험신호처럼 우리가 함수를 아는 경우에 대해 접근해 보겠습니다.

바로 고등학교때 배운 도함수의 정의식으로 접근합니다. 내가 원함수를 알고 있다면, 그러면 위 수식에서 함수 f()를 알수 있으니, 위의 정의식을 적용할 수 있습니다. 그걸 그냥 코드로 작성하면, 

def derivative(f, a, h=0.01):
    return (f(a + h) - f(a))/h

이러합니다. 여기서 h값만 지정되지 않을 경우 0.01로 두라고만 해두었습니다. (이 함수와 그래프 그리는 방식은 math.ubc.ca라는 사이트의 좋은 글에서 가져왔습니다.)

dy1dx = derivative(np.sin, t)

plt.figure(figsize=(12,5))
plt.plot(t, dy1dx, 'r.', label='Calculated Difference')
plt.plot(t, np.cos(t), 'b', label='True Value')
plt.plot(t, y1, 'k', label='Original Value')
plt.legend(loc='best')
plt.show()

이제 derivative 함수를 이용해서 원함수 y1과 미분한 함수와 sin을 미분하면 cos인걸 아니까 그걸 비교해 보도록 하겠습니다.

넵... 도함수의 정의식으로 미분한 것이 실제 미분한 것과 큰 차이가 없네요. (사실 차이가 있지만, 미미하다고 해 두겠습니다) 그런데 만약 원함수를 알 수 없다면 어떻게 할까요. 사실 이런 경우가 더 많겠죠.

y1_d = np.r_[0, np.diff(y1)]/0.1

이럴때는 numpy가 제공하는 diff함수를 사용해서 차분을 구현합니다. numpy의 diff 함수는 현재값에서 이전값을 뺀것을 반환합니다. 거기에 시간간격(0.1)으로 나눠주면 되는 거죠. (차분의 정의식대로 입니다.) 여기서 사용된 np.r_ 함수는 numpy 배열에서 행이나 열을 추가하는 함수입니다. 이 함수를 사용한 이유는 diff를 사용하면 크기가 하나 줄기 때문입니다. 그래서 처음에 0을 하나 추가해서 크기를 맞췄습니다.

plt.figure(figsize=(12,5))
plt.plot(t, y1_d, 'r.', label='Calculated Difference')
plt.plot(t, np.cos(t), 'b', label='True Value')
plt.plot(t, y1, 'k', label='Original Value')
plt.legend(loc='best')
plt.show()

그리고 그려보죠

이렇습니다. 약간 차이가 나네요. 차분 오차입니다.

살짝 이야기를 뒤로 돌려서 numpy가 제공하는 함수가 아니라, 일반함수라면 도함수의 정의식을 구현했던 derivative 함수를 어떻게 사용할까요. lambda 함수를 이용하면 됩니다.

x = np.linspace(0,6,100)
f = lambda x: x**2
y = f(x)
dydx = derivative(f,x)

plt.figure(figsize=(12,5))
plt.plot(x, y, label='y=f(x)')
plt.plot(x, dydx, label="Central Difference y=f'(x)")
plt.legend()
plt.grid(True)
plt.show()

이렇게 인데요.

이렇게 사용할 수 있습니다. 이제, 도함수든 차분방정식이든 하나의 함수(def)에 넣어서 함께 사용하고 싶습니다. 그러면 입력이 데이터인지 함수인지 알아야합니다.

데이터 타입으로 결정할 수도 있습니다만, callable 한지 아닌지로 판단해도 됩니다.

def derivative(f, a, h=0.01):
    if callable(f):
        return (f(a + h) - f(a))/h
    else:
        return np.r_[0, np.diff(f)]/h

이렇게 함수를 만들어 두면 됩니다. callable하면 도함수의 정의식을 아니면 차분방정식을 쓰겠다는 거죠^^

ts = 0.01
t = np.arange(0, 2*np.pi, ts)
y = np.sin(t)
y_d = derivative(y, t, h=ts)
y_df = derivative(np.sin, t, h=ts)

plt.figure(figsize=(12,5))
plt.plot(t, y_d, 'r.', label='Differential Method')
plt.plot(t, y_df, 'b', label='Numerical Method')
plt.plot(t, y, 'k', label='Original Value')
plt.legend()
plt.grid(True)
plt.show()

이제 테스트~

이렇게 됩니다. 결과가 괜찮은듯 합니다.^^

확대해보면 원값과 둘 다 오차를 가지지만, 이는 샘플링타임을 조절해도 되고, 이 정도 오차는 제가 사용하는 범위에서는 큰 문제가 없어서 그냥 쓰면 될 것 같습니다.^^