FORTRANプログラミング 謔R回の演習 - ax

木村拓馬
FORTRAN プログラミング
–第３回の演習–
木村拓馬
2014 年 10 月 15 日 09:28
FORTRAN プログラミング,–第３回の演習– ( 2014 年 10 月 15 日 09:28 )
1/12
.
演習 3.1（pr3.1.f90）
.
演習 1.2（第一回講義）のプログラム（pr1.2.f90）を，
木村拓馬
n をキーボードから入力（1000 より大きくてもＯＫ）
配列の動的割付け（allocate）を用いる
ように書き換えなさい（自分で作った方は自作のを，作ってない方は
.
.
演習 1.2（pr1.2.f90） .. Download .. Download .. Download .. Download
.
.. Download
を書き換え）
n を 3 以上 1000 以下の整数定数とする（n は READ しなくてよい，PARAMETER は OK）．
7 つの実数 a, b, c, d, e, f, g が与えられたとき，
1. 以下の n × n 三重対角行列 M を２次元配列に記憶させ（空白は零），
 a

 c



M = 




b
d
..
.
e
..
..
.
..
.
.
..
.
c
..
.
d
f






 ,



e 
g

 1
 −1

例えば n = 5 とき M =  0
 0

0
−1
2
−1
0
0
0
−1
2
−1
0
0
0
−1
2
0
0
0
0
−1
1




 ,


. 画面と”pr1.2.dat”なるファイルの両方に出力する，
2
プログラムを作成しなさい．
3. （終わった人は，
「n > 20 ならば画面には出力しない」ように書き換えてみよう）
4. （終わった人は， 1. の行列を作成する部分をサブルーチンにしてみよう）
.
FORTRAN プログラミング,–第３回の演習– ( 2014 年 10 月 15 日 09:28 )
2/12
木村拓馬
.
演習 3.2（pr3.2.f90）
.
ガウスの消去法を用いて n 元連立一次方程式を解くプログラムを作成しよう．
とりあえず，11 ページにある例題を解いてみる．
解けたら演習 3.3 の問題を解いてみる．
係数行列・右辺ベクトルを作る部分では，配列の動的割付けを用いること．
終わった方は，ガウスの消去法の部分を副プログラムにしてみよう．
（演習 3.3 が終わった方は残差の無限大ノルムを計算するプログラムを作成）
線型方程式 Ax = b の解 x に対し，残差は Ax − b．
ふつう数値計算の計算結果には誤差が混入，Ax − b = 0 なる x は得られない．
残差はベクトルで全ての値を表示すると大変なので，残差のノルム ∥Ax − b∥ の値
などを見る．
ベクトルの無限大ノルムは
∥x∥∞ = max {|xi |}
i=1,··· ,n
.
配列の最大値は MAXVAL(x), 絶対値は ABS(x), 行列ベクトル積 Ax は
MATMUL(A,x) を利用できる．
FORTRAN プログラミング,–第３回の演習– ( 2014 年 10 月 15 日 09:28 )
3/12
木村拓馬
.
演習 3.3（pr3.3.f90）
.
1. 整数 n ≥ 3 と正数 T が与えられたとき，以下の n × n 行列 A と n 次元ベクトル b に対
し，連立一次方程式 Ax = b を解くプログラムを作成しよう．

 1
 −1



A = 




−1
2
..
.
−1
..
.
..
.
..
..
.
.
−1
..
.
2
0


− cos(0.0) · h




sin(1 · h) · h2




sin(2 · h) · h2


 , b = 
..




.


 sin((n − 2) · h) · h2
−1 


1
sin((n − 1) · h)






 , h = T π / (n − 1).





n はキーボードから入力，配列の動的割付け（allocate）を用いること
（行列 A を作るのは pr1.2.f90 を参照）．
. T = 2 とおいたときの解 x と，
2
yi := sin(T π (i − 1) / (n − 1)),
i = 1, 2, · · · , n,
で定義される n 次元ベクトル y とを比べてみよう（ x − y，その絶対値，その最大値）
. 終わった人は n = 10, 20, 40 について x と y のグラフを重ねて描いてみよう．
.
3
FORTRAN プログラミング,–第３回の演習– ( 2014 年 10 月 15 日 09:28 )
4/12
木村拓馬
参考文献
[1] JIS X 3001-1:2009 (プログラム言語 Fortran – 第 1 部:基底言語)
[2] 戸川隼人: 「ザ・Fortran90/95」, サイエンス社 (1999)
[3] 伊理正夫: 岩波講座 応用数学 線形代数 I・II, 岩波書店 (1994)
[4] 川久保 勝夫: 線形代数学（新装版）, 日本評論社 (2010)
FORTRAN プログラミング,–第３回の演習– ( 2014 年 10 月 15 日 09:28 )
5/12
木村拓馬
Appndix
ガウス消去法のプログラム作成において，ウェブ検索して見つけた
「特別な場合」を鵜呑みにして採用することがあるとイヤなので，オマケを追加．
以後の内容に合わせて作成する必要は全く無い．
むしろ一年次の線形代数で習った手順を実装してみてほしい．
線形代数や数値解析の教科書に載っている手順でもかまわない．
FORTRAN プログラミング,–第３回の演習– ( 2014 年 10 月 15 日 09:28 )
6/12
木村拓馬
.
線形連立方程式 Ax = b の，拡大係数行列を用いた表現
.
定義（拡大係数行列）
n × n 行列 A と n 次列ベクトル b が与えられたとき，
Ax = b
なる関係が成り立つような n 次列ベクトル x を求めるとする．
 

 

b1 
 x1 
a11 a12 · · · a1n 
b 
 x 
a21 a22 · · · a2n 
2
 2 

 


 . 
,
b
=
,
x
=
A =  .
.
.
.



..
 . 

..
.
 .. 

 ..
.
. 
 . 
 


bn
xn
an1 an2 · · · ann
このとき，A を係数行列，b を右辺ベクトル，x を解ベクトル，係数行列 A と
右辺ベクトル b を横にくっつけた行列を「拡大係数行列」と呼ぶことにするa ．
本講義では，拡大係数行列を (A|b) のように表現する．


 a11 a12 · · · a1n b1 
 a21 a22 · · · a2n b2 


(A|b) =  .
..
..
.. 
..
 ..

.
.
.
.


an1 an2 · · · ann bn
.
a
呼び方は教科書によっていろいろ．右辺，未知数，未知ベクトルなど．
FORTRAN プログラミング,–第３回の演習– ( 2014 年 10 月 15 日 09:28 )
7/12
木村拓馬
.
線形連立方程式 Ax = b と拡大係数行列と行基本変形
.
定義（行基本変形）
行列に対して，
1.
”ある行 ”と ”他の行 ”を入れ替える
2.
”ある行 ”を非零倍する
3.
”ある行の非零倍 ”を他の行に加える
の三つの操作を行基本変形という．
行基本変形と線型連立方程式
Ax = b を拡大係数行列 (A|b) で表現し，(A|b) を行基本変形して (C|d) な
る行列が得られたとする．このとき，
Ax = b
⇐⇒
Cx = d
がいえる．つまり，拡大係数行列を行基本変形すれば，同じ解 x を持つ線
形連立方程式に変形できるa ．
よって，単位行列を I と表すことにして，(A|b) を行基本変形して (I|d) な
る行列を導出すれば，
Ax = b
⇐⇒
Ix = d
⇐⇒
x=d
より，解 x を求めたことになる．
.
a
理由は省略．大抵の線形代数の教科書に書いてあるので，そちらを参照してほしい．例えば [3, 4]
FORTRAN プログラミング,–第３回の演習– ( 2014 年 10 月 15 日 09:28 )
8/12
木村拓馬
.
部分軸選択を伴うガウスの消去法（Gaussian elimination with partial pivoting）
.
Ax = b の拡大係数行列 (A|b) に対して行基本変形を繰り返し，解 x を導出する手法
に ”ガウスの消去法 ”がある．

 a11
 a
 21
(A|b) =  .
 ..

an1
a12
a22
..
.
an2
···
···
..
.
···
a1n
a2n
..
.
ann
b1
b2
..
.
bn







前進消去（forward elimination）
：i = 1, 2, · · · , n − 1 について，
. i 行目以降で i 列目の要素のうち絶対値が最も大きいものを探し，
1
「その要素がある行と 第 i 行を入れ替える」．
.
このとき，A が正則ならば aii , 0 となる．
.2 「 ”第 i 行の非零倍 ”を ”他の行 ”に加えて」，
第 i 列の i + 1 行目以降の要素が全て 0 になるように計算する．
FORTRAN プログラミング,–第３回の演習– ( 2014 年 10 月 15 日 09:28 )
9/12
木村拓馬
.
A が正則ならば，前進消去により (A|b) は以下のような行列 (C|d) に変形される．
このとき，c11 , · · · , cnn は非零となる．

 c11

(C|d) = 

0
···
..
.
c1n
..
.
cnn
d1
..
.
dn





後退代入（backward substitution）a
1. まず第 n 行より， x = d /c
n
n nn が得られる．
2. 次に第 n − 1 行より，
cn−1,n−1 xn−1 + cn−1,n xn = dn−1
である．これを変形して，
xn−1 = (dn−1 − cn−1,n xn )/cn−1,n−1
と書ける． 1. で得た xn を代入し， xn−1 が得られる．
. 同様に，i = n − 2, · · · , 1 について，


3
n
∑


xi = di −
ci j x j  /cii
j=i+1
を計算すれば解 x が得られる．
.
a
行基本変形チックに書いていないが，後退代入の操作は全て行基本変形に対応
FORTRAN プログラミング,–第３回の演習– ( 2014 年 10 月 15 日 09:28 )
10/12
木村拓馬
.
計算例（前進消去）
.

(A|b)
=
 0
 −1

1
−2
1
1
2
−1
1
−4
1
1




⇓ 第 1 行以降，第 1 列で絶対値最大は (2,1) 成分の”-1”
⇓ よって第 1 行と第 2 行を入れ替える


1
−1
1 
 −1

 0
−2
2
−4 


1
1
1
1
⇓ 第

 −1
 0

0
1 列の第 2 行以降を全て零にするには，第 1 行の 1 倍を第 3 行に加える
1
−2
2
−1
2
0
1
−4
2




⇓ 第 2 行以降，第 2 列で絶対値最大は (2,2) 成分の”-2”
⇓ よって行の入れ替えを行う必要はない
⇓ 第

 −1
 0

0
2 列の第 3 行以降を全て零にするには，第 2 行の 1 倍を第 3 行に加える
1
−2
0
−1
2
2
1
−4
−2


 =: (C|d)

注）実際は，
「行のリスト」の１次元配列を作成し，行を入れ替える代わりに，リストを入れ替えて，
参照するように作成する人も多い．
.
FORTRAN プログラミング,–第３回の演習– ( 2014 年 10 月 15 日 09:28 )
11/12
木村拓馬
.
計算例（後退代入）
.
(C|d)
=

 −1

 0
0
1
−2
0
−1
2
2
1
−4
−2




後退代入（backward substitution）
1. まず第 3 行より，次のようにして x が計算できる．
3
x3 = d3 /c3,3 = −2/2 = −1
. 次に第 2 行より，
2
−2 x2 + 2 x3 = −4
である．これを次のように変形して， x2 が計算できる．
x2 = (d2 − c2,3 x3 )/c2,2 = (−4 + (−2) · (−1))/(−2) = 1
. 同様にして， x1 が計算できる．
3
(
)
x1 = d1 − c1,2 x2 − c1,3 x3 /c1,1 = (1 − 1 · 1 − (−1) · (−1)) /(−1) = 1


 1 

よって解 x =  1  が得られる．
.
−1
FORTRAN プログラミング,–第３回の演習– ( 2014 年 10 月 15 日 09:28 )
12/12