プログラミング演習２（課題３，第３回）

ベクトル，行列モジュールを用いて連立方程式を解く

練習問題１）以下のサンプルコードを参考にして，LU分解を行う関数(LU_decomp())と， 連立方程式を解く関数(LU_solve())を作成してください．
但し，関数のプロトタイプは変更してもかまいません．

#ifndef LU_H
#define LU_H

#include "array1d.h"
#include "array2d.h"

/*******************************************************************************
 * n×n行列AをLU分解する
 * (LU分解後に行列Aは左下三角行列Lと右上三角行列Uによって上書きされる)
 * @param A LU分解対象の正方行列
 * @param n 正方行列Aの行、列数
 * @return 行列式
 ******************************************************************************/
double
LU_decomp(matrix A, int n);

/*******************************************************************************
 * 元のn×n行列をLU_decompを用いてLU分解した後に、
 * L行列とU行列によって上書きされた行列Aを用いて、
 * n×1行列xとbとするときに、Ax=bを解く
 * (右辺ベクトルとして与えたbは，解ベクトルxによって上書きされる)
 * @param A 元の行列をLU分解後、その左下三角成分をL、その右上三角成分をUとする行列
 * @param b 方程式の右辺ベクトル（解ベクトルによって上書きされる）
 * @param n 行列Aの行、列数、ベクトルbの要素数
 ******************************************************************************/
void
LU_solve(matrix A, vector b, int n);

#endif //LU_H

#include "lu.h"

double
LU_decomp(matrix A, int n) {
    /* LU分解関数を実装せよ */
}

void
LU_solve(matrix A, vector b, int n) {
    /* LU_decompによってLU分解されたL、U行列を用いて連立方程式を解くプログラムを実装せよ */
}

#include <stdio.h>
#include <stdlib.h>
#include "array1d.h"
#include "array2d.h"
#include "lu.h"

int
main(void) {
    int N; /* 行列のサイズ */
    double det; /* 行列式 */
    matrix A; /* 元の行列 */
    vector b; /* 右辺ベクトル */

    N = 4; /* 行列サイズの指定 */

    /* 元の行列の初期化 */
    A = ARRAY2D_new(N, N);
    A[0][0] = 8;
    A[0][1] = 16;
    A[0][2] = 24;
    A[0][3] = 32;
    A[1][0] = 2;
    A[1][1] = 7;
    A[1][2] = 12;
    A[1][3] = 17;
    A[2][0] = 6;
    A[2][1] = 17;
    A[2][2] = 32;
    A[2][3] = 59;
    A[3][0] = 7;
    A[3][1] = 22;
    A[3][2] = 46;
    A[3][3] = 105;

    /* 右辺ベクトルの初期化 */
    b = ARRAY1D_new(N);
    b[0] = 160;
    b[1] = 70;
    b[2] = 198;
    b[3] = 291;

    /* ベクトル・行列の表示*/
    /* 【解答例】 元の行列A =
     * [8.000000,16.000000,24.000000,32.000000,
     * 2.000000,7.000000,12.000000,17.000000,
     * 6.000000,17.000000,32.000000,59.000000,
     * 7.000000,22.000000,46.000000,105.000000]
     * 右辺ベクトルb =
     * [160.000000,70.000000,198.000000,291.000000] */
    printf("元の行列A = \n");
    ARRAY2D_print(A, N, N, "%f");
    printf("右辺ベクトルb = \n");
    ARRAY1D_print(b, N, "%f");

    printf("\nAx=bを解く\n");

    /* LU分解実行 */
    /* 【解答例】 行列式 = 768.000000 */
    det = LU_decomp(A, N);
    printf("行列式 = %f\n", det);

    /* LU分解後の行列Lと行列Uの表示 */
    /* 【解答例】 LU分解後の行列
     * [8.000000,2.000000,3.000000,4.000000,
     * 2.000000,3.000000,2.000000,3.000000,
     * 6.000000,5.000000,4.000000,5.000000,
     * 7.000000,8.000000,9.000000,8.000000] */
    printf("LU分解後の行列\n");
    ARRAY2D_print(A, N, N, "%f");

    /* LU分解後の行列を用いて連立方程式を解く */
    LU_solve(A, b, N);

    /* 解ベクトル表示 */
    /* 【解答例】 解ベクトル
     * [4.000000,3.000000,2.000000,1.000000] */
    printf("解ベクトル\n");
    ARRAY1D_print(b, N, "%f");

    return (EXIT_SUCCESS);
}

LU分解ルーチンの別のプロトタイプ例

#ifndef LU_H
#define LU_H

#include "array1d.h"
#include "array2d.h"

/*******************************************************************************
 * n×n行列AをLU分解する
 * @param A LU分解対象の正方行列
 * @param L LU分解対象の左下三角行列
 * @param U LU分解対象の右上三角行列
 * @param n 正方行列Aの行、列数
 * @return 行列式
 ******************************************************************************/
double
LU_decomp(const matrix A, matrix L, matrix U, int n);

/*******************************************************************************
 * 元のn×n行列をLU_decompを用いてLU分解した後に、
 * L行列とU行列を用いて
 * n×1行列xとbとするときに、LUx=bを解く
 * @param L LU分解後の左下三角行列
 * @param U LU分解後の右上三角行列
 * @param b 方程式の右辺ベクトル
 * @param x 方程式の解ベクトル
 * @param n 行列Aの行、列数、ベクトルbの要素数
 ******************************************************************************/
void
LU_solve(const matrix L, const matrix U, const vector b, vector x, int n);

#endif //LU_H