파이썬 Scipy 수학 라이브러리
파이썬 Scipy 수학 라이브러리
파이썬 라이브러리¶
선형대수 연산¶
scipy 는 거의 항상 numpy 와 같이 사용된다
선형 대수 연산을 지원하는 scipy.linalg 모듈을 불러온다
In [1]:
import numpy as np from scipy.linalg import * # 선형대수 라이브러리
행렬은 np.array() 함수를 이용하여 만든다.
In [2]:
arr = np.array([[1, 2],[3, 4]]) arr
Out[4]:
array([[-2. , 1. ], [ 1.5, -0.5]])
Out[5]:
array([[1.0000000e+00, 0.0000000e+00], [8.8817842e-16, 1.0000000e+00]])
행렬로 연립 방정식 풀기
$$\begin{eqnarray} 3x+2y & = & 2 \\ x-y & = & 4 \\ 5y+z & = & -1 \end{eqnarray}$$
In [6]:
A = np.array([[3, 2, 0], [1, -1, 0], [0, 5, 1]]) b = np.array([2, 4, -1]) # 열벡터 x = solve(A, b) x
비선형 방정식의 해 구하기¶
맷플롯립 패키지를 불러오고, 그래프를 노트북 문서 안에 나타낼 수 있도록 설정한다
In [8]:
import matplotlib.pyplot as plt %matplotlib inline
다음 방정식의 해를 구한다. $$x=e^{-x}$$
좌우변을 그래프로 그려보면
In [9]:
x = np.linspace(-1,1,101) y = np.exp( -x ) plt.plot( x, y ) plt.plot( x, x ) plt.show()
방정식을 풀기 위하여 scipy.optimize 모듈을 불러온다
In [10]:
from math import * # 일반적인 수학 함수를 사용하려면 필요하다. from scipy.optimize import fsolve # fsolve( ) 함수를 도입한다.
풀고자 하는 방정식을 $ f(x)=0 $ 형태로 두고, 함수로 정의한다
In [11]:
def f(x): return x - exp(-x)
해를 구한다: fsolve( 함수, 초기값 )
함수의 최소값¶
scipy.optimize 모듈의 fmin_bfgs() 함수로 최소값을 구한다.
In [13]:
from scipy.optimize import fmin_bfgs
$f(x)= x^2 + 10 \sin (x)$ 의 최소값을 구한다.
파이썬 함수로 정의하고, 함수의 형태를 그래프로 그려본다.
In [14]:
def f(x): return x**2 + 10*np.sin(x) x = np.arange(-10, 10, 0.1) plt.plot(x, f(x)) plt.show()
fmin_bfgs() 함수의 인자로서 최소화할 함수 이름과 x의 시작점을 넘겨준다.
Optimization terminated successfully. Current function value: -7.945823 Iterations: 5 Function evaluations: 18 Gradient evaluations: 6
fmin_bfgs() 함수는 국부적인 최소점을 찾는 기능만 있으므로, 초기값에 따라 결과가 달라질 수 있다.
만약, x = 3 을 시작점으로 하면
Optimization terminated successfully. Current function value: 8.315586 Iterations: 6 Function evaluations: 21 Gradient evaluations: 7
수치 적분¶
scipy.integrate 모듈의 quad() 함수를 이용하여 수치 적분을 계산한다.
In [17]:
from scipy.integrate import quad fx = lambda x: x**2 quad( fx, 0, 1 )
Out[17]:
(0.33333333333333337, 3.700743415417189e-15)
첫번째 반환값이 수치 적분으로 얻어진 결과이며, 두번째의 값은 추정 오차이다.
피적분함수를 함수로서 정의하여 quad() 함수의 인자로 넘겨줄 수 있으며, 수학식을 표현할때 numpy 함수를 사용한다.
$$ \int _0 ^{\infty} e^{-x} dx $$
In [18]:
def f(x) : return np.exp(-x) F, err = quad( f, 0, np.inf ) print( F )
선형 회귀법 (linear regression): 최소제곱법¶
scipy.stats 모듈의 linregress( x, y ) 함수를 사용하여 회기 직선의 기울기와 절편을 구한다.
In [19]:
import numpy as np from scipy.stats import linregress x = np.linspace(0, 4, 5) y = [ 0, 2.1, 4.2, 5.9, 8.3 ] slope, intercept, r_value, p_value, std_err = linregress(x, y) print("slope = %8.3f" % slope) print("intercept = %8.3f" % intercept)
slope = 2.040 intercept = 0.020
In [20]:
yfit = intercept + slope * x plt.plot(x, y, 'o', label='data') plt.plot(x, yfit, 'r', label='fitted line') plt.legend() plt.xlabel('x') plt.ylabel('y') plt.show()
비선형 회귀법 (nonlinear regression): Levenberg-Marquardt 알고리즘¶
다음 그림에 주어진 데이터를 잘 표현하는 수식을 비선형 회귀법으로 얻으려고 한다.
In [21]:
ndata = 101 x = np.linspace(-5,5,ndata) y = 0.6*x + 5* np.exp(-0.2* x**2 ) + np.random.rand(ndata) plt.plot(x,y,'o')
Out[21]:
[<matplotlib.lines.Line2D at 0xb00add0>]
Levenberg-Marquardt 알고리즘을 이용하여 최적의 곡선을 구한다.
scipy.optimize 모듈의 curve_fit() 함수를 이용한다.
원하는 형태의 곡선에 대한 수식을 함수로 정의한다. 함수의 수식을 표현할때, numpy 수학 함수를 사용한다.
$$ y = a x + b { e }^{ -c{ x }^{ 2 } } $$
계수 $a, b, c$ 는 최적화에 의해 결정되어질 값이다.
In [22]:
# non-linear least square method using Levenberg-Marquardt algorithm from scipy.optimize import curve_fit def func(x, a, b, c): return a*x + b*np.exp(-c* x**2 )
매개변수의 초기값을 설정하고, curve_fit 함수를 호출한다. 첫번째로 반환되는 변수 popt 가 최적화된 매개변수에 대한 리스트이다.
In [23]:
p0ini = [0, 0, 0] # initial values of parameters for fitting function popt, pcov = curve_fit( func, x, y, p0 = p0ini ) print( popt )
[0.62760057 5.27337533 0.14737132]
구해진 매개변수로 최적화된 곡선을 예측하고, 원래의 주어진 데이터와 비교한다.
In [24]:
yfit = func( x, *popt ) # 매개변수를 함수의 인자 a, b, c 로 전달하여 계산한다. plt.plot( x, y, 'o', label='data') plt.plot( x, yfit, 'r', label='fitted curve') plt.legend() plt.xlabel('x') plt.ylabel('y') plt.show()
상미분 방정식 - 초기값 문제¶
scipy.integrate 모듈의 solve_ivp() 함수를 이용하여 상미분방정식의 초기값 문제에 대한 수치해를 구할 수 있다.
In [25]:
from scipy.integrate import solve_ivp
일계 미분 방정식의 수치해를 구하는 경우를 고려해 보자.
$$ \frac {dy} {dt} = - 0.5 y $$
도함수를 함수 형식으로 정의한다.
In [26]:
def dydt(t, y) : return -0.5 * y
초기값 $y(0)=1$ 에 대하여, $t=10$ 까지 수치해를 구한다.
In [27]:
t0 = 0 tend = 10 y0 = 1 sol = solve_ivp( dydt, (t0, tend), (y0,), t_eval=np.linspace(t0,tend,21) )
초기값은 한 개일 경우에도, 튜플의 형식으로 전달한다.
인자 t_eval 은 선택사항으로서, 수치해가 계산될 $t$ 의 간격을 명시할 수 있다.
결과가 반환된 변수 sol 에서 $t$ 와 방정식의 해 $y$ 는 각각 다음과 같이 얻어진다.
Out[28]:
array([ 0. , 0.5, 1. , 1.5, 2. , 2.5, 3. , 3.5, 4. , 4.5, 5. , 5.5, 6. , 6.5, 7. , 7.5, 8. , 8.5, 9. , 9.5, 10. ])
Out[29]:
array([[1. , 0.77877121, 0.60652683, 0.47231798, 0.36766997, 0.28642722, 0.22325718, 0.17381699, 0.13533366, 0.10546992, 0.08217585, 0.06396748, 0.04982315, 0.03883629, 0.03024224, 0.0235437 , 0.01834532, 0.01429687, 0.01113201, 0.00867035, 0.0067547 ]])
결과를 그래프로 그린다. sol.y 는 2차원 배열이므로 첫번째 행을 취한다.
In [30]:
plt.plot( sol.t, sol.y[0] ) plt.xlabel('t') plt.ylabel('y') plt.show()
이계 미분방정식¶
$$ y'' +2y'+2y = \cos (2t), \qquad y(0)=0, \quad y'(0)=0$$
이계 미분방정식은 연립 일계 미분방정식으로 변환할 수 있다. $v (t) = y'(t)$ 로 두면
$$ v ' + 2 v + 2y = \cos (2t), \qquad y(0)=0, \quad v(0)=0$$
In [31]:
def dYdt( t, Y ) : # Y[0] = y(t), Y[1] = v(t) return Y[1], -2*Y[1]-2*Y[0]+np.cos(2*t) # y'(t), v'(t) # 초기값 Y0 = ( 0, 0 ) # y(0) v(0) t0 = 0 tend = 10 sol = solve_ivp( dYdt, (t0, tend), Y0, t_eval=np.linspace(t0,tend,101) )
In [32]:
plt.plot( sol.t, sol.y[0] ) plt.xlabel('t') plt.ylabel('y') plt.show()
Numeric Analysis 4 - Numeric Linear Algebra
Numeric Analysis 4 - Numeric Linear Algebra Numeric Linear Algebra ¶ ...
파이썬 Sympy 기호수학 - 기초 $\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad...
파이썬 Scipy 수학 라이브러리 파이썬 라이브러리 ¶ 참고 자료 : Scipy Lecture Notes ...
파이썬 Sympy 기호수학 - 1. 상미분방정식 SymPy 사용법 ¶ SymPy Tutorial sympy 라...