파이썬 적분 라이브러리 - paisseon jeogbun laibeuleoli

파이썬 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 라...