Excelで定積分・微分方程式

８．２数値積分
（１）どんなときに数値積分を行うのか？
• 関数の形が分かっていないとき
（実際の測定データや観測データ）
• 解析的に不定積分を求めることが
不可能／困難なとき
工学的には，
不定積分の関数形が分からなくても，
多少の誤差があっても
数値的な定積分を求めれば良いことが多い。
（２）定積分の定義
 f ( x)dx  F ( x)
のとき

b
a
f ( x)dx  F (b)  F (a) または
F ( x)
b
a
いま， f (x) や F (x) を求めることができないケースであるから，
次のような定義を使う。

ここで，
b
a
n
f ( x)dx  lim  f (i )xi
n 
i 1
n
 0  a,  n  b,  xi  b  a
i 1
ba
等間隔で分割することにすれば， xi 
n
定積分の定義の図
y
ξi
・・・
a
Δxi
b
x
ただし，ｎを無限大にできない。
• したがって，以下の式で我慢する
（多少の誤差は無視する）
n
 f (a  i  x)x
i 1
ただし，
b  a  n  x
例題
• 次の定積分の近似値を求めよ。

0.35
0
sin x
dx
x
e
ただし，近似式は以下のとおりとする。
sin(a  i  x)
x

( a  ix )
e
i 1
n
ここで
ba
x 
n
操作方法（１）
1行目に見出し，「ΔＸ」，「Ｘ」，「Ｙ」，「単純積分」を設定
Ａ２にΔＸの値（ここでは0.05），
Ｂ２にＸの初期値（ここでは0）を設定
Ｂ３に“=”をキーインし，Ｂ２を選択
ついで“＋”をキーインし，Ａ２を選択
②
①
最後に[ENTER]キーを押す
③，④，⑤
操作方法（２）
• Ｂ３の「＝Ｂ２＋Ａ２」の「Ａ２」を絶対指定，
すなわち「Ａ＄２」に変更して，[ＥＮＴＥＲ]キーを押す
＄入力
[ENTER]キー
操作方法（３）
• Ｂ３を選択して，右ボタンクリック
• 「コピー」を選択する。
操作方法（４）
• コピー先を選択して，右ボタンクリック
• 「貼付け」を選択する。
操作方法（５）
• Ｂ列の小数点以下桁数を合わせ，
• 無駄な行を削除（0.00から0.35までを有効）
操作方法（６）
• Ｃ２を選択し，「=sin(B2)/exp(B2)」を入力し，
• 最後に[ENTER]キーを押す
操作方法（７）
• Ｃ２を，Ｃ３からＣ９にコピー
Ｃ２で右ボタンクリック
コピー先を選択して
右ボタンクリック
操作方法（８）
• Ｄ２を「＝Ｃ２」とする
操作方法（９）
• Ｄ３を「＝Ｄ２＋Ｃ３＊Ａ＄２」とする
操作方法（１０）
• Ｄ３をＤ４からＤ９にコピー
操作方法（１１）
• 小数点以下桁数を揃える。
• 以下の結果から約0.05であることが分かる。
ｎを大きく（ΔＸを小さく）すると
精度が良くなる
（３）台形の公式
ｎを変えないで精度を良くするひとつの方法
次のように定義しよう
結果例
トーマス・シンプソンってどんな人？
(Thomas Simpson, 1710.8.20-1761.5.14)
イギリスのライチェスターシャー、マーケット・ボスワース生まれ
ウーリッジ陸軍大学校教授。
ロンドン王立協会会員(1746)。
［業績］
数値積分のシンプソン則
ド・モアヴル流の確率論
πの計算．
（４）シンプソンの公式：放物線で近似する方法（１）
点 ( x2i 1, y2i 1 ), ( x2i , y2i ), ( x2i1, y2i 1 ) に対して，次の式で近似する。
y  px2  qx  r
x2i 1  h, x2i  0, x2i 1  h とすると
y2i 1  ph2  qh  r
y
y2i  r
y2i 1  ph2  qh  r
これらを解いて
1
 y 2i 12 y2i  y2i 1 
2
2h
1
 y 2i 1 y2i 1 
q
2h
r  y2i
P
x2i 1
x2i
x2i 1
x
放物線で近似する方法（２）
したがって
x
2 i 1
1 3 1 2

2
x2i1 ( px  qx  r )dx   3 px  2 qx  rx  x
2 i 1
x2 i 1
1
1
1
  1

  ph3  qh2  rh     ph3  qh2  rh 
2
2
3
  3

2 3
2 1

 ph  2rh   2  y 2i 12 y2i  y2i 1 h 3  2hy2i
3
3  2h

h
  y 2i 14 y2i  y2i 1 
3
放物線で近似する方法（３）
n
h
Si  y0  4( y1  y3  y5    y2 n 1 )  2( y2  y4  y6    y2 n  2 )  y2 n 

3
i 1
n
n 1
h

  y0  y2 n  4 y2i 1  2 y2i 
3
i 1
i 1

区間１
f(x)
区間2
放物線
近似
y0
放物線
近似
y1
y2
y3
S1
y4
S2
x
x0
区間1のt
－h
x2
x3
x4
0
h
区間2のt
－h
0
h
x1
比較するための式定義
ΔX
0.05
X
0.00
0.05
0.10
0.15
0.20
0.25
0.30
0.35
Y
0.000000
0.047542
0.090333
0.128623
0.162657
0.192678
0.218927
0.241636
単純積分台形の公式
0.000000
0.000000
0.002377
0.001189
0.006894
0.004635
0.013325
0.010109
0.021458
0.017391
0.031092
0.026275
0.042038
0.036565
0.054120
0.048079
合計
Y0, Y2n
0.000000
奇数項
偶数項
シンプソンの公式
0.047542
0.090333
0.128623
0.162657
0.192678
0.218927
0.241636
0.241636 0.368843
0.471916
0.044347
VBAによる表現①ボタンのClickイベントハンドラ，関数値設定
Private Const B_max = 1
Private Const B_min = 0
Sub ボタン2_Click()
MsgBox "単純積分
= " & 単純積分()
MsgBox "台形の公式 = " & 台形の公式()
MsgBox "シンプソンの公式 = " & シンプソンの公式()
End Sub
Function Fx(X) As Double
Fx = 1 / (1 + X)
End Function
VBAによる表現③単純数値積分，台形の公式
Function 単純積分() As Double
DX = (B_max - B_min) / 20
S = 0: X = B_min
For i = 1 To 20
S = S + DX * Fx(X)
X = X + DX
Next
単純積分 = S
End Function
Function 台形の公式() As Double
DX = (B_max - Bmin) / 20
S = 0: X1 = B_min
For i = 1 To 19
X2 = X1 + DX
S = S + DX * (Fx(X1) + Fx(X2)) / 2
X1 = X2
Next
台形の公式 = S
End Function
VBAによる表現④シンプソンの公式
Function シンプソンの公式() As Double
DX = (B_max - B_min) / 20
S1 = Fx(B_min) + Fx(B_max)
X1 = B_min + DX
S2 = 0
For i = 1 To 10
S2 = S2 + Fx(X1)
X1 = X1 + DX * 2
Next
S3 = 0
X1 = B_min
For i = 1 To 9
X1 = X1 + DX * 2
S3 = S3 + Fx(X1)
Next
S = DX * (S1 + 4 * S2 + 2 * S3) / 3
シンプソンの公式 = S
End Function
（５）ガウスの数値積分法
ルジャンドルの多項式
ルジャンドル（Legendre）の多項式は次の漸化式で定義される
p0 ( x)  1
p1 ( x)  x
Pn ( x) 
2n  1
n 1
xPn 1 ( x) 
Pn  2 ( x) ( x  2)
n
n
p0 ( x)  1
p1 ( x)  x
P2 ( x)  (3x 2  1) 2
P3 ( x)  (5 x 3  3x) 2
P4 ( x)  (35x 4  30x 2  3) 8
P5 ( x)  (63x 5  70x 3  15x) 8
ガウス・ルジャンドル多項式の性質
①次数が異なる２つのルジャンドル関数は相関がない。したがっ
て，次の式が成立する。
1

Pi ( x) Pj ( x)dx  0 (i  j )
1
②n 次のルジャンドル多項式 Pn (x) は，n -１次の一般的な
多項式 Q (x) と相関がない。したがって，次の式が成立する。
1

1
Pi ( x)Q( x)dx  0 (i  j )
③ルジャンドル方程式の根は，すべて実数である。
ちょっと横道にそれて
ルジャンドルってどんな人？
フランスの数学者。
主な功績：統計学，数論，代数学，解析学
ルジャンドルの仕事が元になって，
アーベルの楕円関数論，
ガウスによる統計学や数論の研究に
受け継がれている。
1752/9/18 - 1833/1/10
Adrien-Marie Legendre
ガウス・ルジャンドルの公式
数式が分かっている場合，
分割点は最適な結果が得られるように決定すればよい。
この分割点を積分点またはガウス点という。
基本式（ガウス・ルジャンドルの公式）
1

1
n
f ( x)dx   Ci f (ti )
i 0
積分区間の拡張
積分区間は t   1, 1だが，次の変換を行うことで，
任意の積分区間に拡張できる。
ab ba
x

t
2
2

b
a
ba 1  ab ba 
f ( x) d x 
f

t dt

2 1  2
2 
ガウス・ルジャンドルの公式による数値積分
積分点としてルジャンドル多項式が
0（根）となるような点を選択し，
右表の係数を用いて計算する。
ab ba
xi 

ti
2
2
b
ab n
Ci f ( xi )
a f ( x)dx  2 
i 0
n i
ti
Ci
１
2
 13
13
１
１
１
３ 2
3
１
４ 2
3
4
１
2
５ 3
4
5
 35
0
２
35
-0.8611363116
-0.3399810436
0.3399810436
0.8611363116
-0.9061798459
-0.5384693101
0
0.5384693101
0.9061798459
5/9
8/9
5/9
0.3478548451
0.6521451549
0.6521451549
0.3478548451
0.2369268851
0.4786286705
0.5688888889
0.4786286705
0.2369268851
ＶＢＡによる記述①
Private Const M = 5
Private Const B_max = 1
Private Const B_min = 0
Private Ci(M) As Double
Private ti(M) As Double
Sub 積分係数設定()
ti(1) = -0.9061798459:
ti(2) = -0.5384693101:
ti(3) = 0:
ti(4) = -ti(2):
ti(5) = -ti(1):
End Sub
Function Fx(X) As Double
Fx = 1 / (1 + X)
End Function
Ci(1)
Ci(2)
Ci(3)
Ci(4)
Ci(5)
=
=
=
=
=
0.2369268851
0.4786286705
0.5688888889
Ci(2)
Ci(1)
ＶＢＡによる記述②
Function ガウス積分() As Double
M1 = (B_max - B_min) / 2
M2 = (B_min + B_max) / 2
S = 0
For i = 1 To M
X = M1 * ti(i) + M2
S = S + Ci(i) * Fx(X)
Next
ガウス積分 = S * M1
End Function
Sub ボタン1_Click()
積分係数設定
MsgBox ガウス積分()
End Sub
二重積分への拡張（１）
積分領域が長方形や直方体であるような二重積分や三重積分に
拡張することができる。
1

二重積分の式
二重積分 V  
1
1 1
d
c

b
a
n
f (u, v)dudv  Ci C j f (ti , t j )
j 0 i 0
f ( x, y)d xdy を考える。
積分区間が t   1, 1 になるように変換する
ba
u  1  a
x
2
d c
u  1  c
y
2
n
ba
dx 
du
2
d c
dy 
dv
2
二重積分への拡張（２）
以上から次のように変換できる。
(b  a)(d  c)
V    g (u, v) f (u, v)
dudv
1 1
4
(b  a)(d  c) n n

Ci C j g (ui , v j )

4
j 0 i 0
1
1
d c
ba

g (u, v)  f 
(u  1)  a,
(v  1)  c 
2
 2

以下，三重積分も同様に変換できる。
演習
（１）次の積分について台形の公式で数値積分を行うプログラム
を作成せよ。ただし，分割数は２０とする。
S
2
1
1
dx
2
x
（２）上記（１）の積分についてシンプソンの公式で数値積分を行う
プログラムを作成せよ。同様に分割数は２０とする。
（３）上記（１）の積分についてガウスルジャンドルの公式を用いて
数値積分を行うプログラムを作成せよ。ただし，Ｎ＝３とする。
（４）翼の長さに沿って均一に浮力が加わっている時の
翼の変位は，次の式で表現される。
L Mx
wx 2
 
dx, M 
0 EI
2
w ：浮力
E ：翼の弾性係数
演習
I ：翼の慣性モーメント
I が次の式で多項式で与えられるものとして，翼の変位を求めよ。
I ( x)  3109 x 4  5 106 x3  3104 x 2  8 103 x  0.1
9
2
3
E

2

10
kg
/
m
,
w

3

10
kg / m, L  12m
ただし
とする。
L


Download Report