Numeric Analysis 1 - 기초

Numeric Analysis¶

In [1]:

from math import *

Overflow¶

수치 값이 대략 $10^{308}$ 이 넘으면 오버플로우(overflow) 에러가 발생한다.

In [2]:

1.79e308

Out[2]:

1.79e+308

위의 값이 컴퓨터가 수치적으로 허용하는 최대값이다. 더 큰 값은 무한대로 취급한다.

In [3]:

1.80e308

Out[3]:

inf

Underflow¶

In [4]:

1.23e-308

Out[4]:

1.23e-308

수치 값이 대략 $10^{-308}$ 보다 작으면 언더플로우(underflow)가 되어 0으로 처리된다.

In [5]:

1.0e-324

Out[5]:

0.0

Roundoff¶

반올림(rounding)에 의하여 수치적인 오차가 생긴다.

In [6]:

round( pi, 2 )

Out[6]:

3.14

컴퓨터는 중간 계산 과정에서는 반올림을 하지 않으므로 큰 문제는 없다.

Truncation error¶

유효 숫자의 갯수가 16자리를 초과하면, 그 뒷자리의 숫자는 잘려나간다. 이로 인하여 수치적 오차가 발생한다.

In [7]:

1.0e10 + 1

Out[7]:

10000000001.0

In [8]:

1.0e15 + 1

Out[8]:

1000000000000001.0

지금까지는 문제가 없다. 하지만 다음의 경우

In [9]:

1.0e16 + 1

Out[9]:

1e+16

아주 큰 숫자와 아주 작은 숫자는 덧셈이 불가능하게 되는 오류가 생기는 것이다.

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)$$

In [10]:

# initial guess
x = 1  

for i in range(10) :
    xnew = 1/3 *(x**2+1)
    x = xnew
    print("%9.5f" % x )

  0.66667
  0.48148
  0.41061
  0.38953
  0.38391
  0.38246
  0.38209
  0.38200
  0.38197
  0.38197

In [11]:

# initial guess
x = 3  

for i in range(10) :
    xnew = 1/3 *(x**2+1)
    x = xnew
    print("%9.5f" % x )

  3.33333
  4.03704
  5.76589
 11.41516
 43.76863
638.89754
136063.68691
6171108965.24800
12694195286988072960.00000
53714197994730060054154112125389242368.00000

$$x = \frac 1 3 (x^2+1)$$

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$

In [12]:

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)

  1.89549

소스 코드

In [13]:

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 )

  1.89549

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$

In [14]:

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) )

approx =   2.21884
ln 9.2 =   2.21920  (exact value)

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$

In [15]:

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 )

 0.746211

Scipy 라이브러리를 사용하여 정확한 적분값을 구해보면

In [16]:

from scipy.integrate import quad

J, err = quad( f, 0, 1 )
print("%9.6f" % J )

 0.746824

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$

In [17]:

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 )

 0.746825

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 \,\,$

In [18]:

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 )

 0.746815

With $n=4$

In [19]:

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 )

 0.746824

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$.

In [20]:

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 forward  difference  f' = -0.732056
the backward difference  f' = -0.739413
the central  difference  f' = -0.735734

the exact value          f' = -0.735759

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}$$