두 고유 벡터는 서로 직교한다. 즉 대칭 행렬의 서로 다른 고유 벡터들은 서로 직교한다. 이와 같이 대칭 행렬에는 그 자신만의 특이한 성질이 있다.
- 실수로 원소로 가지는 대칭 행렬의 고유치들은 모두 실수이다. 따라서 그 고유 벡터들도 실수로 이루어져 있다
- 서로 다른 고유치에 해당하는 고유 벡터들은 서로 직교한다
- 실수로 원소를 가지는 nxn의 대칭 행렬 S가 p중근의 고유치를 가지면 rank값은 n-p가 된다
- 앞의 성질3에 의해서 대칭 행렬 S를 대각화하는 닮음 변환은 항상 존재한다
- 위에서 찾아지는 닮음 변환에 해당하는 행렬 P는 직교 행렬이 된다
대칭 행렬의 제곱근
양의 정부호나 양의 반 정부호인 행렬들의 제곱근은 쉽게 구할 수 있다. 첫 번째 방식은 행렬의 대각화를 이용하는 방식이고, 다른 한 가지 방식은 삼각 행렬 분해법을 이용하는 것이다.
촐레스키법을 이용한 행렬의 제곱근
촐레스키법에 의한 행렬의 삼각 분해라고 할 수 있다.
<수식>
#include <stdio.h>
#include <math.h>
void Choleski_method(float original[20][20], float L[20][20], float U[20][20], int dim) {
int row, column, row2;
float ftemp;
U[0][0] = (float)sqrt((double)original[0][0]);
for (column = 1; column < dim; column++)
U[0][column] = original[0][column] / U[0][0];
for (row = 1; row < dim; row++) {
ftemp = 0;
for (row2 = 0; row2 < row; row2++)
ftemp += U[row2][row] * U[row2][row];
U[row][row] = (float)sqrt((double)(original[row][row] - ftemp));
for (column = row + 1; column < dim; column++) {
ftemp = 0;
for (row2 = 0; row2 < row; row2++)
ftemp += U[row2][row] * U[row2][column];
U[row][column] = (original[row][column] - ftemp) / U[row][row];
}
}
for (row = 0; row < dim; row++)
for (column = 0; column < dim; column++)
L[row][column] = U[column][row];
}
void main() {
int column, row, dim, mul_index;
float original[20][20], L[20][20], U[20][20];
float multiplied[20][20];
printf("entered the dimension fo matrix:");
scanf_s("%d", &dim);
for (row = 0; row < dim; row++)
for (column = 0; column < dim; column++)
L[row][column] = U[row][column] = 0.0;
printf("entered the elements of matrix\n");
for (row = 0; row < dim; row++)
for (column = 0; column < dim; column++) {
printf("A[%d,%d]=", row + 1, column + 1);
scanf_s("%f", &original[row][column]);
}
printf("\n\n*************************************\n");
printf("entered matirx A is: \n");
for (row = 0; row < dim; row++) {
for (column = 0; column < dim; column++) {
printf("%10.3f", original[row][column]);
}
printf("\n");
}
Choleski_method(original, L, U, dim);
for (row = 0; row < dim; row++)
for (column = 0; column < dim; column++) {
multiplied[row][column] = 0;
for (mul_index = 0; mul_index < dim; mul_index++)
multiplied[row][column] += L[row][mul_index] * U[mul_index][column];
}
printf("\n\n****************************************\n");
printf("\n L matirx of A is :\n");
for (row = 0; row < dim; row++) {
for (column = 0; column < dim; column++)
printf("%10.3f", L[row][column]);
printf("\n");
}
printf("\n\n*************************************\n");
printf("U matrix of A is : \n");
for (row = 0; row < dim; row++) {
for (column = 0; column < dim; column++)
printf("%10.3f", U[row][column]);
printf("\n");
}
printf("\n\n**************************************\n");
printf("multiplication of matrix L and matrix u is : \n");
for (row = 0; row < dim; row++) {
for (column = 0; column < dim; column++)
printf("%10.3f", multiplied[row][column]);
printf("\n");
}
}
최소 자승법
주어진 연립 방정식에서 주어진 방정식의 개수가 미지수의 개수보다 많고, 그 방정식들이 서로 독립적이라고 가정한다. 이 경우 해가 존재하지 않고 따라서 가장 근접한 해를 구하는 방법밖에 없다. 실험을 하다 보면 측정 시의 오차 등에 의해서 이러한 경우가 많이 생기는데 이 경우 사용하는 것이 최소 자승법에 의한 해를 구하는 방식이다
가중 최소 자승법
W(가중치)는 대칭이면서 양의 정부호성을 갖는 행렬이다. 만약 어떤 측정값이 다른 값보다 더 확실하다거나 불확실하면 W를 통해 계산상에 가중치를 주어 결과의 신뢰도를 높인다.
최소 자승법에 의한 곡선의 근사
주어진 데이터의 개수가 구하려는 곡선의 매개변수보다 많은 경우 정확한 곡선은 구할 수가 없지만 구해진 곡선과 주어진 데이터와의 오차가 가장 적은 곡선은 최소 자승법을 이용하여 구한다. 여기서 최소 자승법을 이용하여 계수를 구하면 이 계수들로 이루어진 곡선은 주어진 데이터와 가장 적은 오차를 가지게 된다.
최소 자승법은 작은 오차가 여러 개 보이는 것보다 큰 오차 하나에 의해서 더 많은 영향을 받는다. 따라서 되도록이면 큰 오차에 가능성이 있는 데이터는 피하는 것이 좋다. 사전에 데이터에 대한 정보가 있으면 가중 최소 자승법을 사용하고, 아니면 큰 오차의 가능성이 있는 데이터는 분리하는 것이 좋다
<수식>
#include <stdio.h>
void matrix_inversion_using_elementary_operation(float original[20][20],
float inverse[20][20], int dim) {
int row, column, pivot_column, max_index;
float max_value, ftemp1, ftemp2, pivot_value;
for (row = 0; row < dim; row++)
for (column = 0; column < dim; column++) {
if (row == column)
inverse[row][column] = 1;
else inverse[row][column] = 0;
}
for (pivot_column = 0; pivot_column < dim; pivot_column++) {
max_index = original[0][column];
max_value = 0;
for (row = pivot_column; row < dim; row++)
if (original[row][pivot_column] * original[row][pivot_column] > max_value *
max_value) {
max_index = row;
max_value = original[row][pivot_column];
}
if (pivot_column != max_index)
for (column = 0; column < dim; column++) {
ftemp1 = original[pivot_column][column];
ftemp2 = inverse[pivot_column][column];
original[pivot_column][column] = original[max_index][column];
inverse[pivot_column][column] = inverse[max_index][column];
original[max_index][column] = ftemp1;
inverse[max_index][column] = ftemp2;
}
pivot_value = original[pivot_column][pivot_column];
for (column = 0; column < dim; column++) {
original[pivot_column][column] /= pivot_value;
inverse[pivot_column][column] /= pivot_value;
}
for (row = 0; row < dim; row++)
if (row != pivot_column) {
ftemp1 = original[row][pivot_column];
for (column = 0; column < dim; column++) {
original[row][column] -= ftemp1 * original[pivot_column][column];
inverse[row][column] -= ftemp1 * inverse[pivot_column][column];
}
}
}
}
void main() {
int i, j, k, num_eq, dim;
float coeff[20][20], known[20];
float original[20][20], inverse[20][20], Atb[20], var[20];
printf("enter the number of equations.:");
scanf_s("%d", &num_eq);
printf("enter the number of variables.:");
scanf_s("%d", &dim);
for (i = 0; i < num_eq; i++) {
printf("**************************************\n");
printf("enter the coefficients of eq. %d. \n", i + 1);
for (j = 0; j < dim; j++) {
printf("coefficient of variable %d. :", j + 1);
scanf_s("%f", &coeff[i][j]);
}
printf("enther the known value of eq. %d :", i + 1);
scanf_s("%f", &known[i]);
}
printf("************************************\n");
printf("\n\n enter equations : \n");
for (i = 0; i < num_eq; i++) {
for (j = 0; j < dim; j++) {
printf("%10.3f*X%d", coeff[i][j], j + 1);
if (j != dim - 1)printf("+");
else printf("=");
}
printf("%10.3f \n", known[i]);
}
for (i = 0; i < dim; i++)
for (j = 0; j < dim; j++) {
original[i][j] = 0.0;
for (k = 0; k < num_eq; k++)
original[i][j] += coeff[k][i] * coeff[k][j];
}
matrix_inversion_using_elementary_operation(original, inverse, dim);
for (i = 0; i < dim; i++) {
Atb[i] = 0.0;
for (j = 0; j < num_eq; j++)
Atb[i] += coeff[j][i] * known[j];
}
for (i = 0; i < dim; i++) {
var[i] = 0.0;
for (j = 0; j < dim; j++)
var[i] += inverse[i][j] * Atb[j];
}
printf("\n\n**************************************************\n");
printf("results : \n");
for (i = 0; i < dim; i++)
printf("X%d=%10.3f\n", i + 1, var[i]);
}
순환 최소 자승법
새로운 데이터가 들어올 때마다 전체의 방정식을 모두 다루어야 하는 일반적인 최소 자승법의 불편함을 없앤 방식이다. 즉 이전까지의 데이터에 의한 최소 자승법의 결과와 새로 들어온 데이터를 이용해서 이 데이터까지 포함한 새로운 최소 자승법에 의한 결과를 구하는 것이다.
먼저 임의의 개수의 방정식을 받아들여서 그 최소 자승 결과를 구한 수 계속 방정식을 받아들여서 결과를 개선하도록 만들었다.
<수식>
#include <stdio.h>
void matrix_inversion_using_elementary_operation(float AtA_matrix[20][20],
float inv_AtA[20][20], int num_var) {
int row, column, pivot_column, max_index;
float max_value, ftemp1, ftemp2, pivot_value;
for (row = 0; row < num_var; row++)
for (column = 0; column < num_var; column++) {
if (row == column)
inv_AtA[row][column] = 1;
else inv_AtA[row][column] = 0;
}
for (pivot_column = 0; pivot_column < num_var; pivot_column++) {
max_index = AtA_matrix[0][column];
max_value = 0;
for (row = pivot_column; row < num_var; row++)
if (AtA_matrix[row][pivot_column] * AtA_matrix[row][pivot_column] >
max_value * max_value) {
max_index = row;
max_value = AtA_matrix[row][pivot_column];
}
if (pivot_column != max_index)
for (column = 0; column < num_var; column++) {
ftemp1 = AtA_matrix[pivot_column][column];
ftemp2 = inv_AtA[pivot_column][column];
AtA_matrix[pivot_column][column] = AtA_matrix[max_index][column];
inv_AtA[pivot_column][column] = inv_AtA[max_index][column];
AtA_matrix[max_index][column] = ftemp1;
inv_AtA[max_index][column] = ftemp2;
}
pivot_value = AtA_matrix[pivot_column][pivot_column];
for (column = 0; column < num_var; column++) {
AtA_matrix[pivot_column][column] /= pivot_value;
inv_AtA[pivot_column][column] /= pivot_value;
}
for (row = 0; row < num_var; row++)
if (row != pivot_column) {
ftemp1 = AtA_matrix[row][pivot_column];
for (column = 0; column < num_var; column++) {
AtA_matrix[row][column] -= ftemp1 * AtA_matrix[pivot_column][column];
inv_AtA[row][column] -= ftemp1 * inv_AtA[pivot_column][column];
}
}
}
}
void main() {
int input_key;
int i, j, k, num_eq, num_var;
float coeff[20][20], known[20];
float AtA_matrix[20][20], inv_AtA[20][20];
float Atb[20], var[20], pa[20], g[20], apa, ftemp;
printf("enter the number of equations. :");
scanf_s("%d", &num_eq);
printf("enter the number of variables. :");
scanf_s("%d", &num_var);
for (i = 0; i < num_eq; i++) {
printf("***********************************\n");
printf("enter the coefficients of eq. %d. \n", i + 1);
for (j = 0; j < num_var; j++) {
printf("coefficient of variable %d. : ", j + 1);
scanf_s("%f", &coeff[i][j]);
}
printf("enter the known value of eq. %d. : ", i + 1);
scanf_s("%f", &known[i]);
}
printf("*************************************\n");
printf("\n\n entered equations : \n");
for (i = 0; i < num_eq; i++) {
for (j = 0; j < num_var; j++) {
printf("%10.3f*X%d", coeff[i][j], j + 1);
if (j != num_var - 1)printf("+");
else printf("=");
}
printf("%10.3f \n", known[i]);
}
for (i = 0; i < num_var; i++)
for (j = 0; j < num_var; j++) {
AtA_matrix[i][j] = 0.0;
for (k = 0; k < num_eq; k++)
AtA_matrix[i][j] += coeff[k][i] * coeff[k][j];
}
matrix_inversion_using_elementary_operation(AtA_matrix, inv_AtA, num_var);
for (i = 0; i < num_var; i++) {
Atb[i] = 0.0;
for (j = 0; j < num_eq; j++)
Atb[i] += coeff[j][i] * known[j];
}
for (i = 0; i < num_var; i++) {
var[i] = 0.0;
for (j = 0; j < num_var; j++)
var[i] += inv_AtA[i][j] * Atb[j];
}
printf("\n\n******************************\n");
printf("result : \n");
for (i = 0; i < num_var; i++)
printf("X%d=%10.3f\n", i + 1, var[i]);
input_key = 1;
while (input_key != 27) {
printf("\n\n********************************\n");
printf("entered the coefficient of the new equation..\n");
for (i = 0; i < num_var; i++) {
printf("coefficient of variable %d. :", i + 1);
scanf_s("%f", &coeff[0][i]);
}
printf("enter the known value of the new equation. :");
scanf_s("%f", &known[0]);
printf("\n\n*************************************\n");
for (i = 0; i < num_var; i++) {
pa[i] = 0;
for (j = 0; j < num_var; j++)
pa[i] += inv_AtA[i][j] * coeff[0][j];
}
apa = 0;
for (i = 0; i < num_var; i++)
apa += coeff[0][i] * pa[i];
for (i = 0; i < num_var; i++)
g[i] = pa[i] / apa;
ftemp = known[0];
for (i = 0; i < num_var; i++)
ftemp -= coeff[0][i] * var[i];
for (i = 0; i < num_var; i++)
var[i] += g[i] * ftemp;
printf("\n\n***************************************\n");
printf("the new results are ............. \n");
for (i = 0; i < num_var; i++)
printf("X%d=%10.3f \n", i + 1, var[i]);
printf("\n\n****************************************\n");
printf("press any key to continue...ESC to terminate \n");
input_key = _getch();
}
}'학교수업, CS > 수치해석' 카테고리의 다른 글
| 난수(변환법, 거부법, 가우시안 분포, 몬테카를로 적분) (0) | 2022.11.30 |
|---|---|
| 보간법(선형, 랑그라제, 네빌레, 뉴튼 다항식, 전/후향 차분법, 스플라인) (0) | 2022.11.30 |
| 고유치(파데브-레브리어의 알고리즘) (0) | 2022.11.30 |
| 행렬&역행렬 (0) | 2022.11.30 |
| 연립방정식(가우스 조단의 피보팅, 가우스 자이달 방법) (1) | 2022.11.30 |