파이썬 라이브러리¶
scipy 라이브러리¶
선형대수 연산¶
scipy 는 거의 항상 numpy 와 같이 사용된다
선형 대수 연산을 지원하는 scipy.linalg 모듈을 불러온다 (또는 numpy.linalg 모듈도 가능하다)
import numpy as np
from scipy.linalg import * # 선형대수 라이브러리
행렬은 np.array() 함수를 이용하여 만든다.
arr = np.array([[1, 2],[3, 4]])
arr
array([[1, 2],
[3, 4]])
행렬식의 계산: det() 함수
det(arr)
np.float64(-2.0)
역행렬 구하기: inv() 함수
iarr = inv(arr)
iarr
array([[-2. , 1. ],
[ 1.5, -0.5]])
행렬의 곱셈 : np.dot( A, B )
np.dot( arr, iarr )
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}$$
A = np.array([[3, 2, 0], [1, -1, 0], [0, 5, 1]])
b = np.array([2, 4, -1]) # 열벡터
x = solve(A, b)
x
array([ 2., -2., 9.])
정답 확인
np.dot( A, x )
array([ 2., 4., -1.])
비선형 방정식의 해 구하기¶
맷플롯립 패키지를 불러오고, 그래프를 노트북 문서 안에 나타낼 수 있도록 설정한다
import matplotlib.pyplot as plt
다음 방정식의 해를 구한다.
$$x=e^{-x}$$
좌우변을 그래프로 그려보면
x = np.linspace(-1,1,101)
y = np.exp( -x )
plt.plot( x, y )
plt.plot( x, x )
plt.show()
0.5 근처에서 해가 존재함을 알 수 있다.
방정식을 풀기 위하여 scipy.optimize 모듈을 불러온다
from math import * # 일반적인 수학 함수를 사용하려면 필요하다.
from scipy.optimize import fsolve # fsolve( ) 함수를 도입한다.
풀고자 하는 방정식을 $ f(x)=0 $ 형태로 두고, 함수로 정의한다
def f(x):
return x - exp(-x)
해를 구한다: fsolve( 함수, 초기값 )
fsolve( f, 0.5 )
C:\Users\Public\Documents\ESTsoft\CreatorTemp\ipykernel_13544\2602480254.py:2: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.) return x - exp(-x)
array([0.56714329])
경고를 피하려면 배열 인수를 지원하는 numpy 수학 함수를 쓴다.
def f(x):
return x - np.exp(-x)
fsolve( f, 0.5 )
array([0.56714329])
lambda 함수를 쓰면, 더욱 간략하게 된다.
fsolve( lambda x: x-np.exp(-x), 0.5 )
array([0.56714329])
함수의 최소값¶
scipy.optimize 모듈의 fmin_bfgs() 함수로 최소값을 구한다.
from scipy.optimize import fmin_bfgs
$f(x)= x^2 + 10 \sin (x)$ 의 최소값을 구한다.
파이썬 함수로 정의하고, 함수의 형태를 그래프로 그려본다.
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의 시작점을 넘겨준다.
fmin_bfgs( f, 0 )
Optimization terminated successfully.
Current function value: -7.945823
Iterations: 5
Function evaluations: 12
Gradient evaluations: 6
array([-1.30644012])
fmin_bfgs() 함수는 국부적인 최소점을 찾는 기능만 있으므로, 초기값에 따라 결과가 달라질 수 있다.
만약, x = 3 을 시작점으로 하면
fmin_bfgs( f, 3 )
Optimization terminated successfully.
Current function value: 8.315586
Iterations: 6
Function evaluations: 14
Gradient evaluations: 7
array([3.83746709])
가까운 국부 최소점을 찾게된다.
수치 적분¶
scipy.integrate 모듈의 quad() 함수를 이용하여 수치 적분을 계산한다.
$$ \int _0 ^1 x^2 dx $$
from scipy.integrate import quad
fx = lambda x: x**2
quad( fx, 0, 1 )
(0.33333333333333337, 3.700743415417189e-15)
첫번째 반환값이 수치 적분으로 얻어진 결과이며, 두번째의 값은 추정 오차이다.
피적분함수를 함수로서 정의하여 quad() 함수의 인자로 넘겨줄 수 있으며, 수학식을 표현할때 numpy 함수를 사용한다.
$$ \int _0 ^{\infty} e^{-x} dx $$
def f(x) :
return np.exp(-x)
F, err = quad( f, 0, np.inf )
print( F )
1.0000000000000002
선형 회귀법 (linear regression): 최소제곱법¶
scipy.stats 모듈의 linregress( x, y ) 함수를 사용하여 회기 직선의 기울기와 절편을 구한다.
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
회귀 직선과 데이터를 그려서 비교한다.
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 알고리즘¶
다음 그림에 주어진 데이터를 잘 표현하는 수식을 비선형 회귀법으로 얻으려고 한다.
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')
[<matplotlib.lines.Line2D at 0x285ffbb14f0>]
Levenberg-Marquardt 알고리즘을 이용하여 최적의 곡선을 구한다.
scipy.optimize 모듈의 curve_fit() 함수를 이용한다.
원하는 형태의 곡선에 대한 수식을 함수로 정의한다. 함수의 수식을 표현할때, numpy 수학 함수를 사용한다.
$$ y = a x + b { e }^{ -c{ x }^{ 2 } } $$
계수 $a, b, c$ 는 최적화에 의해 결정되어질 값이다.
# 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 가 최적화된 매개변수에 대한 리스트이다.
p0ini = [0, 0, 0] # initial values of parameters for fitting function
popt, pcov = curve_fit( func, x, y, p0 = p0ini )
print( popt )
[0.60378896 5.26885959 0.14667054]
구해진 매개변수로 최적화된 곡선을 예측하고, 원래의 주어진 데이터와 비교한다.
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() 함수를 이용하여 상미분방정식의 초기값 문제에 대한 수치해를 구할 수 있다.
from scipy.integrate import solve_ivp
일계 미분 방정식의 수치해를 구하는 경우를 고려해 보자.
$$ \frac {dy} {dt} = - 0.5 y $$도함수를 함수 형식으로 정의한다.
def dydt(t, y) :
return -0.5 * y
초기값 $y(0)=1$ 에 대하여, $t=10$ 까지 수치해를 구한다.
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$ 는 각각 다음과 같이 얻어진다.
sol.t
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. ])
sol.y
array([[1. , 0.77877121, 0.60652683, 0.47231798, 0.36766997,
0.28642722, 0.22325718, 0.17382285, 0.13534555, 0.10547428,
0.08216671, 0.0639639 , 0.04982546, 0.0388334 , 0.03023536,
0.02354182, 0.01834578, 0.01429275, 0.01112766, 0.00866797,
0.00675392]])
결과를 그래프로 그린다. sol.y 는 2차원 배열이므로 첫번째 행을 취한다.
plt.plot( sol.t, sol.y[0] )
plt.xlabel('t')
plt.ylabel('y')
plt.show()
이계 미분방정식¶
다음의 이계 미분방정식을 고려한다.
이계 미분방정식은 연립 일계 미분방정식으로 변환할 수 있다. $v (t) = y'(t)$ 로 두면
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) )
plt.plot( sol.t, sol.y[0] )
plt.xlabel('t')
plt.ylabel('y')
plt.show()
updated in 2025.09
댓글 없음:
댓글 쓰기