2 数値積分

Download Report

Transcript 2 数値積分

Slide 1

8.2 数値積分
(1)どんなときに数値積分を行うのか?
• 関数の形が分かっていないとき
(実際の測定データや観測データ)
• 解析的に不定積分を求めることが
不可能/困難なとき
工学的には,
不定積分の関数形が分からなくても,
多少の誤差があっても
数値的な定積分を求めれば良いことが多い。


Slide 2

(2)定積分の定義


b

のとき  a

f ( x ) dx  F ( x )

f ( x ) dx  F ( b )  F ( a )

または  F ( x ) ba

いま, f ( x ) や F ( x ) を求めることができないケースであるから,

次のような定義を使う。


ここで,

b
a

n

f ( x ) dx  lim

n 



f ( i )  x i

i 1

n

 0  a ,  n  b ,   xi  b  a
i 1

等間隔で分割することにすれば,  x i 

ba
n


Slide 3

定積分の定義の図
y

ξi
・・・

a

Δxi

b

x


Slide 4

ただし,nを無限大にできない。
• したがって,以下の式で我慢する
(多少の誤差は無視する)
n



f (a  i  x)x

i 1

ただし,

b  a  n  x


Slide 5

例題
• 次の定積分の近似値を求めよ。



0 . 35

sin x

0

e

x

dx

ただし,近似式は以下のとおりとする。
n


i 1

sin( a  i   x )
e

( a  i  x )

x

ここで

x 

ba
n


Slide 6

操作方法(1)
1行目に見出し,「ΔX」,「X」,「Y」,「単純積分」を設定
A2にΔXの値(ここでは0.05),
B2にXの初期値(ここでは0)を設定
B3に“=”をキーインし,B2を選択
ついで“+”をキーインし,A2を選択


最後に[ENTER]キーを押す

③,④,⑤


Slide 7

操作方法(2)
• B3の「=B2+A2」の「A2」を絶対指定,
すなわち「A$2」に変更して,[ENTER]キーを押す
$入力

[ENTER]キー


Slide 8

操作方法(3)
• B3を選択して,右ボタンクリック
• 「コピー」を選択する。


Slide 9

操作方法(4)
• コピー先を選択して,右ボタンクリック
• 「貼付け」を選択する。


Slide 10

操作方法(5)
• B列の小数点以下桁数を合わせ,
• 無駄な行を削除(0.00から0.35までを有効)


Slide 11

操作方法(6)
• C2を選択し,「=sin(B2)/exp(B2)」を入力し,
• 最後に[ENTER]キーを押す


Slide 12

操作方法(7)
• C2を,C3からC9にコピー
C2で右ボタンクリック

コピー先を選択して
右ボタンクリック


Slide 13

操作方法(8)
• D2を「=C2」とする


Slide 14

操作方法(9)
• D3を「=D2+C3*A$2」とする


Slide 15

操作方法(10)
• D3をD4からD9にコピー


Slide 16

操作方法(11)
• 小数点以下桁数を揃える。
• 以下の結果から約0.05であることが分かる。


Slide 17

nを大きく(ΔXを小さく)すると
精度が良くなる


Slide 18

(3)台形の公式
nを変えないで精度を良くするひとつの方法


Slide 19

次のように定義しよう


Slide 20

結果例


Slide 21

トーマス・シンプソンってどんな人?
(Thomas Simpson, 1710.8.20-1761.5.14)

イギリスのライチェスターシャー、マーケット・ボスワース生まれ
ウーリッジ陸軍大学校教授。
ロンドン王立協会会員(1746)。
[業績]
数値積分のシンプソン則
ド・モアヴル流の確率論
πの計算.


Slide 22

(4)シンプソンの公式:放物線で近似する方法(1)
点 ( x 2 i 1 , y 2 i 1 ), ( x 2 i , y 2 i ), ( x 2 i 1 , y 2 i 1 ) に対して,次の式で近似する。
y  px  qx  r
2

x 2 i 1   h , x 2 i  0 , x 2 i  1  h

とすると

y 2 i 1  ph  qh  r
2

y

y 2i  r
y 2 i  1  ph  qh  r
2

これらを解いて
P 

1
2h

q 

1
2h

r  y 2i

2

 y 2 i 1  2 y 2 i 

 y 2 i  1  y 2 i 1 

y 2 i 1 

x 2 i 1

x2i

x 2 i 1

x


Slide 23

放物線で近似する方法(2)
したがって
x 2 i 1

1
1

2
3
2
(
px

qx

r
)
dx

px

qx

rx
x 2 i 1
3

2

 x 2 i 1
x 2 i 1

1
1
1
  1

3
2
3
2
  ph  qh  rh     ph  qh  rh 
2
2
3
  3

2 1
 3


 ph  2 rh  
y

2
y

y
h  2 hy 2 i
2 i 1
2i
2 i 1 
2
3
3  2h

2



h
3

3

 y 2 i 1  4 y 2 i 

y 2 i 1 


Slide 24

放物線で近似する方法(3)
n



i 1

Si 

h
3

 y 0  4 ( y1 

y 3  y 5    y 2 n 1 )  2 ( y 2  y 4  y 6    y 2 n  2 )  y 2 n 

n
n 1
h

  y 0  y 2 n  4  y 2 i 1  2  y 2 i 
3
i 1
i 1


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


Slide 25

比較するための式定義

Δ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


Slide 26

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


Slide 27

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


Slide 28

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


Slide 29

(5)ガウスの数値積分法
ルジャンドルの多項式
ルジャンドル(Legendre)の多項式は次の漸化式で定義される
p0 ( x)  1
p1 ( x )  x
Pn ( x ) 

2n  1
n

xPn 1 ( x ) 

n 1
n

Pn  2 ( x )

p0 ( x)  1
p1 ( x )  x
P2 ( x )  ( 3 x  1) 2
2

P3 ( x )  ( 5 x  3 x ) 2
3

P4 ( x )  ( 35 x  30 x  3 ) 8
4

2

P5 ( x )  ( 63 x  70 x  15 x ) 8
5

3

( x  2)


Slide 30

ガウス・ルジャンドル多項式の性質
①次数が異なる2つのルジャンドル関数は相関がない。したがっ
て,次の式が成立する。



1

1

Pi ( x ) P j ( x ) dx  0

(i  j )

②n 次のルジャンドル多項式 Pn ( x ) は,n -1次の一般的な
多項式 Q ( x ) と相関がない。したがって,次の式が成立する。



1

1

Pi ( x ) Q ( x ) dx  0

(i  j )

③ルジャンドル方程式の根は,すべて実数である。


Slide 31

ちょっと横道にそれて
ルジャンドルってどんな人?
フランスの数学者。
主な功績:統計学,数論,代数学,解析学
ルジャンドルの仕事が元になって,
アーベルの楕円関数論,
ガウスによる統計学や数論の研究に
受け継がれている。

1752/9/18 - 1833/1/10

Adrien-Marie Legendre


Slide 32

ガウス・ルジャンドルの公式
数式が分かっている場合,
分割点は最適な結果が得られるように決定すればよい。
この分割点を積分点またはガウス点という。
基本式(ガウス・ルジャンドルの公式)



1

1

n

f ( x ) dx 

C
i0

i

f (t i )


Slide 33

積分区間の拡張
積分区間は t   1, 1 だが,次の変換を行うことで,
任意の積分区間に拡張できる。
x

ab



2



b

a

f ( x)d x 

ba
2

ba

t

2



1

1

ab ba
f

2
 2


t  dt



Slide 34

ガウス・ルジャンドルの公式による数値積分
積分点としてルジャンドル多項式が
0(根)となるような点を選択し,
右表の係数を用いて計算する。

xi 



b

a

ab



2

f ( x ) dx 

ba
2
ab
2

ti
n

C
i0

i

f ( xi )

n

i




2




3 2
3

4 2
3
4

2
5 3
4
5



ti

Ci

13




13
3 5

0
3 5

-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


Slide 35

VBAによる記述①
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)


Slide 36

VBAによる記述②
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


Slide 37

二重積分への拡張(1)
積分領域が長方形や直方体であるような二重積分や三重積分に
拡張することができる。
1

 

二重積分の式

1

1 1

二重積分 V 

d

 
c

n

f ( u , v ) dudv 

b

a

f ( x , y ) d xdy

CC
i

j0 i0

を考える。

積分区間が t   1, 1 になるように変換する
x

ba
2

y

d c
2

u  1   a
u  1   c

n

dx 

ba

du

2
dy 

d c
2

dv

j

f (ti , t j )


Slide 38

二重積分への拡張(2)
以上から次のように変換できる。
V 


1

 

1

1 1

g (u , v ) f (u , v )

( b  a )( d  c )

dudv

4

( b  a )( d  c )
4

n

n

CC
i

j

g (u i , v j )

j0 i0

d c
ba

g (u , v )  f 
( u  1)  a ,
( v  1)  c 
2
 2


以下,三重積分も同様に変換できる。


Slide 39

演習
(1)次の積分について台形の公式で数値積分を行うプログラム
を作成せよ。ただし,分割数は20とする。
S 



2

1

1
x

2

dx

(2)上記(1)の積分についてシンプソンの公式で数値積分を行う
プログラムを作成せよ。同様に分割数は20とする。
(3)上記(1)の積分についてガウスルジャンドルの公式を用いて
数値積分を行うプログラムを作成せよ。ただし,N=3とする。


Slide 40

(4)翼の長さに沿って均一に浮力が加わっている時の
翼の変位は,次の式で表現される。
 



Mx

L

M 

EI

0

w :浮力

dx ,

wx

演習

2

2

E :翼の弾性係数

I

:翼の慣性モーメント

I が次の式で多項式で与えられるものとして,翼の変位を求めよ。
I ( x )  3  10

9

x  5  10
4

ただし E  2  10 kg / m ,
9

2

6

x  3  10
3

4

x  8  10
2

w  3  10 kg / m ,
3

3

x  0 .1

L  12 m

とする。
L