Lec 3 | MIT 18.02 Multivariable Calculus, Fall 2007

　多変数微積分学の第３回目の講義を見た。

外積の応用と行列の基礎の話だった。だいたい知ってる内容だった。

3次元空間に3つの点P1, P2, P3がある。これらの点は平面上にある。

もうひとつの点Pが与えられる。Pは上の3点と同一平面上にあるか？

という問題が面白かった。

• determinantを使ってparallelepipedの体積が0になるかどうかで判定する方法
• cross productを使ってnormal vectorを出した後、それが平面上のvectorと直交するかどうかで判定する方法

の2つが紹介されていた。

Copyright Information

Prof. Denis Auroux, MIT 18.02 Multivariable Calculus, Fall 2007

View the complete course at: http://ocw.mit.edu/18-02F07

License: Creative Commons BY-NC-SA 

また、講義の終盤で逆行列を求める方法（小さいサイズの行列に対してのみ有効）が紹介されていた。プログラムを作って動かしてみた。4*4の行列まで動作確認済。

#include <cstdio>
#include <vector>

using namespace std;

typedef vector<double> vec;
typedef vector<vec> mat;

/*
 * remove row r and colum c from a matrix
 */
mat removeRC(const mat &A, int r, int c) {
    int n = A.size();

    mat ret;

    for (int i = 0; i < n; i++) {
        if (i == r) 
            continue;

        vec row;
        for (int j = 0; j < n; j++) {
            if (j == c)
                continue;
            row.push_back(A[i][j]);
        }

        ret.push_back(row);
    }

    return ret;
}

/*
 * calculate determinant of a matrix.
 */
double det(const mat &A) {

    int n = A.size();

    if (n == 1)
        return A[0][0];

    double ret = 0;

    for (int i = 0; i < n; i++) {
        double tmp = A[0][i] * det(removeRC(A, 0, i));
        if (i % 2 == 0)
            ret += tmp;
        else
            ret -= tmp;
    }

    return ret;
}

/*
 * make adjoint of a matrix.
 */
mat adj(const mat &A) {
    
    int n = A.size();
    mat ret(n, vec(n));

    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            ret[j][i] = det(removeRC(A, i, j));
            if ((i + j) % 2)
                ret[j][i] = -ret[j][i];
        }
    }

    return ret;
}

/*
 * make an inverse of a matrix.
 */
mat inv(const mat &A) {

    int n = A.size();
    double d = det(A);
    mat ret = adj(A);
    
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            ret[i][j] /= d;

    return ret;
}

/*
 * test use only
 */
void disp(const mat &A) {
    int n = A.size();
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++)
            printf("%.8lf ", A[i][j]);
        puts("");
    }
    puts("");
}

/*
 * test use only
 */
mat multiply(const mat &A, const mat &B) {
    int n = A.size();
    mat ret(n, vec(n));
    
    for (int i = 0; i < n; i++)
        for (int k = 0; k < n; k++)
            for (int j = 0; j < n; j++)
                ret[i][j] += A[i][k] * B[k][j];
    
    return ret;
}

int main(int argc, char **argv) {

    // 2*2 matrix    
    mat A1 = {
        {2, 1},
        {5, 3}
    };
    disp(multiply(A1, inv(A1)));

    // 3*3 matrix
    mat A2 = {
        {1, 2, 5},
        {1, -1, 1},
        {0, 1, 2}
    };
    disp(multiply(A2, inv(A2)));

    // 4*4 matrix
    mat A3 = {
        {0, -1, -1, 0},
        {3, 0, 1, 1},
        {-3, 0, 1, -1},
        {0, 1, 2, 3},
    };
    disp(multiply(A3, inv(A3)));
    
    return 0;
}