lecture_1.py


# Elementary Numerical analysis 1
# based on Python
# (C) 2013. Keehwan Nam, Dept. of physics, KAIST.
# snowall@gmail.com / snowall@kaist.ac.kr

# Differentiation
# 테일러 정리와 테일러 전개를 적극 사용하면 된다.
# 수치적으로 미분은 다음과 같이 정의한다
# df/dx = (f(x+h) - f(x))/h
# 이 경우 오차가 h의 제곱에 비례한다. 2점 미분공식이다.

# 오차를 줄이고 싶은 경우 대칭성을 이용해서 오차를 줄일 수 있다.
# df/dx = (f(x+h)-f(x))/h-(f(x-h)-f(x))/h = (f(x+h)-f(x-h))/2h
# 이 경우 오차가 h의 세제곱에 비례한다. 3점 미분공식이다.

# 더 줄이고 싶다면 위의 3점 미분공식을 5점으로 확장할 수 있다.
# df/dx = (f(x+2h)+f(x+h)-f(x-h)-f(x-2h))/6h
# 하지만 이래봐야 오차는 그냥 h의 세제곱에 비례하게 된다.

# 머리를 써 보면 5점 미분공식은 다음과 같다.
# df/dx = (-f(x+2h)+8f(x+h)-8f(x-h)+f(x-2h))/12h
# 이렇게 하면 오차가 h의 4제곱에 비례하게 된다.

# 2차 이상의 미분을 계산해야 한다면, 1차 미분을 반복적으로 계산할 수도 있고, 또는 직접 2차 도함수를 얻을 수도 있다.
# d2f/dx2 = 2(f(x+h)-2f(x)+f(x-h))/(h*h)

# 이와 같이 n차 도함수를 직접 계산하는 공식은 테일러 전개에 의하여 손쉽게 얻을 수 있다.

# 실제 구현은?

# 파이썬이라면, 수치가 리스트로 주어져 있을 것이다.
f = numpy.array([x1, x2, x3, x4, x5])

h = 1.0e-2 # 미분에서 사용한 h값이 된다.

# 첫번째로, 2점 미분공식을 얻어보자.

df = (f[1:]-f[:-1])/h
# 이 공식이 왜 작동하는지는 차근차근 생각해 보면 된다.

# 마찬가지로, 3점 미분공식도 똑같이 얻을 수 있다.
df = (f[2:]-f[:-2])/(2*h)
   
# 5점 미분공식은 직접 고민해 보자.
# 하지만, 항상 리스트로만 주어져 있을리가 없다. 화일에서 직접 값을 한줄씩 읽어오면서 계산해야 하는 경우도 있을 것이다.
# 물론 화일을 리스트에 집어넣고나서 위의 방법을 쓰면 되지만,
# 화일이 메모리에 안 들어가는, 가령 수십 GB용량의 데이터를 미분해야 한다면 그건 그거대로 막막하리라고 본다.

data = open("mydata.txt", "r")
f = [data.readline(), data.readline(), data.readline()]
df = []
while data.EoF():
    df.append((f[0]-f[2])/(2.*h))
    f.append(data.readline())
    f=f[1:]
   
# data.EoF()는 파일의 끝인지 아닌지 판단하는 함수이다. (구현은 셀프!)
# 이게 왜 작동하는지는 직접 생각해 보고, 5점 미분공식으로 바꾸려면 어떻게 해야 할지 생각해 보자. 또는 2차도함수 계산이라면?
# df.append로 리스트에 집어넣는 것이 메모리를 잡아먹을 것 같다면, df를 파일로 열어서 df.writeline함수를 사용해도 된다.
# 다른 프로그래밍 언어라면 구체적으로는 다르게 구현해야 할 것이다. 자신이 사용하는 언어에서는 어떻게 구현할 수 있을까?

# 파이썬의 scpiy 모듈에서 미분을 제공하는데, scipy.misc.derivative(func, x0, dx=1.0, n=1, args=(), order=3) 함수를 사용하면 된다.
# 이 함수는 func로 주어진 함수가 x0위치에서 주어진 dx로 미분한 n차 "미분계수" 값을 알려준다. order는 3점공식인가 5점공식인가 등등이다. '당연히' 홀수만 넣어야 한다.
# 내가 정의한 함수가 수를 넣으면 수가 나오는 함수라면, 위의 함수로 주어진 점에서의 미분계수를 알아낼 수 있다. 그러나 해석적인 도함수를 알려주지는 않는다.


# Integration

# 용어를 먼저 알고 가자.
# Quadrature 또는 Numerical quadrature는 수치해석에서 말하는 수치 적분을 뜻한다. 1차원 이상의 고차원 수치적분을 포함하기는 하지만, 1차원 수치적분의 뜻이 강하다.
# Cubature는 1차원 이상의 고차원 수치적분을 대체로 뜻한다.

# 1차원 적분을 먼저 논의해 보자.

# 가장 간단하게는 사각형으로 근사할 수 있다. 원래 적분은 연속 함수를 잘게 썰어서 각각의 넓이를 구한 후, 더 많은 조각으로 더 잘게 썰어가며 극한을 구하는 과정이다.
# 따라서 컴퓨터에게 수치적으로 적분을 계산을 시킨다면, 충분히 많은 조각으로 잘게 썰어서 각각의 넓이를 구한 후 다 더하도록 하면 된다.
# 가장 간단하게 사각형으로 근사한다면, i번재 함수값을 f[i]라 한다면
# F = sum(f[i], i=0 to N)*(b-a)/N
# 물론 각 조각은 N등분했다고 가정했다. 조각의 길이가 구간마다 다른 경우라면 골치아파지므로 그러지 말자.
# 위의 코드는 파이썬으로 간단히 구현할 수 있다.

f = [f1, f2, f3, f4, ,f5] # f는 함수값들을 수치적으로 저장하고 있는 리스트이다.
F=0
for i in f:
    F+=i
F/=h # h는 구간의 길이이다. h=(b-a)/N 정도로 정의하면 될 듯.
# 대충 이렇게 쓰면 될 것이다. 만약 조금 더 짧게 쓰고 싶다면 reduce와 lambda를 쓸 수 있다. 이 두 구문이 뭔지에 대해서는 자습해보자.

F = reduce(lambda x, y:x+y, f)*h

# 위의 구문이 적분식이다. 이걸 midpoint rule 이라고 한다.

# 사각형으로 적분하는 것보다, 아무래도 사다리꼴로 나누는 것이 조금 더 정확하지 않을까? 하는 마음에 만든 공식이 사다리꼴 규칙이다. Trapezoidal rule (trapezium rule)

F = reduce(lambda x, y:x+y, f)*h - 0.5*(f[0]+f[-1])*h

# 위의 구문이 사다리꼴 적분에 해당한다. 별 차이 없어 보이는데 아주 조금 더 정확하다. 왜그런지는 어렵지 않으니 직접 생각해 보자.
# 이 경우 오차는 대략 구간의 크기인 h의 제곱에 비례한다.

# 보다 정확한 근사식은 Simpson's rule이 제공한다. 사다리꼴 규칙에서 심슨 규칙으로 진화하는 것은 미분에서 2점 공식을 3점 공식으로 바꾸는 과정과 유사하다.
# 즉, 사다리꼴 규칙은 2점을 이용해서 그 사이의 면적을 근사했지만, 심슨 규칙은 3점을 이용해서 그 사이의 면적을 근사한다.
# 즉, 3점으로 이루어진 2차곡선이 원래 주어진 곡선보다 위인 구간과 아래인 구간이 있게 되는데, 넘치는 부분과 모자라는 부분이 상쇄되어 보다 정확한 근사값을 얻는다. 곡선의 근사는 틀리게 되지만, 적분값은 정확해진다는 것이다.
# 심슨 규칙은 다음과 같다.
# F=(f(a)+4f((a+b)/2)+f(b))*(b-a)/6
# 딱 3점, 즉 a, b, (a+b)/2의 3군데 값을 알고 있을 때 그 사이의 면적은 위와 같이 근사된다. 이 경우 오차는 구간의 크기인 b-a의 5제곱에 비례한다.
# 파이썬으로 구현한다면 어떻게 될까?

# 가장 간단하게 구현한다면 다음과 같을 것이다.
f = [f1, f2, f3]
F = (f[0]+4*f[1]+f[2])*h/6

# 하지만 이런 경우는 점이 정말 3개밖에 없는 경우이다. 만약 그보다 많다면? 일단 정확히 홀수개 있다고 하자.
f = [f1, f2, f3, f4, f5] # 리스트는 더 길어질 수도 있다.
F = (2.*reduce(lambda x, y: x+y, map(lambda x:x[0]+x[1]+x[1],f[:-1].reshape(len(f)/2,2)))-f[0]+f[-1])*h/6
# 이 코드는 간단해 보여도 자그마치 3개의 파이썬 기능이 합쳐진 콤비네이션이다.
# 홀수개만큼 있을 때는 위의 코드를 쓰면 된다. 왜 되는지는 다시 잘 생각해 보도록 하자. 파이썬 공부 겸 아이큐테스트라고 하자. (나도 테스트 안해봄.)
# 짝수개일 때는 어떻게 될까? 짜투리 부분의 오차를 어떻게 정리할 수 있을지는 각자 생각해 보자.

# 파이썬에서 구현된 적분은 scipy에서 제공한다.
# scipy.integrate(f, a, b)는 주어진 함수 f를 a에서 b까지 적분해준다. 다중적분도 다룰 수 있다.

# 1차원 적분의 경우 위의 코드로만 적분하더라도 해볼만하다. 문제는 다차원 적분은 적분 구간을 (a, b)처럼 쉽게 쪼개줄 수 없다는 것이다.
# 물론 잘 쪼개서 다중 반복문을 돌리면 된다. 당연히 된다. 만약 함수가 수치적으로, 즉 모든 위치에서의 함수값이 '진짜 값'으로 주어져 있다면 다중반복문을 돌려서 적분해야 할 것이다.
# 하지만, 만약 고차원에서 정의된 매우 이상한 함수가 있는데, 이 함수의 함수값은 구할 수 있지만 함수값이 다 주어진게 아니라 그냥 함수만 주어져 있다면? 그리고 그 해석적인 적분은 구할 수 없는 경우라면?
# 이때는 확률적 방법을 통해서 구할 수 있는데 그것이 바로 Monte-carlo integration이다.
# 다차원 다중적분 및 MC는 다른 기회에 설명하도록 하겠다. 


---

Quadratic은 "사각형의"라는 뜻을 담고 있는데, 그래서 "제곱"이라는 뜻도 있다. 일반적으로 사각형은 2차원 평면에서 정의된 도형이고, 1차원에서 1차원으로 가는 함수는 보통 직사각형으로 근사해서 적분하게 되므로 quadrature는 2차원 적분의 뜻이 강하다. Cube는 정육면체를 뜻하고, Cubic은 그래서 세제곱이라는 뜻이 있다. 거기서 나온 Cubature는 3차원 적분의 뜻이 강하다. 3차원 이상에도 각 차원마다 라틴어 어원을 갖는 형용사들이 있겠지만 그냥 Cubature라고 부르는 듯. Curvature는 "곡률"이라는 뜻으로 cubature와는 아무 관련 없다.

by snowall 2013. 9. 2. 21:54