3차 스플라인 보간

CSAPI 명령

다음 명령은

values = csapi(x,y,xx)

not-a-knot 끝점 조건을 사용하여 주어진 데이터(x,y)에 대한 3차 스플라인 보간을 하여 xx에서의 값을 반환합니다. 여기서의 보간 함수는 절점 시퀀스 x를 갖는 조각별 3차 함수입니다. 이 함수의 3차 조각이 합쳐져 2개의 연속 도함수를 갖는 하나의 함수가 생성됩니다. "not-a-knot" 끝점 조건은 첫 번째 내부 절점과 마지막 내부 절점에서 3계 도함수도 연속적(반올림 오차까지)이라는 조건입니다.

데이터 점을 2개만 지정하면 직선 보간이 생성됩니다.

x = [0 1];
y = [2 0];
xx = linspace(0,6,121);
plot(xx,csapi(x,y,xx),'k-',x,y,'ro')
title('Interpolant to Two Points')

3개의 데이터 점을 지정하면 포물선 보간이 생성됩니다.

x = [2 3 5];
y = [1 0 4];
plot(xx,csapi(x,y,xx),'k-',x,y,'ro')
title('Interpolant to Three Points')

더 일반적으로, 4개 이상의 데이터 점을 지정하면 3차 스플라인 보간이 생성됩니다.

x = [1 1.5 2 4.1 5];
y = [1 -1 1 -1 1];
plot(xx,csapi(x,y,xx),'k-',x,y,'ro')
title('Cubic Spline Interpolant to Five Points')

CSAPI의 출력을 확인하는 방법

위의 보간은 훌륭한 보간처럼 보이지만, csapi가 제대로 동작하는지 확인하려면 어떻게 해야 할까요?

우리는 데이터 점을 플로팅하고 보간 함수가 이러한 점들을 관통하는 것을 보면서 csapi가 보간을 잘 수행했음을 확인했습니다. 하지만 3차 스플라인을 확실히 얻으려면, 예상되는 유형의 3차 스플라인에서 가져온 데이터로 시작하여 csapi가 그 3차 스플라인을 재현하는지, 즉 데이터를 가져온 바로 그 원래 3차 스플라인을 다시 제공하는지 확인해 봐야 할 수 있습니다.

예: 절단 멱함수

확인할 수 있는 3차 스플라인 함수의 간단한 예로 3차 절단 멱함수를 들 수 있습니다. 즉, 다음과 같습니다.

여기서 xi는 절점 중 하나이고 아래 첨자 "+"는 다음과 같이 명령 subplus에 의해 제공되는 절단 함수를 나타냅니다.

help subplus

 SUBPLUS Positive part.
 
                                   x , if  x>=0
    y  = subplus(x) := (x)_{+}  =               ,
                                   0 , if  x<=0
 
   returns the positive part of X. Used for computing truncated powers.

3차 절단 멱함수는 특정 선택 xi = 2에 대해 아래와 같이 플로팅됩니다. 예상한 바와 같이, 2의 왼쪽은 0이고 2의 오른쪽은 (x-2)^3처럼 상승합니다.

plot(xx, subplus(xx-2).^3,'y','LineWidth',3)
axis([0,6,-10,70])

이번에는 데이터 지점 0:6에서 이 특정 3차 스플라인을 보간하고, 스플라인 위에 보간 함수를 검은색으로 플로팅합니다.

x = 0:6;
y = subplus(x-2).^3;
values = csapi(x,y,xx);
hold on
plot(xx,values,'k',x,y,'ro')
hold off
title('Interpolant to ((x-2)_+)^3')

보간 오차

두 함수를 비교할 때는 대개 두 함수의 차이를 플로팅하는 것이 유익합니다.

plot(xx, values - subplus(xx-2).^3)
title('Error in Cubic Spline Interpolation to ((x-2)_+)^3')

그 차이의 크기를 가늠해 보기 위해 최대 데이터 값도 구해 볼 수 있습니다. 그러면 오차가 불가피한 반올림 오차보다는 더 낫다는 것을 알 수 있습니다.

max_y = max(abs(y))

max_y = 64

재현할 수 없는 절단 멱함수

추가적인 테스트로, 지점 0:6에서 csapi로 생성된 보간과 일치할 수 없는 절단 멱함수를 보간합니다. 예를 들어, csapi가 "not-a-knot" 조건을 사용하기 때문에 보간 스플라인의 첫 번째 내부 절점은 실제로 매듭이 아닙니다. 따라서 이 보간은 해당 지점에서 3개의 연속 도함수를 갖습니다. 즉, 3계 도함수가 그 지점에서 불연속이므로 그 지점을 중심으로 3차 절단 멱함수를 재현할 수 없음을 의미합니다.

values = csapi(x,subplus(x-1).^3,xx);
plot(xx, values - subplus(xx-1).^3)
title('Error in Not-a-Knot Interpolant to ((x-1)_+)^3')

1은 첫 번째 내부 매듭이므로, 이 보간에서는 활성화되지 않습니다.

그 차이의 크기는 0.18이지만, 1에서 멀어질수록 빠르게 감쇠됩니다. 위 그림은 3차 스플라인 보간이 본질적으로 국부적임을 보여줍니다.

값 대신 ppform 사용하기

3차 보간 스플라인을 그다음 값의 계산이나 도함수 계산 또는 기타 조작에 적합한 형식으로 유지할 수 있습니다. 이를 위해 다음 형식의 csapi를 호출합니다.

pp = csapi(x,y)

그러면 보간 함수의 ppform이 반환됩니다. 다음 명령으로 새로운 점 xx에서 이 형식의 값을 구할 수 있습니다.

values = fnval(pp,xx)

다음 명령으로 보간 함수를 미분할 수 있습니다.

dpp = fnder(pp)

또는 다음 명령으로 보간 함수를 적분할 수 있습니다.

ipp = fnint(pp)

그러면 각각 미분 또는 적분의 ppform이 반환됩니다.

예: 보간 함수 미분 및 적분하기

보간 함수의 미분을 표시하기 위해 다음 절단 멱함수의 도함수를 플로팅(노란색)한 다음

그 위에 원래 3차 절단 멱함수에 대한 보간 함수의 도함수를 플로팅(검은색)합니다.

plot(xx,3*subplus(xx-2).^2,'y','LineWidth',3)
pp = csapi(x,subplus(x-2).^3);
dpp = fnder(pp);
hold on
plot(xx,fnval(dpp,xx),'k')
hold off
title('Derivative of Interpolant to ((x-2)_+)^3')

여기서도 그 차이를 플로팅하여 비교하는 것이 더 유익하며, 전과 같이 그 차이가 반올림 오차보다 더 크지 않습니다.

plot(xx, fnval(dpp,xx) - 3*subplus(xx-2).^2)
title('Error in Derivative of interpolant to ((x-2)_+)^3')

절단 멱함수의 2계 도함수는 다음과 같습니다.

이 함수와 원래 함수에 대한 보간 함수의 2계 도함수 간의 차이를 플로팅하면 현재 비약이 발생하지만, 이는 여전히 반올림 오차 범위 이내입니다.

ddpp = fnder(dpp);
plot(xx, fnval(ddpp,xx) - 6*subplus(xx-2))
title('Error in Second Derivative of Interpolant to ((x-2)_+)^3')

절단 멱함수의 적분은 다음과 같습니다.

이 함수와 원래 함수에 대한 보간 함수를 적분한 것 간의 차이를 플로팅해도 역시 오차가 반올림 오차 범위 이내라는 점을 알 수 있습니다.

ipp = fnint(pp);
plot(xx, fnval(ipp,xx) - subplus(xx-2).^4/4)
title('Error in Integral of Interpolant to ((x-2)_+)^3')

CSAPE 명령

csapi와 마찬가지로, csape 명령을 실행하면 주어진 데이터에 대한 3차 스플라인 보간을 얻을 수 있습니다. 하지만 이 명령에서는 다양한 끝점 조건이 추가로 허용됩니다. 이 함수의 가장 간단한 형식은 다음과 같습니다.

pp = csape(x,y)

이 명령에는 라그랑주 끝점 조건이 사용되는데, 이 조건은 csapi에서 사용하는 not-a-knot 조건에 대한 일반적인 대안입니다. csape는 보간의 값을 직접 반환하지 않고 ppform만 반환합니다.

예를 들어, 다음 함수에 대한 보간을 다시 생각해 보겠습니다.

csapi가 제대로 재현하지 못하는 함수입니다. csape에서 얻은 보간의 오차(빨간색)와 함께, csapi에서 반환된 not-a-knot 보간의 오차(검은색)를 플로팅합니다.

exact = subplus(xx-1).^3;
plot(xx, fnval(csapi(x,subplus(x-1).^3),xx) - exact,'k')
hold on
plot(xx, fnval(csape(x,subplus(x-1).^3),xx) - exact,'r')
title('Error in Not-a-Knot vs. Lagrange End Conditions')
legend({'Not-a-Knot' 'Lagrange'});
hold off

이 경우 두 보간의 차이는 그다지 크지 않습니다.

기타 끝점 조건: '자연' 스플라인 보간

csape 명령을 사용하면 3차 보간 스플라인에 몇 가지 다른 유형의 끝점 조건도 지정할 수 있습니다. 예를 들어, 다음 명령은

pp = csape(x,y,'variational')

소위 '자연' 끝점 조건을 사용합니다. 이는 2계 도함수가 양 끝 절점에서 0이라는 의미입니다.

이 단계에서는 '자연' 3차 스플라인 보간을 다음 함수에 적용하고

오차를 플로팅하는 방법을 보여줍니다. 아래 코드는 'variational' 인수와 동일한 대체 인수 구문으로 '자연' 스플라인 보간을 계산하는 코드입니다. 'second'를 사용하여, csape가 양 끝 데이터 지점에서 2계 도함수를 디폴트 값인 0으로 설정하도록 지정했습니다.

pp = csape(x,subplus(x-2).^3,'second');
plot(xx, fnval(pp,xx) - subplus(xx-2).^3)
title('Error in ''Natural'' Spline Interpolation to ((x-2)_+)^3')

오른쪽 끝점 근처의 큰 오차에 주목하십시오. 이 같은 오차가 발생한 이유는 '자연' 끝점 조건에서 그곳에 0인 2계 도함수가 있도록 묵시적으로 요구했기 때문입니다.

기타 끝점 조건: 2계 도함수 규정하기

오차를 줄이기 위해 올바른 2계 도함수를 명시적으로 사용할 수도 있습니다. 우선, 끝점에서 절단 멱함수의 올바른 2계 도함수 값을 계산합니다.

endcond = 6*subplus(x([1 end])-2);

그런 다음, 끝점에 있는 2계 도함수가 방금 계산한 2계 도함수 값과 일치하도록 지정하여 보간을 만듭니다. 이를 위해서는 데이터 값과 함께 왼쪽 끝점 조건에는 endcond(1)을 제공하고 오른쪽 끝점 조건에는 endcond(2)를 제공하면 됩니다.

pp = csape(x,[endcond(1) subplus(x-2).^3 endcond(2)], 'second');
plot(xx, fnval(pp,xx) - subplus(xx-2).^3,'r')
title(['Error in Spline Interpolation to ((x-1)_+)^3'; ...
       '  When Matching the 2nd Derivative at Ends  '])

기타 끝점 조건: 기울기 규정하기

csape에서는 또한 끝점 기울기도 지정할 수 있습니다. 이것이 clamped(또는 complete) 3차 스플라인 보간입니다. 다음 명령문은

pp = csape(x,[sl,y,sr],'clamped')

왼쪽 끝 데이터 지점의 기울기는 sl이고 오른쪽 끝 데이터 지점의 기울기는 sr인 3차 스플라인 보간을 데이터(x, y)에 대해 생성합니다.

기타 끝점 조건: 혼합된 끝점 조건

이들 조건을 혼합하는 것도 가능합니다. 예를 들어, 많이 연습해 본 다음 절단 멱함수는

x=0에서 기울기가 0이고 x=6(마지막 데이터 지점)에서 2계 도함수가 30입니다.

따라서 왼쪽 끝점의 기울기와 오른쪽 끝점의 곡률에 그 값을 사용하면 결과로 생성되는 보간에 오차가 없으리라 예상됩니다.

pp = csape(x, [0 subplus(x-1).^3 30], [1 2]);
plot(xx, fnval(pp,xx) - subplus(xx-1).^3)
title(['Error in Spline Interpolation to ((x-1)_+)^3'; ...
       '        with Mixed End Conditions.          '])

기타 끝점 조건: 주기 조건

주기 끝점 조건을 규정하는 것도 가능합니다. 예를 들어, 사인 함수가 2*pi의 주기를 갖고 지점 (pi/2)*(-2:2)에서 [0 -1 0 1 0]의 값을 갖습니다. 이러한 지점에서 사인 함수와 사인 함수의 주기적 3차 스플라인 보간 사이의 차이는 2%에 불과합니다. 나쁘지 않은 수준입니다.

x = (pi/2)*(-2:2);
y = [0 -1 0 1 0];
pp = csape(x,y, 'periodic' );
xx = linspace(-pi,pi,201);
plot(xx, sin(xx) - fnval(pp,xx), 'x')
title('Error in Periodic Cubic Spline Interpolation to sin(x)')

CSAPI 또는 CSAPE에 명시적으로 포함하지 않은 끝점 조건

csapi 또는 csape에 명시적으로 포함하지 않은 끝점 조건은 csape 디폴트 끝점 조건으로 보간을 생성한 다음, 그 보간에 0의 값과 몇 가지 끝점 조건에 대한 보간의 알맞은 스칼라 배수를 더하는 방법으로 처리할 수 있습니다. 만족해야 할 `비표준' 끝점 조건이 두 개 있는 경우 2×2 선형 시스템을 먼저 풀어야 할 수 있습니다.

예를 들어, 다음 데이터에 대한 3차 스플라인 보간 s를 계산하고

x = 0:.25:3;
q = @(x) x.*(-1 + x.*(-1+x.*x/5));
y = q(x);

아래의 조건을

lambda(s) := a * (Ds)(e) + b * (D^2 s)(e) = c

점 e에서 s의 1계 도함수와 2계 도함수에 적용한다고 가정하겠습니다.

데이터는 다음과 같은 특정 파라미터로 이 끝점 조건을 만족하는 4차 다항식에서 생성되었습니다.

e = x(1);
a = 2; b = -3; c = 4;

이 특정 조건을 만족하는 보간을 생성하기 위해, 먼저 디폴트 끝점 조건을 가진 보간 함수를 생성하고,

pp1 = csape(x,y);

위 보간 함수의 첫 번째 다항식 조각의 1계 도함수를 생성합니다.

dp1 = fnder(fnbrk(pp1,1));

또한 아래와 같이 e에서 기울기가 1이라고 지정하고 0의 데이터 값에 대한 3차 스플라인 보간을 생성합니다.

pp0 = csape(x,[1,zeros(size(y)),0], [1,0]);

첫 번째 다항식 조각의 1계 도함수도 생성합니다.

dp0 = fnder(fnbrk(pp0,1));

그런 다음 아래와 같이 pp1과 pp0 모두에 대해 lambda를 계산하고,

lam1 = a*fnval(dp1,e) + b*fnval(fnder(dp1),e);
lam0 = a*fnval(dp0,e) + b*fnval(fnder(dp0),e);

3차 스플라인을 얻기 위해 pp1과 pp0의 올바른 일차 결합을 생성합니다.

s := pp1 + ((c - lambda(pp1))/lambda(pp0)) * pp0

이 3차 스플라인은 오른쪽 끝점에서 디폴트 끝점 조건뿐 아니라 원하는 조건도 만족합니다. 이 일차 결합은 fncmb의 도움을 받아 만듭니다.

s = fncmb(pp0,(c-lam1)/lam0,pp1);

보간 오차의 플롯을 통해, s가 디폴트 조건을 가진 보간 pp1보다 e 근처에서 4차 다항식을 약간 더 낫게 피팅한다는 것을 알 수 있습니다.

xx = (-.3):.05:.7; yy = q(xx);
plot(xx, fnval(pp1,xx) - yy, 'x')
hold on
plot(xx, fnval(s,xx) - yy, 'o')
hold off
legend({'Default conditions' 'Nonstandard conditions'},'location','SE')

보간의 3계 도함수에 아래 조건을 적용하려는 경우(4차 다항식이 이 조건을 만족함),

mu(s) := (D^3 s)(3) = 14.6

0의 값에 대해 보간하고 왼쪽 끝점에 0의 1계 도함수가 있어 pp0으로부터 확실히 독립적일 수 있는 3차 스플라인을 추가로 생성합니다.

pp2 = csape(x,[0,zeros(size(y)),1],[0,1]);

그런 다음, 아래의 일차 결합에서 계수 d0과 d2를 구합니다.

s := pp1 + d0*pp0 + d2*pp2

이 일차 결합은 아래의 선형 시스템을 푸는 일차 결합입니다.

lambda(s) = c

mu(s) = 14.6

pp0과 pp2가 둘 다 모든 보간 지점에서 사라지므로, d0과 d2 중 어떤 것을 선택하든 s가 주어진 데이터와 일치하게 됩니다.

재미 삼아 MATLAB® 인코딩 기능을 사용해 lambda(pp_j)와 mu(pp_j)(이때 j=0:2임)를 계산하는 루프를 작성해 보겠습니다.

dd = zeros(2,3);
for j=0:2
   J = num2str(j);
   eval(['dpp',J,'=fnder(pp',J,');']);
   eval(['ddpp',J,'=fnder(dpp',J,');']);
   eval(['dd(1,1+',J,')=a*fnval(dpp',J,',e)+b*fnval(ddpp',J,',e);']);
   eval(['dd(2,1+',J,')=fnval(fnder(ddpp',J,'),3);']);
end

pp0, pp1, pp2에 대한 lambda와 mu의 값이 주어지므로 올바른 일차 결합을 정의하는 계수를 구합니다.

d = dd(:,[1,3])\([c;14.6]-dd(:,2));
s = fncmb(fncmb(pp0,d(1),pp2,d(2)),pp1);

xxx = 0:.05:3;
yyy = q(xxx);
plot(xxx, yyy - fnval(s,xxx),'x')
title('Error in Spline Interpolant to y = x*(-1 + x*(-1+x*x/5))')

재확인을 위해, 이 함수에 대한 complete 3차 스플라인 보간에서 얻은 오차와 이 오차를 비교합니다.

hold on
plot(xxx, yyy - fnval(csape(x,[-1,y,-7+(4/5)*27],'clamped'),xxx),'o')
hold off
legend({'Nonstandard conditions' 'Endslope conditions'})

두 오차가 끝점 근처에서만 다릅니다(크게 다르지는 않음). 이는 곧 pp0과 pp2가 둘 다 각각의 끝점 근처에서만 상당한 크기라는 사실을 증명해 줍니다.

마지막 확인으로, s가 3에서 원하는 3계 도함수 조건을 만족하는지 확인합니다.

fnval(fnder(s,3),3)

ans = 14.6000