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

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

| 3.1  1.3 -5.7 |   | x1 |   | -1.3 |
| 1.0 -6.9  5.8 | * | x2 | = | -0.1 |
| 3.4  7.2 -8.8 |   | x3 |   |  1.8 |

をLAPACKを用いて解き x1, x2, x3 を求めよ。

回答例：

[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)
でサブルーチンを呼び出す。
> f77 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 |   | x1 |    | x1 |
| 1.0      -6.9+3.2i  5.8+2.2i | * | x2 | = k| x2 |
| 3.4-4.0i  7.2+2.9i -8.8+3.2i |   | x3 |    | x3 |

を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)
と複素行列を指定する。
> f77 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)はファイルで与える。[ AT（注：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個で構成された力学系の縦波成分を考える。

    このような線形バネ系のハミルトニアンは

    で与えられる。また、線形の系なので解として

    の調和振動の形を仮定することができ、固有振動数（ω）と固有振動モード（その振動数での各質点の振幅：α=[α1,α2,...,αN-1]）は以下の固有値方程式で表すことが出来る（注２）。

    N=5、およびN=6の場合について、すべての固有振動数（ω）と固有振動モード（α=[α1,α2,...,αN-1]）を求めよ。


    注１：本学総合人間学部基礎科学科冨田博之先生の講義ノート[PDF]より拝借させていただきました。

    注２：線形バネ系の振動解析については適当な力学の参考書を見よ。上記の冨田博之先生の講義ノートにも、大変詳しく分かりやすい解説がある。

    ヒント：N=4の場合、以下の対称実３重対角行列の固有値問題に帰着する。[Fortran]

      |  2 -1  0 |    | α1 |         | α1 |
      | -1  2 -1 | * | α2 | = λ/K| α2 |
      |  0 -1  2 |    | α3 |         | α3 |

    固有値ω、固有振動モードαともにそれぞれ３個求まる。各モードでのi番目の質点の振幅[α1,α2,...,αN-1]をiに対してプロットすると、振動数の低いモードから順に１倍波、２倍波、…となっていることがわかる。
    　

《レポート》

  1. 今回の演習１〜３で用いたプログラムと、１、２の結果をメールで以下に送ること。

    提出先：　ryoichi@r02.mbox.media.kyoto-u.ac.jp
    締切り：　１１月３０日（火）

    氏名、学籍番号を忘れずに。授業に関する要望・感想があれば遠慮せずに書いてください。

    提出者：　

    0500-16-3529

    0500-15-0173
    0500-15-0389
    0500-15-0404
    0500-15-0638
    0500-15-0647
    0500-15-0656
    0500-15-0807
    0500-15-0852
    0500-15-1279
    0500-15-1528
    0500-15-1537
    0500-15-1804
    0500-15-1921
    0500-15-1958
    0500-15-1976
    0500-15-2033
    0500-15-2991
    0500-15-3077
    0500-15-3246
    0500-15-3596
    0500-15-3676
    0500-15-3845
    0500-15-3999
    0500-15-4191
    0500-15-4271
    0500-15-4252
    0500-15-4539
    0500-15-4557
    0500-15-5142
    0500-15-6211
    0500-15-6220
    0500-15-6337
    0500-15-6892
    0500-15-6991
    0500-15-7236
    0500-15-7281
    0500-15-7451
    0500-15-7782
    0500-15-8091
    0500-15-8233
    0500-15-8289
    0500-15-8313
    0500-15-8725
    0500-15-8805
    0500-15-8823
    0500-15-9179
    0500-15-9419
    0500-15-9464
    0500-15-9473
    0500-15-9571
    0500-15-9581
    0500-15-9894

    0500-14-0721
    0500-14-1030
    0500-14-1245
    0500-14-1746
    0500-14-2350
    0500-14-2476
    0500-14-2529
    0500-14-2645
    0500-14-2734
    0500-14-2841

    0500-13-0420
    0500-13-0466

    0500-11-2154
    0500-11-3071

    6100-15-2893
    6100-14-1218

    ＊受け取りの確認の意味で学生番号を載せています。掲載を希望しない人があれば対応しますので、レポート毎にメールに希望を記載してください。