Interpolation : Part 3

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

목요일에 글을 쓰기 시작하는것도 오랜만이군요

과연 목요일이 끝나기 전에 글을 올릴 수 있느냐가 문제지만 -_-;

지난 두 포스트에 이어서 계속 같은 데이터( 3개의 점: (1, 2), (2, 3), (3, 6) )를 가지고 얘기해보겠습니다

…. We remark that various conclusions can be drawn from the data by using the forward differences; for instance whether it was generated using a polynomial (in which case the differences truncate).
We now construct a polynomial which goes through a set of points which are not necessarily evenly spaced. Let us consider the polynomial

It is worth pausing at this point and checking that this curve goes through each of the points (x_j, f_j). For example for j = 3, we set z = x_3 and only the third term is non-zero and we have

Hence the value of the polynomial at z = x_3 is f_3 (as we would hope). This is an example of a Lagrange polynomial; which could have equally been written as

This is a convenient way to write out the cubic we require (as it is relatively easily extended to higher-order cases) and in order to evaluate it we can use the MATLAB code: …

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

영어로 뭐라고 막 쓰여있지만 가볍게 무시해줍시다 -_-;

주어진 순서쌍의 집합에는 (x_n, f_n)꼴의 점들이 있는데, '주어진 N+1개의 점을 지나는 N차 다항식'을 구하기 위해

첫번째 그림에 있는 것과 같은 모양의 다항식을 생각합니다

그러면 저 다항식은 두번째 그림의 예와 같이 f(x_n) = f_n 이라는 관계를 갖게 됩니다

이러한 다항식을 Lagrange polynomial이라고 부르고 세번째 그림과 같이 줄여서 씁니다

책에서는 4개의 점을 예로 들어서 다항식을 소개하였는데,

이것을 'N개'의 점에 대해서 확장시키면 다음과 같습니다

N개의 점에 대한 N-1차 Lagrange polynomial.

N+1개의 점에 대해서 확장시키려면 위의 N 자리에 N+1을 대입하기만 하면 됩니다

책에서는 N개의 점에 대해서 소개하고 있지만, 저는 보간법 포스팅을 N+1개의 점에 대해서 시작했기 때문에

설명하는 내용과 등장하는 수식 간에 다소 차이가 있을 수 있습니다

이 공식을 우리의 데이터: (1, 2), (2, 3), (3, 6)으로 적용해보면 다음과 같습니다

다소 계산은 복잡합니다만

지난주에 포스팅했던 Newton Forward Differences 와 비교했을 때 훨씬 '단순한' 모양을 하고 있습니다

수업시간에 들은 내용이 잘 기억나지 않지만, 이 다항식의 단점은 딱히 없었던 것 같습니다

굳이 단점이라 할만한 것이 있다면 손으로 계산하기 귀찮다는 것이랄까 -_-;

(개인적인 생각이지만, Vandermonde matrix >> Lagrange polynomial > Newton Forward Differences 순으로 귀찮은 연산인 것 같습니다 -_-;)

코드로 구현해보면 다음과 같습니다

#include <stdio.h> #include <stdlib.h> double combination( double *xs, int n, int exp ); double numer( double *xs, int n, int prs, int exp ); double denom( double *xs, int prs, int exc ); void main() { int pairs; int i; int j; double *xs; // 각 점의 x좌표를 모아놓은 array double *ys; // 각 점의 y좌표를 모아놓은 array double *coef; double temp; 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]); } coef = ( double* )malloc( sizeof( double ) * pairs ); // coef를 계산 for( i = 0; i < pairs; i++ ) { temp = 0; for( j = 0; j < pairs; j++ ) { temp += numer( xs, i, pairs, j ) * ys[j] / denom( xs, pairs, j ); } coef[i] = temp; } // 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); } // 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; } // xs : 각 점의 x좌표를 모아놓은 array // exp : (x^exp)의 계수를 구하기 위한 지수 // prs : 모든 점의 개수 // exc : 특정 연산을 수행하지 않는 index double numer( double *xs, int exp, int prs, int exc ) { double res; double *newxs; int i; int j; newxs = ( double* )malloc( sizeof( double ) * ( prs - 1 ) ); for( i = 0, j = 0; i < prs; i++ ) { if( i == exc ) continue; else { newxs[j] = xs[i]; j++; } } res = combination( newxs, prs-1, exp ); return res; } // xs : 각 점의 x좌표를 모아놓은 array // prs : 모든 점의 개수 // exc : 특정 연산을 수행하지 않는 index double denom( double *xs, int prs, int exc ) { double res = 1; int i; for( i = 0; i < prs; i++ ) { if( i == exc ) continue; else { res *= xs[exc] - xs[i]; } } return res; }

위의 소스에서 combination이라는 함수는 Newton Forward Differences에 썼던 것을 그대로 가져왔습니다

numer라는 함수는 다항식의 계수부분 중 분자(numerator) 부분을 구하는 함수이고,

(정확히 말하면 분자를 모두 구하는 것은 아니고, 위의 공식에서 f_i 를 제외한 나머지부분에 대한 분자를 구합니다)

denom이라는 함수는 다항식의 계수부분 중 분모(denominator) 부분을 구하는 함수입니다

이제 우리의 데이터: (1, 2), (2, 3), (3, 6)을 가지고 프로그램을 실행한 결과는 아래와 같고,

다른 case를 테스트하기 위해 적당히 4개의 점을 입력한 결과는 아래와 같습니다

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

다음주는 시험기간이라 한 회 쉬게 될 것 같습니다

더불어, 요즘 pku 3685: matrix 를 풀어보고있는데, 잘 안됩니다

AC받은 2人은 시간 나는대로 서로 상의를 하든지 해서 나 군입대 하기 전까지 풀이를 올려주길//

Talk카테고리에 미리 얘기했던 것처럼 명확한 이유가 없다면 태클 마구 들어갈겁니다 -_-;

(포스팅 소요시간 : 약 65분 + 서버점검 8시간 + 5분)