４日目「行列・ベクトル」：準備
能書き
コンピュータは連続な変数を扱うことができないので、 物理の問題はしばしば離散化してから扱うことになります。このとき、行列やベクトル（線形代数）は強力なツールとなるため、 コンピュータ上で自在に操れなければなりません。
今のうちに練習して、感覚をつかんでおきましょう。
ファイルを作る
program day4 !------------------------ ! Fortran template !------------------------ implicit none integer, parameter :: double=selected_real_kind(15) !--- begin header --- !--- end header --- !--- begin main --- !--- end main --- end program
Lesson 11. 配列
定義：ベクトル・行列の宣言
n成分のベクトルは『変数名(n)』、m x nの行列は『変数名(m,n)』として宣言する。型 :: ベクトル変数名(indexの上限) 型 :: 行列変数名(index1の上限,index2の上限)例）3成分の倍精度実数型ベクトル『avec』、2 x 2の倍精度実数型行列『Amat』を宣言
real(kind=double) :: avec(3), Amat(2,2)
変数名の大文字小文字は区別しません。なので、行列『A』とスカラー変数『a』を同時に宣言することはできません。
定義：ベクトル・行列を使う
ベクトル・行列は「配列」という変数の集合であって、普通の変数とは扱いが異なる。ベクトル・行列は、『配列名(i,j)』のように要素のindexを指定することで、 普通の変数と同じように使うことが出来る。例）行列『Amat』の (1行, 2列) 成分に『0』を代入
Amat(1,2) = 0
たとえ中身が１成分しか無くても、「値」を持てるのは「配列『a』」ではなく、「配列の１つめの成分『a(1)』」です。
学習
普通の変数との違いに注意しながら下のプログラムをよく読んで、 実行してみましょう。 integer :: i, j
real(kind=double) :: rvec(3), Amat(2,2)
!--- vector r ---
do i=1,3
rvec(i) = i
end do
print *, "x = ", rvec(1)
print *, "y = ", rvec(2)
print *, "z = ", rvec(3)
!--- matrix A ---
do i=1,2
do j=1,2
Amat(i,j) = i+2.0d0*j
end do
end do
do i=1,2
do j=1,2
print *, "A(",i,",",j,") =", Amat(i,j)
end do
end do
x = 1.00000000000000 y = 2.00000000000000 z = 3.00000000000000 A( 1 , 1 ) = 3.00000000000000 A( 1 , 2 ) = 5.00000000000000 A( 2 , 1 ) = 4.00000000000000 A( 2 , 2 ) = 6.00000000000000『do』文が２つ、続けて出てくる部分がありますね。 これも確実に理解してほしいポイントなので、次で詳しく説明します。
Lesson 12. 多重ループ
学習： 多重ループ
ループは多重構造にすることができる。 このとき、内側のループが、外側のループにより繰り返し実行される。ループの内側にループを作ることを「多重ループ」と言います。
このあたりからプログラムの構造が複雑になってきますが、
行列を扱うには必要不可欠なので頑張ってマスターしましょう。
例えば
do i=1,n
do j=1,m
{処理}
end do
end do
だと、『n』x『m』回、{処理}を繰り返します。内側のループ変数、この例の場合だと『j』が先に変化することに注意して下さい。
すなわち、(i,j)=(1,1),(1,2),...,(1,m),(2,1),(2,2),...,(n,m)という順番になります。
「iとjの両方が変わるんだな」と、なんとなく理解するだけではプログラミングは上達しません。 内側と外側の実行の順番を意識して、変数の値の変化を正しく追えるようにしましょう。
実践
単位行列を定義するプログラムを作ってみましょう。 integer :: i, j
real(kind=double) :: Imat(3,3)
do i=1,3
do j=1,3
Imat(i,j) = 0.0d0
end do
end do
do i=1,3
Imat(i,i) = 1.0d0
end do
print *, "Imat"
do i=1,3
print *, Imat(i,1), Imat(i,2), Imat(i,3)
end do
１つめと２つめの（２重）ループで全成分を『0』にしてから、３つめのループで対角成分のみ『1』にしています。Imat 1.00000000000000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 1.00000000000000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 1.00000000000000行列の「成分」を扱う感覚はつかめましたか？
ちなみに、『Imat=1.0d0』は単位行列とは全く別の意味になります。
問題３：行列積を求めるプログラムを作りなさい。
ループを駆使して、行列積を求めてみましょう。 integer :: i, j, k
real(kind=double) :: Amat(3,3), Bmat(3,3), Cmat(3,3)
do i=1,3
do j=1,3
Amat(i,j) = 1.0d0*i
Bmat(i,j) = 3.0d0*i-1.0d0*j
end do
end do
!{ここで C = A B を求める}
print *, "Cmat is"
do i=1,3
print *, Cmat(i,1), Cmat(i,2), Cmat(i,3)
end do
初めての人には結構難問なので、まずは紙の上で手計算してみるといいかもしれません。問題を解く「手順」をしっかり理解しておくことが重要です。
ヒント：Cmatの(1,1)成分は、Amatの(1,1)成分とBmatの(1,1)成分を掛けたもの足す・・・、これをひたすら続ける。
『Cmat=Amat*Bmat』ではダメですよ。念のため。
結果が下のようになれば成功です。
Cmat is 15.0000000000000 12.0000000000000 9.00000000000000 30.0000000000000 24.0000000000000 18.0000000000000 45.0000000000000 36.0000000000000 27.000000000000010分悩んで出来なかったら、一度答えを見てから再挑戦してください。 自分で手順を整理できるよう、理解して練習することが大事です。
>>解答