Interpolation : Part 2

안녕하십니까 Mr. K입니다

지난주에 한회 쉬었음에도 제 날짜에 포스팅을 하지 못해 뭐라 할 말이 없습니다 -_-;

지난 포스트(Interpolation : Part 1)에 이어서

같은 데이터( 3개의 점: (1, 2), (2, 3), (3, 6) )를 가지고 얘기해보도록 하겠습니다

이번에 얘기할 주제는 Newton Forward Differences 입니다

We shall now discuss the operation of fitting an Nth order polynomial through N + 1 points and remark again this will yield a unique answer. We consider the set of data points: (x_0, f_0), (x_1, f_1), …, (x_N, f_N). We now introduce the difference operator Δ such that

or in general

These are called forward differences, since points forward of the current value are used (there are also backward differences and central differences, which we shall meet in our discussion of the solution of differential equations). We consider the composition of the operator whereby

Of course we can now proceed to define Δⁿf_0 in an iterative manner.
We introduce the polynomial

which will ultimately terminate after N terms. Alternatively this can be written as:

… (5.1)

We consider the data values to be equally spaced, so the value of Δx_j is h ∀j. This analysis can be extended to irregularly spaced points, but we shall not attempt that here.
We are now able to construct the polynomial (5.1) for a set of points. We start with a simple …

(위 내용은 AN INTRODUCTION TO Programming and Numerical Methods in MATLAB에서 일부 발췌했음을 알려드립니다)

책에 쓰인 글을 그대로 옮겨온거라 오타가 있을지도.. -_-;

어렵게 쓰여있는 듯 하지만 실은 별거 아닙니다

먼저 difference operator Δ(delta라고 읽습니다)를 정의하는데, 이것은 주어진 순서쌍의 집합에서 '다음값 - 현재값'으로 정의합니다

그리고 Δ의 중첩은 세번째 그림와 같이 정의합니다

그리고나서 네번째 그림과 같은 다항식을 생각하게 되는데,

이것이 바로 우리가 구하고자 하는 '주어진 N+1개의 점을 지나는 N차 다항식'입니다

그것을 줄여서 쓴 것이 다섯번째 그림입니다

이때 h라는 변수는

주어진 순서쌍의 집합에서 x좌표를 기준으로 오름차순 정렬을 했다고 가정하고

'다음 x값 - 현재 x값'이 일정한 상수값을 가질 때 그 상수값을 h에 대입합니다

지난 포스트에서 미리 잡아둔 데이터: (1, 2), (2, 3), (3, 6)를 가지고 이것을 적용해보도록 하지요

우선 세 점의 x좌표는 1, 2, 3인데 이것은 오름차순으로 정렬되어있고, 따라서 '다음 x값 - 현재 x값'은 1이라는 고정된 값을 주게 됩니다

따라서 h = 1이고,

f_0와 Δf_0, 그리고 Δ²f_0는 위의 표에 계산이 되어있으니 공식에 적용해보도록 하겠습니다

간단하게 나오지요?

다만 단점이라면, h라는 변수를 정의하기 위해 주어진 데이터의 x값은 균일해야 한다는 것입니다

(위 지문의 마지막 부분을 보면 균일하지 않은 점에 대해서도 확장해서 적용할 수 있으나 여기서는 하지 않겠다라며 넘어갑니다)

이것을 코드로 쓰면 다음과 같습니다

#include <stdio.h> #include <stdlib.h> double difference( double *f, int degree, int prs ); double combination( double *xs, int n, int exp ); double ipow( double m, int n ); double fact( int n ); void main() { int pairs; int i; int j; double *xs; // 각 점의 x좌표를 모아놓은 array double *ys; // 각 점의 y좌표를 모아놓은 array double *coef; double temp; double h; printf("점의 개수를 입력 > "); scanf("%d", &pairs); printf("\n"); xs = ( double* )malloc( sizeof( double ) * pairs ); ys = ( double* )malloc( sizeof( double ) * pairs ); for( i = 0; i < pairs; i++ ) { printf("%2d번째 점의 x좌표와 y좌표를 순서대로 입력 > ", i+1); scanf("%lf %lf", &xs[i], &ys[i]); } h = xs[1] - xs[0]; coef = ( double* )malloc( sizeof( double ) * pairs ); // coef를 계산 for( i = 0; i < pairs; i++ ) { temp = 0; for( j = 1; j < pairs; j++ ) { temp += difference( ys, j, pairs ) * combination( xs, j, i ) / ( ipow( h, j ) * fact(j) ); } coef[i] = temp; } coef[0] += ys[0]; // coef를 출력 printf("\n주어진 점을 지나는 %d차 다항식의 계수는 다음과 같습니다\n", pairs-1); for( i = 0; i < pairs; i++ ) { if( !i ) printf("상수항 = "); else printf("%2d차항 = ", i); printf("%10.4lf\n", coef[i]); } free(xs); free(ys); free(coef); } double difference( double *f, int degree, int prs ) { if( degree == 1 ) { return ( f[1] - f[0] ); } else { int i; double *diff = ( double * )malloc( sizeof( double ) * ( prs - 1 ) ); double res; for( i = 0; i < degree; i++ ) { diff[i] = f[i+1] - f[i]; } res = difference( diff, degree-1, prs ); free( diff ); return res; } } // xs : 각 점의 x좌표를 모아놓은 array // n : 배열의 일부를 사용하는데, 그 size // exp : (x^exp)의 계수를 구하기 위한 지수 double combination( double *xs, int n, int exp ) { double res = 0; double temp; int *dummy; // 반복문에 쓰일 n-exp개의 index int i; int j; int r; // 조합 nCr에서의 r r = n - exp; if( r == 0 ) return 1; if( r < 0 ) return 0; dummy = ( int* )calloc( r, sizeof( int ) ); for( i = 0; i < r; i++ ) { dummy[i] = i; } while( dummy[0] < exp+1 ) { temp = 1; for( i = 0; i < r; i++ ) { temp *= -xs[dummy[i]]; } res += temp; dummy[r-1]++; for( i = r-1, j = 0; i > 0; i--, j++ ) { if( dummy[i] > n-1-j ) { dummy[i-1]++; dummy[i] = dummy[i-1] + 1; } } } return res; } double ipow( double m, int n ) { if( m == 1 || n == 0 ) return 1; else { double res = 1; int i; for( i = 1; i <= n; i++ ) { res *= m; } return res; } } double fact( int n ) { if( n <= 1 ) return 1; else { double res = 1; int i; for( i = 2; i <= n; i++ ) { res *= i; } return res; } }

이번 포스트는 이것으로 마치고,

다음주에 쓰게 될 포스트에서는 데이터의 x값이 균일하지 않아도 쓸 수 있는 Lagrange polynomial에 대해서 얘기하겠습니다

(포스팅 소요시간 : 약 75분)