## 數值微分

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]

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

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

### 練習題

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，理由為何？