Fortran
Cygwin
numpy
初心者
python3

情報処理の問題をfortranとpythonで書く!#01 -2次正方行列の固有値、固有ベクトル

はじめに

大学で出されたfortranの課題をpythonで書いていきます。ここに書いているfortranのコードは自分で書いているため、少し長かったり、あまりよろしくない書き方もあるかもしれません。

※問題文にある規格化ベクトルを求める際、規格化をしていないことに気付いたため、コードを後日、変更します 2018/2/23

以下に、それぞれの環境を記しておきます

fortran90

windows10
cygwin2.8.1

python

windows10
python3.6.2

課題内容 4-13

正則な 2 次正方行列の固有値,規格化された固有ベクトルを求めるプログラムを作れ

それでは、早速それぞれのコードを書いておきます

fortran90

toi4-13.f90
program toi4_13
implicit none
integer::i1
real::r,sol1,sol2
real,dimension(4)::a
real,dimension(2)::b1,b2

a(1)=8
a(2)=1
a(3)=4
a(4)=5

r=(a(1)+a(4))**2-4*(a(1)*a(4)-a(2)*a(3))

if (r<0) then
    write(*,*) 'imaginary solutions!!'
    stop
end if
sol1=(a(1)+a(4)+r**0.5)/2
sol2=(a(1)+a(4)-r**0.5)/2

write(*,'(A,f7.4,A,f7.4)')'eigenvalue=  ', sol1,',',sol2

b1(1)=a(1)-sol1
b1(2)=a(3)

b2(1)=a(2)
b2(2)=a(4)-sol2

write(*,'(A,f5.2,f5.2,A,f5.2,f5.2,A)')  'Eigenvectors=  (', b1,'),(',b2 ,')'

end program toi4_13

出力結果
eigenvalue= 9.0000, 4.0000
Eigenvectors= (-1.00 4.00),( 1.00 1.00)

今回の課題の行列は2次正方行列のため、2次方程式を解いて固有値を求める方針でいきました。

   \begin{vmatrix}
\mathbf{λ-a(1)} & \mathbf{a(2)} \\
\mathbf{a(3)} &  \mathbf{λ-a(4)}  \\\end{vmatrix}=0

より

λ^2-(a(1)+a(4))λ+a(1)*a(4)-a(2)*a(3)=0

少し躓いたところは、最後(固有ベクトル)のwrite(*,'()')の書式指定が配列に関しては、それぞれの要素ごとに指定しないといけないことを知らなくて少し戸惑いました
→f5.3,f5.3,f5.3など連続する場合は3f5.3と略せます

python

4-13.py
#-*-coding:utf-8-*-
import numpy as np

def main():
    matrix=np.array([[8,1]
                   ,[4,5]])
    la,v=np.linalg.eig(matrix)
    print(u"固有値:"+str(la))
    print(u"固有ベクトル:\n"+str(v))

if __name__=='__main__':
    main()

出力結果
固有値:[ 9. 4.]
固有ベクトル:
[[ 0.70710678 -0.24253563]
[ 0.70710678 0.9701425 ]]

pythonに関してはnumpyという拡張モジュールをインポートしておけばnp.linalg.eig関数が使えて、固有値と固有ベクトルがすぐ求められるらしいのでこれだけ短くなりました。さすがオブジェクト指向...
ちなみにnumpyの読みは勝手に「ナンピー」だと思っていたのですが、「ナムパイ」or「ナンパイ」らしいです(wikipedia)

参考にしたサイト
numpyで線形代数