數值微分

 

4. 誤差還不賴

經過先前的介紹與分析,我們知道數值微分的限制以及如何透過分析尋找到最佳的數值解。我們可以用不同的差分公式來計算數值微分,但對應係數的可能性實在太多樣了,下列程式可以幫我們計算數值微分的對應係數。

 

 

 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
import numpy as npimport numpy.linalg as ladef difapx(N,points):
    """    difapx(N,points) will get the difference approximation for the Nth derivative    N: the order of derivative    points: the inteval of points    c: the output coefficient    err: the coefficient of error term    eoh: the order of error term     For example, difapx(1,[0,-2]) : derivative based on {f_0, f_{-1}, f_{-2}}    """
    l = max(points)
    L = abs(points[1]-points[0])+1
    if L < N+1:
        print 'More points are needed!'
        return
    A = np.zeros([L+2,L])
    A = np.mat(A)
    for n in range(1,L+1):
        A[0,n-1] = 1;
        for m in range(2,L+3):
            A[m-1,n-1] = A[m-2,n-1]*float(l)/(m-1)
        l -= 1
    b = np.zeros([L,1])
    b[N,0] = 1
    a = A[0:L,:]
    c = la.solve(a,b)
    err = float(A[L,:]*c)
    eoh = L-N #order of error term
    if abs(err)< 10**(-15):
        err = float(A[L+1,:]*c);
        eoh = L-N + 1
    c = c.transpose()
    c = c[0]
    if points[0]<points[1]:
        c = list(c)
        c.reverse()
        c = np.array(c)
    return c, err, eoh

使用方法是

>>> print difapx.__doc__

    difapx(N,points) will get the difference approximation for the Nth derivative

    N: the order of derivative

    points: the inteval of points

    c: the output coefficient

    err: the coefficient of error term

    eoh: the order of error term

    For example, c, err, eoh = difapx(1,[0,-2]) : derivative based on {f_0, f_{-1}, f_{-2}}

>>> c, err, eoh = difapx(1,[0,-2])

>>> c

array([ 1.5, -2. ,  0.5])

>>> err

-0.3333333333333333

>>> eoh

2

我們將以下程式儲存成diff_exam.py。透過這個程式要驗證我們剛才的分析是否與實際的數值解狀況相符。

 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
#!/usr/bin/python

import sys

sys.path.append('/home/glophy/python/toolbox')
from pynum_int import *
import numpy as np
#This example compute the relationship between step size and the error

def test_f(x):
    return np.sin(x)

def test_df(x):
    return np.cos(x)
st = 5 
et = 12 
H = [10**(-i) for i in range(st,et)]


result = list()
for h in H:
    r = np.pi/4
    x = [r,r+h]
    x = np.array(x)
    y = test_f(x)
    dy = test_df(x)
    c,err,eoh = difapx(1,[0,1])
  
    result.append(abs(sum(c*y)/h-dy[0]))

from pylab import *
fig = figure(1)
plot(range(st,et),result)
fig.show()
title('error between dy and num_dy')
xlabel(r'$h = 10^i$: i ')
fig.savefig('diff_error.png')

將此檔案變成可執行檔在linux 環境下執行,或將第一行刪除後直接在Python IDLE 中執行,可以製作出誤差與step size 的關係圖 。

difference error

 

我們看到,這結果與我們的分析是一致的,並且在 \(h_i\) 控制得宜時,數值解與真解的誤差還不賴。希望透過這樣的展示能讓同學體會到數值分析的精神。


練習題

  1. 請寫一個 Python 程式,給定函數 \(y(x) = sin(x)\),\(x_0 = \frac{\pi}{4}\),\(h = 1/8, 1/16,\dots, 1/2^{10} \),同時比較 forward , backward 和 central difference 的差別。

  2. 在固定間距 \(h\) 和階數時,forward difference 和 backward difference 有相同的truncation error,試問在何種場合你會選擇 backward difference method 甚於 forward difference method,理由為何?