# 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:x[0]+x[1]+x[1], f[:-1].reshape(len(f)/2,2)))-f[0]+f[-1])*h/6

# 홀수개만큼 있을 때는 위의 코드를 쓰면 된다. 왜 되는지는 다시 잘 생각해 보도록 하자. 파이썬 공부 겸 아이큐테스트라고 하자. (나도 테스트 안해봄.)
# 짝수개일 때는 어떻게 될까? 짜투리 부분의 오차를 어떻게 정리할 수 있을지는 각자 생각해 보자.

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

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