4−５．行列演算のまとめと演習

例題１．連立方程式 [SGESV]

| 3.1  1.3 -5.7 |   | x₁ |   | -1.3 |
| 1.0 -6.9  5.8 | * | x₂ | = | -0.1 |
| 3.4  7.2 -8.8 |   | x₃ |   |  1.8 |

をLAPACKを用いて解き x₁, x₂, x₃ を求めよ。

回答例：

[Fortran]
A(1,1)=3.1
A(1,2)=1.3
A(1,3)=-5.7
A(2,1)=1.0
A(2,2)=-6.9
A(2,3)=5.8
A(3,1)=3.4
A(3,2)=7.2
A(3,3)=-8.8
と行列を指定し
call SGESV(3, 1, A, 3, pivot, b, 3, ok)
でサブルーチンを呼び出す。
> gfortran linear.f -o linear -L. -llapack -lblas
> ./linear
1.00000012
1.
1.

[C]
AT[1]=3.1;
AT[2]=1.0;
AT[3]=3.4;
AT[4]=1.3;
AT[5]=-6.9;
AT[6]=7.2;
AT[7]=-5.7;
AT[8]=5.8;
AT[9]=-8.8;
と列連続な１次元配列で行列を指定し、
sgesv_(&c1, &c2, AT, &c1, pivot, b, &c1, &ok);
でサブルーチンを呼び出す。サブルーチン名の後の"_"を忘れずに。
> gcc linear.c -o linear -L. -llapack -lblas -lm -lg2c
> ./linear
1.000000e+00
1.000000e+00
1.000000e+00

　

例題２．固有値問題 [ZGEEV]

| 3.1-1.8i  1.3+0.2i -5.7-4.3i |   | x₁ |    | x₁ |
| 1.0      -6.9+3.2i  5.8+2.2i | * | x₂ | = k| x₂ |
| 3.4-4.0i  7.2+2.9i -8.8+3.2i |   | x₃ |    | x₃ |

をLAPACKを用いて解き、すべての固有値kを求めよ。

回答例：

[Fortran]
complex*16 A(3,3)
A(1,1)=(3.1, -1.8)
A(1,2)=(1.3, 0.2)
A(1,3)=(-5.7, -4.3)
A(2,1)=(1.0, 0)
A(2,2)=(-6.9, 3.2)
A(2,3)=(5.8, 2.2)
A(3,1)=(3.4, -4.0)
A(3,2)=(7.2, 2.9)
A(3,3)=(-8.8, 3.2)
と複素行列を指定する。
> gfortran eigen.f -o eigen -L. -llapack -lblas
> ./eigen
(1.47067795,-3.19468549)
(-0.843010085,7.51121797)
(-13.2276682,0.283467657)

[C]
double AT[2*3*3];
AT[1]=3.1;  AT[2]=-1.8;
AT[3]=1.0;  AT[4]=0.0;
AT[5]=3.4;  AT[6]=-4.0;
AT[7]=1.3;  AT[8]=0.2;
AT[9]=-6.9; AT[10]=3.2;
AT[11]=7.2; AT[12]=2.9;
AT[13]=-5.7;AT[14]=-4.3;
AT[15]=5.8; AT[16]=2.2;
AT[17]=-8.8;AT[18]=3.2;
とまず実虚が連続、次に列が連続な１次元配列で複素行列を指定する。
> gcc eigen.c -o eigen -L. -llapack -lblas -lm -lg2c
> ./eigen
1.470678 -3.194685
-0.843010 7.511218
-13.227668 0.283468
　

《演習》

1. 次の連立１次方程式の解をLAPACKを使って求めよ。ただし実数一般行列A(100x100)、実数ベクトルb(100)はファイルで与え る。[ A^T（注：Cで連続して読み込めるように 行と列を転置させた形でデータが入っています 。）, b ]　
   
   A * x = b
   
   ヒント１：正しくプログラミング出来ていれば、すべての解の値が誤差内で同じ値になる
   
   ヒント２：ファイルの読み込み方法の具体例
   
   [Fortran]
   Real*4 A(100,100), b(100)
   （中略）
   open(10,file="a.dat")
   open(11,file="b.dat")
   
   do j=1,100
   read(10,*)(a(i,j),i=1,100)
   enddo
   read(11,*)(b(i),i=1,100)
   
   close(10)
   close(11)
   
   ただし今回のようにデータが列連続で与えられている場合のみ以下でも可！！
   read(10,*)a
   read(11,*)b
   （中略）
   上でAT, bを単精度(Real or Real*4)で定義しているので、ここではSGESVを使う。もしAT, b を倍精度(Real*8 or Double)で定義するならDGESVを使えばよい。
   
   [C]
   #include <stdio.h>
   #define N 100
   （中略）
   double AT[N*N];
   double b[N];
   FILE *fpa,*fpb;
   （中略）
   fpa=fopen("a.dat","r");
   fpb=fopen("b.dat","r");
   for(i=0;i<N*N;i++){
   fscanf(fpa,"%lf",&AT[i]);
   }
   for(i=0;i<N;i++){
   fscanf(fpb,"%lf",&b[i]);
   }
   fclose(fpa);
   fclose(fpb);
   （中略）
   上でAT, bを倍精度(double)で定義しているので、ここではdgesv を使う。もしAT, b を単精度(float)で定義するならsgesvを使えばよい。
   　
2. 上の行列Aの固有値をすべて求めよ。（うまくいかない人は行列Ａ、ベクトルbともに複素数で定義 してやってください。その場合、ＺＧＥＥＶを用います。）
   
   
3. 下図（注１）のような定数Kの線形バネN個と質量mの質点N-1個で構成された力学系の縦波成分を考える。
   
   このような線形バネ系のハミルトニアンは
   
   で与えられる。また、線形の系なので解として
   
   の調和振動の形を仮定することができ、固有振動数（ω）と固有振動モード（その振動数での各質点の振幅：α=[α₁,α₂,...,α_N-1]） は以下の固有値方程式で表すことが出来る（注２）。
   
   N=5、およびN=6の場合について、すべての固有振動数（ω）と固有振動モード（α=[α₁,α₂,...,α_N-1]） を求めよ。
   
   
   注１：本学総合人間学部基礎科学科冨田博之先生の講義ノート[PDF] より拝借させていただきました。
   
   注２：線形バネ系の振動解析については適当な力学の参考書を見よ。上記の冨田博之先生の講義ノートにも、大変詳しく分かりやすい解説がある。
   
   ヒント：N=4の場合、以下の対称実３重対角行列の固有値問題に帰着する。[Fortran]
   
     |  2 -1  0 |    | α₁ |         | α₁ |
     | -1  2 -1 | * | α₂ | = λ/K| α₂ |
     |  0 -1  2 |    | α₃ |         | α₃ |
   
   固有値ω、固有振動モードαともにそれぞれ３個求まる。各モードでのi番目の質点の振幅[α₁,α₂,...,α_N-1] をiに対してプロットすると、振動数の低いモードから順に１倍波、２倍波、…となっていることがわかる。