Numeric Analysis¶
from math import *
Overflow¶
수치 값이 대략 $10^{308}$ 이 넘으면 오버플로우(overflow) 에러가 발생한다.
1.79e308
위의 값이 컴퓨터가 수치적으로 허용하는 최대값이다. 더 큰 값은 무한대로 취급한다.
1.80e308
Underflow¶
1.23e-308
수치 값이 대략 $10^{-308}$ 보다 작으면 언더플로우(underflow)가 되어 0으로 처리된다.
1.0e-324
Roundoff¶
반올림(rounding)에 의하여 수치적인 오차가 생긴다.
round( pi, 2 )
컴퓨터는 중간 계산 과정에서는 반올림을 하지 않으므로 큰 문제는 없다.
Truncation error¶
유효 숫자의 갯수가 16자리를 초과하면, 그 뒷자리의 숫자는 잘려나간다. 이로 인하여 수치적 오차가 발생한다.
1.0e10 + 1
1.0e15 + 1
지금까지는 문제가 없다. 하지만 다음의 경우
1.0e16 + 1
아주 큰 숫자와 아주 작은 숫자는 덧셈이 불가능하게 되는 오류가 생기는 것이다.
Solution of Equations by Iteration¶
to find solutions of a single equation numerically
$$ f(x) = 0 $$
Fixed-Point Iteration¶
After transforming into the form $$ x = g(x) $$
We start with a guessed value $x_0$, and compute
$\qquad x_1 = g(x_0)$, $x_2 = g(x_1)$, $\cdots$
It is also called the method of successive substitution.
[ Example 1 ] $\,\,$ Find a solution of $\, f(x)=x^2-3x+1=0 $
$\qquad$ Rewrite the equation as $$ x = \frac 1 3 (x^2+1)$$
# initial guess
x = 1
for i in range(10) :
xnew = 1/3 *(x**2+1)
x = xnew
print("%9.5f" % x )
# initial guess
x = 3
for i in range(10) :
xnew = 1/3 *(x**2+1)
x = xnew
print("%9.5f" % x )
Newton’s Method for Solving Equations ƒ(x) = 0¶
$$ x_{n+1} = x_n - \frac {f(x_n)}{f'(x_n)} $$
[ Example 4 ] $\,\,$ Find a solution of $\,\, 2 \sin x =x $
from scipy.optimize import newton
def func(x):
return 2*sin(x)-x # f(x) = 0
xsol = newton( func, 2 ) # newton(함수이름,초기값)
print("%9.5f" % xsol)
소스 코드
def func(x) :
return 2*sin(x)-x # f(x) = 0
def newton( x0 ) :
x = x0
dx = x * 1e-5
for i in range(100) : # maximum 100 iterations
f = func( x )
if abs( f ) < 1e-5 : break
dfdx = ( func( x+dx ) - f ) / dx
x += - f / dfdx
return x
xsol = newton( 2 )
print("%9.5f" % xsol )
Interpolation $\quad$ (내삽, 보간)¶
Linear interpolation¶
interpolation by the straight line through $ (x_0,f_0)$, $(x_1,f_1)$
[ Example 1 ] $\,\,$ Linear Lagrange Interpolation
$\qquad$ For $f(x)=\ln x$, compute $\ln 9.2$ from $\,\, \ln 9.0= 2.1972$ and $ \ln 9.5 = 2.2513$
def interpol( x, x0,f0, x1,f1 ) :
p = (x-x1)/(x0-x1)*f0 + (x-x0)/(x1-x0)*f1
return p
p = interpol( 9.2, 9.0, 2.1972, 9.5, 2.2513 )
print("approx = %9.5f" % p )
print("ln 9.2 = %9.5f (exact value)" % log(9.2) )
Numeric Integration¶
Rectangular rule¶
Piecewise constant approximation of $f(x)$
Trapezoidal rule¶
Piecewise linear approximation of $f(x)$
$$ h = \frac {b-a} n $$
[ Example 1 ] $\,\,$ Trapezoidal Rule
$\qquad \,\,$ Evaluate $\,\, J= \int _0 ^1 e^{-x^2} dx $
def f(x) :
return exp(-x**2)
n = 10
a = 0
b = 1
h = (b-a)/n
sum = f(a)/2 + f(b)/2
for i in range(1,n) :
xi = a + i*h
sum += f(xi)
J = h * sum
print("%9.6f" % J )
Scipy 라이브러리를 사용하여 정확한 적분값을 구해보면
from scipy.integrate import quad
J, err = quad( f, 0, 1 )
print("%9.6f" % J )
Simpson’s Rule¶
Piecewise quadratic approximation of $f(x)$
인접한 3개의 점들을 연결하는 2차 곡선 (포물선)을 만들고, 이러한 2차 곡선들을 적분함으로써 함수 $f(x)$의 정적분을 계산한다.
1차 곡선(직선)으로 가정하는 사다리꼴 공식(trapezoidal rule)보다 정확도가 더 높다.
$$ h = \frac {b-a} {2m} $$
where $2m$ is the number of subintervals, $n=2m$
[ Example 3 ] $\,\,$ Simpson’s Rule
$\qquad \,\,$ Evaluate $\,\, J= \int _0 ^1 e^{-x^2} dx \,\,$ with $\,n=2m=10$
def f(x) :
return exp(-x**2)
n = 10
a = 0
b = 1
h = (b-a)/n
s0 = f(a) + f(b)
s1 = 0
s2 = 0
for i in range(1,n) :
xi = a + i*h
if i%2 == 1 : # 홀수
s1 += f(xi)
else : # 짝수
s2 += f(xi)
J = h/3 *( s0 + 4*s1 + 2*s2 )
print("%9.6f" % J )
Only the last digit differs from the exact value of 0.746824
Adaptive Integration¶
스텝의 크기 $h$ 를 구간 마다 가변적으로 적용하여 적분 값의 오차를 줄이는 방법
Gauss Integration Formula (Gauss Quadrature)¶
미리 정해진 노드 $t_j$ 와 계수 (또는 가중치) $A_j$ 를 이용하여 적분 값을 계산한다.
$$ x = \frac 1 2 [-a(t-1) + b(t+1)]$$
변수 변환에 의하여
$$ \int _a ^b f(x)dx = \frac {b-a} 2 \int _{-1}^1 f[x(t)] dt \approx \frac {b-a} 2 \sum _{j=1} ^n A_j f_j $$
[ Example 7 ] $\,\,$ Gauss Integration Formula with $n=3$
$\qquad \,\,$ Evaluate $\,\, J= \int _0 ^1 e^{-x^2} dx \,\,$
def f(x) :
return exp(-x**2)
n = 3
t = [ -0.7745966692, 0, 0.7745966692]
A = [ 0.5555555556, 0.8888888889, 0.5555555556]
a = 0
b = 1
sum = 0
for i in range(n) :
xi = 1/2*(-a*(t[i]-1) + b*(t[i]+1) )
sum += A[i] * f(xi)
J = (b-a)/2 * sum
print("%9.6f" % J )
With $n=4$
def f(x) :
return exp(-x**2)
n = 4
t = [ -0.8611363116, -0.3399810436, 0.3399810436, 0.8611363116]
A = [ 0.3478548451, 0.6521451549, 0.6521451549, 0.3478548451]
a = 0
b = 1
sum = 0
for i in range(n) :
xi = 1/2*(-a*(t[i]-1) + b*(t[i]+1) )
sum += A[i] * f(xi)
J = (b-a)/2 * sum
print("%9.6f" % J )
The result is almost the same with the exact value of 0.746824
Numeric Differentiation¶
Differentiation is defined by
$$ f'(x) = \lim _{h \rightarrow 0} \frac {f(x+h)-f(x)} h $$
Numeric differentiation is approximated by the forward difference
$$ f'(x) \approx \frac {f(x+h)-f(x)} h \qquad \,\,\, \text {with a small } h $$or given by the backward difference
$$ f'(x) \approx \frac {f(x)-f(x-h)} h \qquad \,\,\, $$The central difference gives more accure approximation
$$ f'(x) \approx \frac {f(x+h)-f(x-h)} {2h} $$[ Example ] find the value of the derivative of $ e^{-x^2}$ at $x=1$.
def f(x) :
return exp(-x**2)
def dfdx(x) :
return -2*x*exp(-x**2)
x = 1
h = 0.01
# the forward difference
fp = ( f(x+h) - f(x) )/h
print("the forward difference f' = %9.6f" % fp )
# the backward difference
fp = ( f(x) - f(x-h) )/h
print("the backward difference f' = %9.6f" % fp )
# the central difference
fp = ( f(x+h) - f(x-h) )/(2*h)
print("the central difference f' = %9.6f" % fp )
# the exact value
print()
print("the exact value f' = %9.6f" % dfdx(x) )
The second-order derivative¶
The forward difference
$$ f''(x) \approx \frac {f(x)-2f(x+h)+f(x+2h)} {h^2} $$
The backward difference
$$ f''(x) \approx \frac {f(x-2h)-2f(x-h)+f(x)} {h^2} $$
The central difference
$$ f''(x) \approx \frac {f(x-h)-2f(x)+f(x+h)} {h^2} $$
댓글 없음:
댓글 쓰기