スライド タイトルなし

8.数値微分・積分・微分方程式
工学的問題においては
解析的に微分値や積分値を求めたり,
微分方程式を解くことが難しいケースも多い。
このような場合は数値的に解かざるを得ない。
8.1 数値微分
(1)差分法を用いる
関数
f (x) の x  a における微分係数 f (x ) の定義
f ( a  h)  f ( a )
f (a)  lim
h 0
h
h
: 刻み幅
差分近似
[前進差分近似]
[中心差分近似]
[後退差分近似]
f ( a  h)  f ( a )
f (a) 
h
f ( a  h)  f ( a  h)

f (a) 
2h
f ( a )  f ( a  h)
f (a) 
h
(2)差分法の誤差
①刻み幅が大きいと,
誤差が大きい(打切り誤差)。
②刻み誤差が小さ過ぎると,
丸め誤差による誤差が大きくなる。
①前進差分法の打切り誤差
前進差分では,以下のテーラ展開から導く
h2
h3
f ( x  h)  f ( x)  hf ( x) 
f ( x) 
f ( x)  O(h 4 )
2!
3!
x  a を代入して,整理すると。
f ( a  h)  f ( a ) h
h2
f (a) 
 f (a) 
f (a)  O(h 3 )
h
2!
3!
f ( a  h)  f ( a )

 O ( h)
h
h に比例した打切り誤差
②後退差分法の打切り誤差
後退差分では,以下のテーラ展開から導く
h2
h3
f ( x  h)  f ( x)  hf ( x) 
f ( x) 
f ( x)  O(h 4 )
2!
3!
x  a を代入して,整理すると。
f ( a )  f ( a  h) h
h2
f (a) 
 f (a) 
f (a)  O(h 3 )
h
2!
3!
f ( a  h)  f ( a )

 O ( h)
h
h に比例した打切り誤差
③中心差分法の打切り誤差
後退差分では,以下の2つのテーラ展開から導く
h2
h3
f ( x  h)  f ( x)  hf ( x) 
f ( x) 
f ( x)  O(h 4 )
2!
3!
h2
h3
f ( x  h)  f ( x)  hf ( x) 
f ( x) 
f ( x)  O(h 4 )
2!
3!
2h 3
 f ( x  h)  f ( x  h)  2hf ( x) 
f ( x)
3!
x  a を代入して,整理すると。
f ( a  h)  f ( a  h) 2h 2
f (a) 

f (a)  O(h 3 )
2h
3!
f ( a  h)  f ( a  h)

 O(h 2 )
2h
h 2 に比例した打切り誤差
前進差分,後退差分,中心差分をExcelで式定期
誤差の比較(1)
前進差分や後退差分に比べ,中心差分の誤差が少ない!!
誤差の比較(2)
微分値は,ほとんど中心差分と重なっているが…
細かくみると
誤差はある…
誤差の比較(3)
刻み幅が小さくなると誤差が減るが
刻み幅による誤差の比較
小さすぎると丸め誤差が大きくなる
誤
差
刻み幅 h
刻み幅による誤差データの生成プログラム(1)
Function FX(X) As Single
X2 = X * X: X3 = X2 * X: X4 = X3 * X
FX = X4 + X3 + X2 + X + 1
End Function
Function DFX(X) As Double
X2 = X * X: X3 = X2 * X
DFX = 4 * X3 + 3 * X2 + 2 * X + 1
End Function
誤差の影響を大きくするために,あえて関数値はSingle,
微分値はDoubleにしてある。
次のページの変数値も同様である。
[参考]
[参考]
刻み幅による誤差データの生成プログラム(2)
Sub ボタン1_Click()
Dim DF1 As Single: Dim DF2 As Single
With Worksheets("Sheet2")
h = 0.1
For i = 1 To 100
DF = DFX(1)
DF1 = (FX(1 + h) - FX(1)) / h
DF2 = (FX(1 + h) - FX(1 - h)) / (2 * h)
G1 = Abs(DF1 - DF)
G2 = Abs(DF2 - DF)
.Cells(i + 1, 1) = h
.Cells(i + 1, 2) = G1
.Cells(i + 1, 3) = G2
h = h / 1.2
Next
End With
End Sub
‘
‘
‘
‘
真値
前進差分
中心差分
誤差の計算
VBAによる表現 ①Excelシートの定義
A~D列はExcelでの式定義における定義と同じにする。
EFG列は,2行目以降は空白にして,
差分の結果を入れるセルとする。
VBAによる表現 ②データ宣言とデータ設定
Private DX As Double
Private Y() As Double
Private DY() As Double
Private N As Double
Sub データ設定()
With Worksheets("Sheet4")
i = 2
‘
入っているデータをカウントする。
Do While .Cells(i, 2) <> "": i = i + 1: Loop
N = i - 2
ReDim Y(N), DY(N)
DX = .Cells(2, 1)
For i = 1 To N
Y(i) = .Cells(i + 1, 3)
Next
End With
End Sub
VBAによる表現 ③結果設定とボタンのClickイベントハンドラ
Sub 結果設定(DY, Ist, N, ID)
With Worksheets("Sheet4")
For i = Ist To N
.Cells(i + 1, ID) = DY(i)
Next
End With
End Sub
Sub Sheet4_ボタン1_Click()
データ設定
前進差分 DX, Y, DY, N
結果設定 DY, 1, N - 1, 5
後退差分 DX, Y, DY, N
結果設定 DY, 2, N, 6
中心差分 DX, Y, DY, N
結果設定 DY, 2, N - 1, 7
End Sub
VBAによる表現 ④各差分の定義
Sub 前進差分(DX, Y, DY,
For i = 1 To N - 1
DY(i) = (Y(i + 1) Next
End Sub
Sub 後退差分(DX, Y, DY,
For i = 2 To N
DY(i) = (Y(i) - Y(i
Next
End Sub
Sub 中心差分(DX, Y, DY,
For i = 2 To N - 1
DY(i) = (Y(i + 1) Next
End Sub
N)
Y(i)) / DX
N)
- 1)) / DX
N)
Y(i - 1)) / (2 * DX)
(3)高次の差分法
f ( x  2h) をテーラ展開し
( 2h) 2
( 2 h) 3
f ( x  2h)  f ( x)  2hf ( x) 
f ( x) 
f ( x)  O(h 4 )
2!
3!
h
1 
( 2h) 3
4 






f ( x)   f ( x  2h)  f ( x)  2hf ( x) 
f ( x)  O(h ) 
2!
4h 
3!

(3)高次の差分法
h
f ( x ) を1次の式に代入し
2!
f ( a  h)  f ( a ) h
h2
f (a) 
 f (a) 
f (a)  O(h 3 )
h
2!
3!
4 f (a  h)  4 f (a)  f (a  2n)  f (a) f (a)


 O(h 2 )
4h
2
したがって,
f (a) f (a  h)  f (a) h
h2

 f (a) 
f (a)  O(h 3 )
2
h
2!
3!
4 f ( a  h)  3 f ( a )  f ( a  2n)

 O(h 2 )
4h
2次の前進差分
f (a) 
4 f ( a  h)  3 f ( a )  f ( a  2n)
 O(h 2 )
2h
(3)高次の差分法
同様に,2次の後退差分,中心差分も求めることができる。
2次の後退差分
f (a) 
 4 f ( a  h)  3 f ( a )  f ( a  2n)
 O(h 2 )
2h
2次の中心差分
f (a) 
f ( a  2h)  8 f ( a  h)  8 f ( a  h)  f ( a  2h)
 O(h 4 )
12 h
前進差分における1次と2次の誤差
同様に,2次の後退差分,中心差分でも確認してみよう。