1/12の講義資料

情報基礎演習Ⅱ
第12回演習
担当光武雄一
授業のWebsite: http://web.me.saga-u.ac.jp/~mitutake/info2/info2.html
メールアドレス： [email protected]
．
第12回演習内容




REPORT1222_1，REPORT1222_2 解答
内部サブルーチン・内部関数の使い方
CONTAINS文（教科書p.132）
二分探索法・Newton法を用いた方程式の求解
台形則を用いた数値積分
REPORT1222_1の解答例
関数副プログラムFACTの掲載は省略（第11回演習資料の解答例参照）
PROGRAM REPORT1222_1
IMPLICIT NONE
INTEGER :: I, N
INTERFACE
FUNCTION NCR(N, R) ! 整数型関数NCRの仮引数の定義
IMPLICIT NONE
INTEGER :: NCR
INTEGER, INTENT(IN) :: N, R
END FUNCTION NCR
END INTERFACE
! 二項係数nCrの計算
FUNCTION NCR(N, R)
IMPLICIT NONE
INTEGER :: NCR
INTEGER, INTENT(IN) :: N, R
INTERFACE
FUNCTION FACT(N)
IMPLICIT NONE
INTEGER :: FACT
INTEGER, INTENT(IN) :: N
END FUNCTION FACT
END INTERFACE
READ(*,*) N
! 整数Nの入力
DO I = 0, N
! パスカルの三角形を出力
! 二項係数の計算
WRITE(*,*) ‘N = ‘,I, (NCR(I,J), J=0, I) !二項係数iC0～iCiを出力 NCR = FACT(N)/FACT(R)/FACT(N-R)
END DO
STOP
END PROGRAM REPORT1222_1
RETURN
END FUNCTION NCR
REPORT1222_2の解答例
関数副プログラムFACTは省略
PROGRAM REPORT1222_2
WRITE(*,*) ‘Number of terms = ‘, N
IMPLICIT NONE
WRITE(*,*) ‘ Summation S = ‘, S
REAL, PARAMETER :: EPS = 1.0E-5 ! 級数の打ち切り誤差
WRITE(*,*) ‘Error = ‘, COS(X) - S
REAL :: A, X, S
STOP
INTEGER :: N
END PROGRAM REPORT1222_2
INTERFACE
FUNCTION FACT(N) ! 整数型関数FACTの仮引数の定義
IMPLICIT NONE
• 関数ＦＡＣＴを用いると級数展開の一般項の式が
INTEGER :: FACT ! N! の計算値
そのまま記述できる
INTEGER, INTENT(IN) :: N
END FUNCTION FACT
• 数式の記述の間違い × -1.0**REAL(I)
END INTERFACE
• 実数計算式中の定数，変数，関数の型の不一致
READ(*,*) X
! 独立変数ｘの入力、単位rad
A = 1.0 ! 級数の一般項変数Aに初項a0 =1 の値を初期設定
S=A
! 級数の合計値変数Sに初項の値を初期設定
N=0
! 級数の項番Nに初項の０を初期設定
DO WHILE (ABS(A)> EPS) ! 一般項の大きさがEPS以下になるまで繰り返す
N=N+1
A = (-1.0)**REAL(I)*X**REAL(2*I)/REAL(FACT(2*I))
漸化式で求めていた級数の一般項
S=S+A
が階乗関数を用いて代入式で記述可能となる．
END DO
Fortran90/95のプログラム単位の分類


メインプログラム
サブルーチン

外部サブルーチン（広義のサブルーチン）学習済
他のプログラムと独立して定義されたサブルーチン副プログラム

内部サブルーチン ← 本日の演習で学ぶ
プログラムの内部で定義されたサブルーチン副プログラム

関数

外部関数（広義の関数）学習済
他のプログラムと独立して定義された関数副プログラム

内部関数 ← 本日の演習で学ぶ
プログラムの内部で定義された関数副プログラム
内部サブルーチン・内部関数とは
教科書p.121




独立したプログラム単位（メインプログラム，外部サブルーチン，外部関数）
の中でのみ通用するサブルーチン・関数をそのプログラム単位の中で記述
したもの．
内部サブルーチン・内部関数での処理に用いられるローカル変数は，全て
上位のプログラム単位のローカル変数として管理される．
副プログラムで使用する変数は，呼び出し側のプログラムの変数として共
有されるグローバル変数となる．仮引数や関数名を除く全ての変数は上位
プログラムにおいて型宣言を行う必要がある．
副プログラムとその上位のプログラムの変数管理が独立していないため，
勝手に変数の内容が副プログラムで書き替えられてしまう危険性がある
→ 外部サブルーチン・外部関数の使用が推奨される．しかし，内部関数は、プログラム中で
関数をINTERFACE文無しで手軽に定義して使えるので便利（旧規格FORTRANの文関数に
相当する）
内部サブルーチン・内部関数の記述方法
メインプログラムではSTOP～END PROGRAMの間、
副プログラムでは、RETURN～END SUBROUTINE/FUNCTIONの間
にそれぞれCONTAINSに続いてサブルーチン/関数を記述する
内部関数の記述方法
内部サブルーチンの記述方法
PROGRAM ○○○
PROGRAM ○○○
IMPLICIT NONE
変数の型宣言文
IMPLICIT NONE
変数の型宣言文
INTEGER :: ○○○
内部関数で使う変数もここで宣言
INTEGER :: ○○○
内部サブルーチンで使う変数もここで宣言
REAL :: ○○○
型宣言文の後のINTERFACE文
REAL :: ○○○
型宣言文の後のINTERFACE文
・・・・・
は書かない
・・・・・
は書かない
・・・・・
・・・・・
STOP
STOP
CONTAINS
これより下に内部サブルーチン・関数を記述
SUBROUTINE △△△(○, ・・・)
IMPLICIT NONE
INTEGER, INTENT(○○) :: ○ 仮引数以外の型宣言は
REAL, INTENT(○○) :: ・・・
書くことができない
具体的なサブルーチンでの処理内容の記述
・・・・・
RETURN
END SUBROUTINE △△△
END PROGRAM ○○○
CONTAINS これより下に内部サブルーチン・関数を記述
FUNCTION △△△(○, ・・・)
IMPLICIT NONE
REAL :: △△△
関数名と仮引数の型宣言
INTEGER, INTENT(IN) :: ○○ は、外部関数と同じ
REAL, INTENT(IN) :: ・・・
関数名と仮引数以外
具体的な関数の計算内容の記述
の型宣言は書けない
・・・・・
△△△ = ・・・
RETURN
END FUNCTION △△△
END PROGRAM ○○○
内部関数のサンプルプログラム
PROGRAM ASF
IMPLICIT NONE
REAL :: A, B, X, Y
INTEGER :: I
READ(*,*) A, B
DO I = 1, 5
READ(*,*) X
Y = F(X)
WRITE(*,*) ‘ Y =‘, Y
END DO
STOP
CONTAINS
FUNCTION F(X)
IMPLICIT NONE
REAL :: F
REAL, INTENT(IN) :: X
F = A*X + B
RETURN
END FUNCTION F
END PROGRAM ASF
変数A,Bは，プログラム単位ASF内の変数で
あるため呼び出し側のメインプログラムと
内部関数副プログラムから変数の参照と代入が可能．
このため，意図しない変数の値の変更が行われ
てしまう可能性がある．これを避けるには，
他の変数にその値をコピーして計算に用いる工夫が必要．
! 変数Xは、仮引数としてメインプログラムより渡される変数A, Bはメインプログラムで
! 定義された変数A, BをメインプログラムASFのグローバル変数として、参照している
! もちろんA, Bは引数として渡すことも可能であるが，この方がシンプルである．
内部サブルーチン・内部関数の
記述上の注意






呼び出し側のプログラムにおけるINTERFACE～END INTERFACE文
の記述は不要
内部サブルーチンでは仮引数のみ型宣言を行う。
内部関数では、関数名と仮引数のみ型宣言を行う。
副プログラム内部での処理で用いる仮引数・関数名以外のすべての変
数は、呼び出し側のプログラムにおいて型宣言を行う。
副プログラムの仮引数以外のメインプログラムの変数も参照可能（グロ
ーバル変数）（注：外部サブルーチン、外部関数では、引数による受け
渡し以外、他のプログラム単位の変数の値の参照は不可能）
メインプログラムのすべての変数が副プログラム側で参照・代入可能な
ため、副プログラム側でのグローバル変数への代入については十分注
意する。必要に応じて他の変数に値をコピーして使うこと。
課題１ ENSYU0112_1




教科書p.134の例題6.5 二分探索法による超越方程式の根
f (x) = 0 を求めるプログラムを入力して実行しなさい．
内部関数を用いた方程式の定義の記述方法を確認せよ
実行結果を確認したら，EXIT文或いはDO WHILE文を用い
て，許容誤差EPS=10-5以下になるまで，二分探索を繰り返
す処理に書き直せ．(条件付繰り返し処理は演習済み）
二分探索法のアルゴリズムは，第４回演習で取り扱った.
忘れた者は、次ページからの解説或いは第４回演習のスラ
イドを読むこと．
２分探索法（Binary Search)の
アルゴリズムの解説１
中間値の定理により，方程式f(x)=0の解の存在範囲が分か
ると，解の存在区間を絞り込み，はさみうち法により解を求め
ることができる．
この時，２分探索法を用いると解の存在範囲の絞り込みが
速いので，方程式の収束解が素早く得られる．
２分探索法は，解の存在範囲を２分割し,どちら側に解が存
在するかを調べて，解の存在範囲の絞り込みを行なう．
中間値の定理（微分積分学Ⅰで習った）
区間[a,b]において関数f(ｘ)は連続とする．f(a)とf(b)の間の実数kに対して，
f(c) = k となる cが a < c < bの範囲に少なくとも１つ存在する．
２分探索法（Binary Search)の
アルゴリズムの解説２
f(x)= x – cos(x)は，全ての実数に対して連続．
f(0)= -1<0, f(/2)=/2 > 0なので，
中間値の定理より，x=0とx=/2の間にf(x)=0の
解が少なくとも１つ存在する．
そこで，解が存在する範囲を絞り込むために，
１）区間(0, /2)を２分割し，どちらの区間に解が
存在するかを調べる．
２） f(/4) > 0ならば，左側の区間(, /4) ，
f(/4) < 0ならば，右側の区間(/4, /2)
に解は存在する．
３） f(/4) > 0なので，新しい区間(, /4) に対して，
同じ手順で範囲の絞り込みを行う．
絞り込みを繰り返すと区間の幅は，(1/2)nで小さ
くなる．その幅が必要な精度以下になったら繰り返し
計算を止める．
f(x) = x - cos(x)
f(XR) > 0
解の存在範囲１
1
解の存在範囲２
f(XM) > 0
XL=0
0
0
0.5
1
XM
-1 f(XL) < 0
1.5
XR
=/2
= (XL+XR)/2
= /4
２分探索法のプログラミング
解の存在区間を(XL,XR)とする．ただし，XLとXRの初期値は，
F(XL)・F(XR) < 0を満足するように与えられるとする．
|XL-XR| > の間処理を繰り返す
DO WHILE文
f(x) = x - cos(x)
f(XR) > 0
解の存在範囲１
1
XM = (XL+XR)/2.0
解の存在範囲２
IF THEN ELSE文
0
0
F(XM)> 
No
XL=XM
f(XM) > 0
XL=0
0.5
Yes
XM
XR=XM
-1 f(XL) < 0
END IF文
END DO文
収束解(XL+XR)/2.0と誤差を表示
1
1.5
XR
=/2
= (XL+XR)/2
= /4
f(XR)<0の時，この条件判断式では区間の絞り込みが正しく行われ
ない．f(XL),f(XR)の符号に関わらず正しくプログラムが動くために
は，F(XM)>0の条件式をどのようにしたらよいだろう
は，END DO, END IFに相当する合流点
課題２ REPORT0112_1


三乗根 x  3 a の解を二分探索法で求めるプログ
ラムを作成し，a = 2に対するxの値を求めなさい．
ただし，解の上界をa(=2.0)，下界は1.0とする．
ただし、課題１で作成した二分探索法のプログラム
を用いて、方程式f(x)=0のf(x)を定義する内部関数
副プログラムのみ書き換えること．
注意）当然のことであるが，1/3乗を使わずに解くこと．
f(X) = 0 の方程式をどのように定義すれば良いか考
えよ．2の三乗根は，1.25992である．
課題３台形公式を用いた定積分
の計算 REPORT0112_2
1

4
台形公式を用いて定積分
dx を計算するプログラムを内部関数を用いて作成
2
0
せよ．
1 x
この定積分の厳密解は，４（tan-1(1) - tan-1(0)) = πとなる．
ただし，

積分区間［X0，XN] （実数），分割数N（整数）の値は，キーボードから入力．

被積分関数f(x)は内部関数として定義すること．

積分値とともにATAN関数（教科書P.33）を用いて定積分の厳密解πとの誤差を表示．

積分区間の分割数n が４，１６，６４，２５６の４つの場合について，定積分の誤差がどのよ
5
うに変化したか調べよ．

xn
x0
n 1
h

f ( x)dx   f ( x0 )  2 f ( xi )  f ( xn )
2
i 1

ただし，x1～xn-1は，積分区間 [x0 , xn]をn等分
した分割点のx座標，hは分割幅=(xn – x0)/nである．
分割数n = 10の場合
f(x0)
4
f(x4)
2
f(x) = 4/ (1+x )
ヒント
★台形公式（右図参照）
f(x5)
3
S1 S2
S3 S 4
S5
2
h
f(x10)
S6
S7
S8
S9 S10
1
0
0 x1 0.2 x3 0.4 x5 0.6 x7 0.8 x9
x0
x2
x4
x6
x8
1
x10
台形の面積
S5 = h(f(x4)+f(x5))/2
シンプソンの公式を用いた定積分の計算
教科書p.136 例題６．６
4
0 1  x 2 dx
1
シンプソンの公式を用いて定積分
を計算するプログラムを関数副プログラムを用いて作成せよ．
ただし，上記定積分の解は，４（tan-1(1) - tan-1(0)) = πとなる．プログラム中でATAN関
数（教科書P.33）を用いて計算結果の誤差を表示させること．また，積分区間の分割数nを
４，１６，６４，２５６と変化させたとき，定積分の誤差がどのように変わるか確認せよ．
★シンプソンの公式（右図参照）
x0
m
m 1
h

f ( x)dx   f ( x0 )  4 f ( x2i 1 )  2 f ( x2i )  f ( xn )
3
i 1
i 1

f(x0)
f(x) = 4/ (1+x2)

xn
5
ただし，xi は，積分区間 [x0 , xn]を2m等分
した分割点のx座標，hは分割幅=(xn – x0)/(2m)である．
★関数副プログラムは被積分関数f (x)の計算
に用いる．
(xi-1, yi-1), (xi, yi), (xi+1, yi+1)の3点を通る２次式は，
2
x’= x – xi の座標変換を施して， y  ax  bx  c
ただし，2(ah2+c) = yi-1+ yi+1 , c = yi
h
S i   ( ax  2  bx   c ) dx  
h
分割数2n = 10の場合
4
f(x4)
f(x5)
3
2
S1
f(x10)
S2
h
f(x6)
S3
S4
1
S5
0
0 x1 0.2 x3 0.4 x5 0.6 x7 0.8 x9
x0
x2
x4
x6
x8
1
x10
２次曲線の面積
S3 = h(f(x4)+4f(x5)+f(x6))/3
2
h
h(ah 2  4c)  ( yi 1  4 yi  yi 1 )
3
3
今週の課題〆切１月１８日１９時
1/15(金) 情報基盤センターはセンター試験のため17時で閉館です．
メール送付先：[email protected]
レポート課題（Subject：report0112_1） 3点
課題2のソースプログラムを作成し、課題で指示された入力データに対する実行結
果とともに送付せよ．なお，結果として１）近似解，２）方程式の誤差，３）繰り返し
回数，４）厳密解と近似解との差を出力させること
ソースファイル名：report0112_1.f95、実行結果：result0112_1.txt
レポート課題（Subject：report0112_2） 3点
課題３のソースプログラムを作成し、課題で指示された入力データに対する実行結
果とともに送付せよ．厳密解との誤差が小さいことを確認すること。
ソースファイル名：report0112_2.f95、実行結果：result0112_2.txt
注意


プログラム名、ソースファイル名、実行結果ファイル名は，各課題での指示通りとすること。
字下げ、コメントがないプログラムは減点する．

Download Report