情報処理概論
今日の授業内容(第9回)
情報処理概論(第9回)
„数値計算
‹定積分の計算
„台形近似法
„シンプソン法
„SSL Ⅱ(科学用サブルーチンライブラリ)の利用
情報基盤センター 井上 仁
台形近似法(1)
定積分
y = f (x)
h = (b-a)/m
区間[xi,xi+1]の積分siを
台形で近似
f(x)
oa
∫
b
a
b
Δx
n
f ( x)dx =　lim ∑
n →∞
k =1
b−a
=　f ( x )Δx Δx　n
k
∫
a
b
x0
m等分
xi
xm
xk= a + kΔx
f (x) の原始関数を F (x)
b
a
f(xi+2)
f(xi+1)
f(xi)
si si+1
(i = 0,…,m)
f ( x)dx　= [　F ( x)]　a　= F (b) − F (a)
b
台形近似法(2)
Xi = a + I*h
xi
xi+1 xi+2
台形近似法(3)
si = (h/2){f(xi)+f(xi+1)}
S = s0 + s1 + ・・・ + sm-1
= (h/2)[{f(x0)+f(x1)}+
(h/2){f(x1)+f(x2)}+
･･･
(h/2){f(xm-1)+f(xm)}]
h = (b-a)/m
s1 = 0.0D0
s0
s1
sm-1
= (h/2)[f(a)+f(b)+２{f(x1)+f(x2)+・・・+ｆ(xm-1)}]
DO i =1, m-1
s1 = s1 + f(a+i*h)
END DO
s = h*(f(a)+f(b)+2*s1)/2
DO構文を用いて計算
第9回
1
情報処理概論
シンプソン法(1)
台形近似法(4)
h = (b-a)/2m
PROGRAM trapezoid
・・・・・・・・・・・・
END PROGRAM trapezoid
FUNCTION f(x)
REAL :: f,x
f = 4/(1+x**2)
END FUNCTION f
区間[x2i,x2i+2]の積分siを
(x2i,f(x2i)),(x2i+1,f(x2i+1)),
(x2i+2,f(x2i+2))を通る
二次式の積分で近似
f(x)
a
a+ih
b
f(x2i+1)
2m等分
x0
xi
x2m
f(x2i)
(i = 0,…,2m)
シンプソン法(2)
x2i
si = (h/3)[f(x2i) + 4f(x2i+1) + f(x2i+2)]
si =
f(x2i+2)
(h/3)[f(x2i) + 4f(x2i+1) + f(x2i+2)]
si
x2i+1 x2i+2
シンプソン法(3)
f(x2i+1)
f(x2i)
x2i
f(x2i+2)
x2i+1 x2i+2
S = s0 + s1 + ・・・ + sm-1
= (h/3)[f(x0)+４ｆ（x1)+f(x2)]
＋(h/3)[f(x2)+４ｆ（x3)+f(x4)]
・・・・・・
＋(h/3)[f(x2m-2)+４ｆ（x2m-1)+f(x2m)]
S = s0 + s1 + ・・・ + sm-1
f(x)
s0a
x0
si
sm-1
xi
b
x2m
(i = 0,…,2m)
シンプソン法(4)
h = (b-a)/(2*m)
s1 = f(a+h)
s2 = 0.0D0
f(a) f(b)
= (h/3)[f(x0)+f(x2m)
+4{f(x1)+f(x3)+f(x5)+・・・+ｆ(x2m-1)} S1
+2{
f(x2)+f(x4)+・・・+f(x2m-2)}] S2
SSLⅡの利用(simpson2.f90)
„科学用サブルーチンライブラリ
CALL SIMP2(a, b, f, eps, s, icon)
DO i =1, m-1
s1 = s1 + f(a+(2*i+1)*h)
s2 = s2 + f(a+2*i*h)
END DO
s = h*(f(a)+f(b)+4*s1+2*s2)/3
第9回
‹関数f(x)の区間[a,b]での積分値をシンプソン法で
求める
‹epsは(絶対)誤差で、s0を真の積分値としたとき、
｜s - s0 ｜≦ eps を満たす s を求める
2
情報処理概論
コンパイルの流れ
ライブラリのリンク
% frt prog.f90
prog.f90
% frt prog.f90 -lfssl2
(ライブラリfssl2をリンクする)
Fortran90ソースファイル
prog.f90
コンパイル(翻訳)
prog.o
オブジェクト
ファイル
ライブラリ
prog.o
標準ライブラリ
fssl2
リンク(結合編集)
a.out
実行可能ファイル
a.out
演習9
次回の予定(第10回)
„台形近似法とシンプソン法、SSL Ⅱを用いた方
法で、計算結果がどうなるか確認しなさい。
„f(x)=4/(1+x**2)とし、積分区間は[0,1]で求
めなさい。
1
4
∫ 1+ x
0
第9回
„数値計算
‹方程式の解
„ニュートン法
„丸め誤差と桁落ち
dx　= ４[arctan( x)]0
1
2
3

付表5 力タログドプロシジャー覧

PowerPoint プレゼンテーション

[3分プレゼン]シュレーディンガー方程式の数値解