선형 보간법
선형 보간법은 여러 기본적인 수치해석 방법의 기초가 된다. 선형 보간법을 이용하여 사다리꼴 법이라는 한 적분 방법도 유도해 낼 수 있으며 미분 방정식을 푸는데도 사용된다. 선형 보간법은 단순히 주어진 두 점을 이은 직선으로 구하는 방법이다.
<수식>
#include <stdio.h>
#include <math.h>
#include <conio.h>
#include <windows.h>
void clrscr(void)
{
COORD Cur = { 0, 0 };
unsigned long dwLen;
FillConsoleOutputCharacter(GetStdHandle(STD_OUTPUT_HANDLE), ' ', 80 * 25, Cur, &dwLen);
}
double linear_interpol(double* x1, double* x2, double* y1, double* y2, double* x) {
double result;
result = (*y2 - *y1) / (*x2 - *x1) * (*x - *x1) + *y1;
return result;
}
int main(void) {
double x, x1, x2, y1, y2, output;
x1 = 1.0;
x2 = 2.0;
y1 = exp(x1);
y2 = exp(x2);
clrscr();
printf("\n------\t Interpolated value ------ Real value -------------\n");
for (x = 1.0; x <= 2.0; x += 0.1) {
output = linear_interpol(&x1, &x2, &y1, &y2, &x);
printf("\n\t f(%lf)=%lf\t exp(%lf)=%lf", x, output, x, exp(x));
}
printf("\n");
_getch();
}
랑그라제 다항식을 이용한 보간법
점과 점을 직선이 아닌 어떤 곡선으로 연결하여 구하는 방법
<수식>
#include <stdio.h>
#include <math.h>
#include <conio.h>
#include <windows.h>
#define DATA 4
void clrscr(void)
{
COORD Cur = { 0, 0 };
unsigned long dwLen;
FillConsoleOutputCharacter(GetStdHandle(STD_OUTPUT_HANDLE), ' ', 80 * 25, Cur, &dwLen);
}
double linear_interpol(double* x1, double* x2, double* y1, double* y2, double* x) {
double result;
result = (*y2 - *y1) / (*x2 - *x1) * (*x - *x1) + *y1;
return result;
}
int main(void) {
int i, j, k;
double result, numerator, denominator;
double x[] = { 1.5,2.7,3.4 };
double in[] = { 1.,2.,3.,4. };
double out[] = { 0.0,0.30103,0.47712,0.60206 };
clrscr();
printf("\n-------\t interpolated value ----- real value --------------\n");
for (i = 0; i < 3; i++) {
result = 0.0;
for (j = 0; j < DATA; j++) {
numerator = 1.0;
denominator = 1.0;
for (k = 0; k < DATA; k++) {
if (k == j)continue;
numerator *= (x[i] - in[k]);
denominator *= (in[j] - in[k]);
}
result += numerator / denominator * out[j];
}
printf("\n\tf(%lf)=%lf\t log(%lf)=%lf", x[i], result, x[i], log10(x[i]));
}
printf("\n");
_getch();
}
네빌레의 반복 보간법
랑그라제 보간 다항식의 단줌중 하나가 데이터 중 새로운 점이 하나만 더 추가되어도 그 앞에서 구했던 계산 결과를 이용할 수 없다는 점이다. 새로운 점이 추가되면 새로운 식을 처음부터 다시 구해야 하는 단점이 있는데, 네빌레의 반복보간법은 앞에서 구했던 계산이나 결과를 다음 단계에서 활용함으로써 보다 편하게 보간 다항식을 구할 수 있다.
<수식>
#include <stdio.h>
#include <math.h>
#include <conio.h>
#include <windows.h>
#define DATA 5
void clrscr(void)
{
COORD Cur = { 0, 0 };
unsigned long dwLen;
FillConsoleOutputCharacter(GetStdHandle(STD_OUTPUT_HANDLE), ' ', 80 * 25, Cur, &dwLen);
}
int main(void) {
int i, j;
float in[5], out[5][5], x;
/*initialization*/
for (i = 0; i < DATA; i++) {
in[i] = (float)(i + 1);
out[0][i] = sqrt((float)(i + 1));
}
clrscr();
printf("\n Enter the x-axis value:");
scanf_s("%f", &x);
/*nevile alogorithm*/
for (j = 0; j < DATA-1; j++)
for (i = 0; i < DATA - j - 1; i++)
out[j + 1][i];
((x - in[i]) * out[j][i + 1] - (x - in[i + j + 1]) * out[j][i]) /
(in[i + j + 1] - in[i]);
/*print out the result*/
printf("\n the interpolated out put of input %f is %f\n", x, out[0][1]);
printf("the real value of input %f is %f\n", x, sqrt((float)x));
}
뉴튼 다항식에 의한 보간법
랑그라제 보간법의 단점
- 하나의 보간에 필요로 하는 계산량이 많다
- 데이터의 수가 증가할 때 그 전의 사용 결과가 사용되지 않는다
- 에러 계산이 용이하지 않다
뉴튼 보간법은 이 단점을 극복한다.
- 뉴튼 보간법은 주어진 데이터를 기초로 하여 차분표를 적성한다
- 차분표가 작성되면 보간공식을 쉽게 구할 수 있다. 따라서 보간 다항식의 차수는 데이터가 추가되어도 쉽게 그 차수를 늘리는 장점이 있다.
<수식>
#include <stdio.h>
#include <math.h>
#include <conio.h>
#include <windows.h>
#define DATA 5
void clrscr(void)
{
COORD Cur = { 0, 0 };
unsigned long dwLen;
FillConsoleOutputCharacter(GetStdHandle(STD_OUTPUT_HANDLE), ' ', 80 * 25, Cur, &dwLen);
}
int main(void) {
int i, j;
double temp, result, x;
double in[] = { 3.,3.3,3.5,3.7,4. }, out[DATA][DATA] =
{ {1.09861,1.19392,1.25276,1.30833,1.335}, };
clrscr();
printf("\n enter the x-axis value: ");
scanf_s("%lf", &x);
for (j = 1; j < DATA; j++)
for (i = 0; i < DATA - j; i++)
out[j][i];
(out[j - 1][i + 1] - out[j - 1][i]) / (in[i + j] - in[i]);
result = 0.;
for (j = 0; j < DATA; j++) {
temp = 1.;
for (i = 0; i < j; i++)
temp *= (x - in[i]);
result += temp * out[j][0];
}
printf("\n the interpolated output of input %lf is %lf\n", x, result);
}
전향 차분법
<수식>
#include <stdio.h>
#include <math.h>
#include <conio.h>
#include <windows.h>
#define DATA 5
void clrscr(void)
{
COORD Cur = { 0, 0 };
unsigned long dwLen;
FillConsoleOutputCharacter(GetStdHandle(STD_OUTPUT_HANDLE), ' ', 80 * 25, Cur, &dwLen);
}
double binomial(double s, int j) {
int m, n;
double temp, bino;
if (j == 0)bino = 1.;
else if (j == 1)bino = s;
else {
m = j;
temp = s;
for (n = 1; n < m; n++) {
temp *= (s - (double)n);
j *= n;
}
bino = temp / (double)j;
}
return bino;
}
int main(void) {
int i, j;
double result, x, s=0;
double in[] = { 10.,11.,12.,13.,14. }, out[DATA][DATA] =
{ {3.162278,3.316626,3.464102,3.605551,3.741657}, };
printf("\n enter the x-axis value:");
scanf_s("%lf", &x);
/*algorithm implementation*/
/*finding h and s
h=in[1]-in[0]
s=(x-in[0]/h */
/*finding result*/
result = 0.;
for (j = 0; j <= 4; j++) {
for (i = j + 1; i <= 4; i++)
out[j + 1][i] = out[j][i] - out[j][i - 1];
result += binomial(s, j) * out[j][j];
}
printf("\n the interpolated output of input %lf is %lf\n", x, result);
}
후향 차분법
<수식>
#include <stdio.h>
#include <math.h>
#include <conio.h>
#include <windows.h>
#define DATA 5
void clrscr(void)
{
COORD Cur = { 0, 0 };
unsigned long dwLen;
FillConsoleOutputCharacter(GetStdHandle(STD_OUTPUT_HANDLE), ' ', 80 * 25, Cur, &dwLen);
}
double binomial(double s, int j) {
int m, n;
double temp, bino;
if (j == 0)bino = 1.;
else if (j == 1)bino = s;
else {
m = j;
temp = s;
for (n = 1; n < m; n++) {
temp *= (s - (double)n);
j *= n;
}
bino = temp / (double)j;
}
return bino;
}
int main(void) {
int i, j;
double result, x, s=0;
double in[] = { 1.0,1.3,1.6,1.9,2.2 }, out[DATA][DATA] =
{ {0.7651977,0.6200860,0.4554022,0.2818186,0.1103623}, };
clrscr();
printf("\n enter the value of x-domain :");
scanf_s("%lf", &x);
/*algorithm implementation*/
/*finding h and s
h=in[1]-in[0]
s=(x-in[DATA-1])/h */
/*finding result*/
result = 0.;
for (j = 0; j <DATA; j++) {
for (i = j + 1; i <DATA; i++)
out[j + 1][i] = out[j][i] - out[j][i - 1];
result += binomial(-s, j) * out[j][DATA-1]*pow(-1.0,j);
}
printf("\n the interpolated output of input %lf is %lf\n", x, result);
}
3차원 스플라인 보간법
종종 많은 데이터들이 하나의 매끄러운 단일 곡선의 점들일 수가 있다. 만일 이러한 경우에 주어진 데이터들을 라그랑제 보간법이나 뉴튼의 방법으로 보간한다면 차수가 늘어남에 따라 더욱 큰 오차가 발생하게 될 것이다. 이러한 경우에는 3차원 스플라인 보간법이 최적의 방법이다
3차원 스플라인 보간법에는 두 개의 이웃하는 데이터 사이에 3차 다항식이 사용된다. 이 다항식은 4개의 미정 계수가 있으므로 4가지 조건을 필요로 한다. 4개의 미정 계수들 중 2개는 이 다항식이 두 점을 반드시 통과해야 한다는 조건에서 구할 수 있고, 나머지 2개는 1차 미분, 2차 미분한 값이 주어진 두 점에서 연속해야 한다는 사실로부터 구할 수 있다.
<수식>
#include <stdio.h>
#include <math.h>
#include <conio.h>
#include <windows.h>
#define DATA 9
void clrscr(void)
{
COORD Cur = { 0, 0 };
unsigned long dwLen;
FillConsoleOutputCharacter(GetStdHandle(STD_OUTPUT_HANDLE), ' ', 80 * 25, Cur, &dwLen);
}
void tridia(double element_1[], double element_2[], double element_3[],
double element_4[], double *s) {
int i;
double temp;
for (i = 1; i < DATA - 2; i++) {
temp = element_1[i] / element_2[i - 1];
element_2[i] -= temp * element_3[i - 1];
element_4[i] -= temp * element_4[i - 1];
}
s[DATA - 3] = element_4[DATA - 3] / element_2[DATA - 3];
for (i = DATA - 4; i >= 0; i--)
s[i] = (element_4[i] - element_3[i] * s[i + 1] / element_2[i]);
}
void interpol(double in[], double out[], double s[], double h[], double inv) {
int i;
double temp, a[DATA], b[DATA], c[DATA], d[DATA], result;
for(i = 0; i < DATA - 1; i++) {
a[i] = (s[i + 1] - s[i]) / (6.0 * h[i]);
b[i] = s[i] / 2.0;
c[i] = (out[i + 1] - out[i]) / h[i] - (2.0 * h[i] * s[i] + h[i] * s[i + 1]) /
6.0;
d[i] = out[i];
}
temp = 0.0;
do {
for (i = 0; i < DATA - 1; i++) {
if ((temp >= in[i]) && (temp < in[i + 1])) {
result = d[i] + (temp - in[i] * (c[i] + (temp - in[i]) * (b[i] + a[i] *
(temp - in[i]))));
printf("\n F[%lf]=%lf", temp, result);
break;
}
}
temp += inv;
} while (temp < in[DATA - 1]);
}
int main(void) {
int i, cond;
double h[DATA - 1], element_1[DATA], element_2[DATA], element_3[DATA], element_4[DATA],
s[DATA], inv;
double in[] = { 0.0, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 1.0 }, out[] =
{ 0.302,0.185,0.106,0.093,0.240,0.579,0.561,0.468,0.302 };
clrscr();
printf("\n type the output interval of x-domain that you wanna find: ");
scanf_s("%lf", &inv);
for (i = 0; i < DATA - 1; i++)
h[i] = in[i + 1] - in[i];
for (i = 1; i < DATA - 3; i++) {
element_1[i] = h[i];
element_2[i] = 2. * (h[i] + h[i + 1]);
element_3[i] = h[i + 1];
}
for (i = 0; i < DATA - 2; i++)
element_4[i] = 6.0 * ((out[i + 2] - out[i + 1]) / h[i + 1] -
(out[i + 1] - out[i]) / h[i]);
for (cond = 0; cond < 3; cond++) {
if (cond == 0) {
s[0] = s[DATA - 1] = 0.0;
element_1[DATA - 3] = h[DATA - 3];
element_2[0] = 2.0 * (h[0] + h[1]);
element_2[DATA - 3] = 2.0 * (h[DATA - 3] + h[DATA - 2]);
element_3[0] = h[1];
tridia(element_1, element_2, element_3, element_4, s);
printf("\n for boundary condition 1 -----------------------");
interpol(in, out, s, h, inv);
}
else if (cond == 1) {
element_1[DATA - 3] = h[DATA - 3];
element_2[0] = 3.0 * (h[0] + h[1]);
element_2[DATA - 3] = 3.0 * (h[DATA - 2] + 2.0 * h[DATA - 3]);
element_3[0] = h[1];
tridia(element_1, element_2, element_3, element_4, s);
s[0] = s[1], s[DATA - 1] = s[DATA - 2];
printf("\n for boundary condition 2 -----------------------");
interpol(in, out, s, h, inv);
}
else {
element_2[0] = (h[0] + h[1]) * (h[0] + 2. * h[1]) /
h[1] * (h[1] * h[1] - h[0] * h[0]) / h[1];
element_2[DATA - 3] = (h[DATA - 2] + h[DATA - 3]) *
(h[DATA - 2] + 2. * h[DATA - 3] )/ h[DATA - 3];
element_1[DATA - 3] = (h[DATA - 3] * h[DATA - 3] - h[DATA - 2] * h[DATA - 2]) /
h[DATA - 3];
element_3[0]= ( h[1] * h[1] - h[0] * h[0]) / h[1];
tridia(element_1, element_2, element_3, element_4, s);
s[0]= ((h[0] + h[1]) * s[1] - h[0]*s[2]) / h[1];
s[DATA-1]= ((h[DATA - 3] + h[DATA - 2]) * s[DATA - 2] - h[DATA - 2]*s[DATA-3]) /
h[DATA - 3];
printf("\n for boundary condition 2 -----------------------");
interpol(in, out, s, h, inv);
}
}
}'학교수업, CS > 수치해석' 카테고리의 다른 글
| 푸리에 변환(DFT, FFT, 파워스펙트럼 추정, 자기상관도) (0) | 2022.11.30 |
|---|---|
| 난수(변환법, 거부법, 가우시안 분포, 몬테카를로 적분) (0) | 2022.11.30 |
| 최소 자승법(촐레스키법을 이용한 행렬의 제곱근) (0) | 2022.11.30 |
| 고유치(파데브-레브리어의 알고리즘) (0) | 2022.11.30 |
| 행렬&역행렬 (0) | 2022.11.30 |