パソコンで数値計算

正方行列の積

行列の積

2つのn次正方行列A,Bの成分をそれぞれと表記すると、その積は
   (1.1)
と表せます。は積の計算結果となる行列Cの成分です。行列の掛け算は、手計算でもできるほど複雑な計算ではありませんので、わざわざ数値計算させる必要はありませんが、次数が大きい行列の計算等では数値計算の力が役立ちます。

プログラムソース

以下のプログラムでは、掛け算を行う行列が書き込まれた2つのファイル(a.dat,b.dat)を用意する必要があります。fscanf関数は膨大な数値が書き込まれたファイルから数値を読み込むのに大変便利です。ここでは実対称行列の対角化の計算結果で出た直行行列(a.datに保存)
0.577350 0.707107 -0.408248
-0.577350 0.707107 0.408248
0.577350 0.000000 0.816497
とその転置行列(b.datに保存)
0.577350 -0.577350 0.577350
0.707107 0.707107 0.000000
-0.408248 0.408248 0.816497
との積を計算します。計算結果は単位行列になるはずです。

プログラムはエラーチェックを行わない等簡単に書かれていますので、行列の掛け算ができるようにするなど、まだまだ改良の余地はありますので、自分が使いやすいように改良してみると良いでしょう。

#include <stdio.h>
#include <stdlib.h>

#define N 10            //N次の正方行列まで扱えるようにする

void matrixmultiply(int n,double a[N][N],double b[N][N],double c[N][N]);

int main()
{
        int i,j,n;
        double A[N][N],B[N][N],C[N][N];

                FILE *readin1,*readin2;

/*行列の値が書き込まれたファイルを開く*/
                if((readin1=fopen("a.dat","r"))==NULL) {
                        printf("a.datを開けません\n");
                        exit(1);
                }
                if((readin2=fopen("b.dat","r"))==NULL) {
                        printf("b.datを開けません\n");
                        exit(1);
                }

                printf("行列の次数を入力してください\n");
                scanf("%d",&n);
                printf("%d次の正方行列の掛け算を行います\n\n",n);

/*ファイルから数値を読み込み、配列に代入する*/
                for(i=0;i<n;i++) {
                        for(j=0;j<n;j++) {
                                fscanf(readin1,"%lf",&A[i][j]);
                                fscanf(readin2,"%lf",&B[i][j]);
                        }
                }

        matrixmultiply(n,A,B,C);        //関数を呼び出し行列の掛け算を行う。

/*結果を表示する*/
                printf("計算結果\n");
        for(i=0;i<n;i++) {
                for(j=0;j<n;j++) {
                        printf("%f ",C[i][j]);
                }
                printf("\n");
        }

                fclose(readin1);
                fclose(readin2);

        return 0;
}

/*掛け算を行う行列2つと、結果を入れる行列を引数として受け取る。*/
void matrixmultiply(int n,double a[N][N],double b[N][N],double c[N][N])
{
        int i,j,k;

/*受け取った2つの行列の掛け算を行う。*/
        for(i=0;i<n;i++) {
                for(j=0;j<n;j++) {
                        for(k=0;k<n;k++) {
                                c[i][j]+=a[i][k]*b[k][j];
                        }
                }
        }
}

数値計算実行結果

プログラムを実行すると、計算を行う正方行列の次数の入力を求められます。その数値を正しく入力すれば、計算結果が表示されます。結果は単位行列になるはずですが、厳密にはなっていません。これは近似、誤差の影響によるものでしょう。

行列の次数を入力してください
3
3次の正方行列の掛け算を行います

計算結果
1.000000 0.000001 -0.000000
0.000001 1.000000 0.000000
-0.000000 0.000000 1.000000

2009/04/19 更新