Interpolacja splajnami kubicznymi

Przybliżanie, metoda najmniejszych kwadratów, wielomiany interpolacyjne i inne.
janusz47
Użytkownik
Użytkownik
Posty: 7910
Rejestracja: 18 mar 2009, o 16:24
Płeć: Mężczyzna
Podziękował: 30 razy
Pomógł: 1670 razy

Interpolacja splajnami kubicznymi

Post autor: janusz47 »

Interpolacja splajnami kubicznymi


Splajn kubiczny to najbardziej popularny rodzaj splajnu, W matematyce obliczeniowej jest wygodnym narzędziem do

gładkiej interpolacji i numerycznego różniczkowania.

Dla każdego przedziału \(\displaystyle{ [x_{i}, x_{i+1}], \ \ i = 1,2,...n }\) konstruujemy splajn postaci:

\(\displaystyle{ S_{i}(x) = a_{i} + b_{i}\cdot (x-x_{i}) +c_{i}\cdot (x-x_{i})^2 +d_{i}\cdot (x-x_{i})^3 = S(x), \ \ i=1,2,..., n-1 \ \ (1)}\)

gdzie:

\(\displaystyle{ a_{i}, b_{i}, c_{i}, d_{i} }\) są liczbami, które musimy znaleźć.

Niech \(\displaystyle{ h_{i} = x_{i+1} -x_{i}.}\)

W każdym punkcie węzłowym musi zachodzić warunek dopasowania:

\(\displaystyle{ f(x_{i}) = S(x_{i}), \ \ i = 1,2, ..., n.}\)

Stąd

\(\displaystyle{ a_{i}= f(x_{i}), \ \ i=1,2,...,n. }\)

Ponadto musi być spełniony warunek ciągłości:

\(\displaystyle{ S_{i+1}(x_{i+1}) = S_{i}(x_{i+1}),}\)

\(\displaystyle{ f(x_{i+1}) = a_{i} +b_{i} (x_{i+1}-x_{i}) +c_{i}(x_{i+1}-x_{i})^2 + d_{i} (x_{i+1}-x{i}). }\)

\(\displaystyle{ a_{i+1} = a_{i} + b h_{i} + c_{i} h^2_{i} +d_{i} h^3_{i}, \ \ i=1,2,..., n-1 \ \ (2)}\)

Różniczkując równanie \(\displaystyle{ (1) }\), otrzymujemy

\(\displaystyle{ S^{'}(x) = b_{i} + 2c_{i} \cdot (x-x_{i}) +3d_{i}\cdot (x-x_{i})^2,}\)

\(\displaystyle{ S^{''}(x) = 2c_{i} + 6d_{i}\cdot (x-x_{i}). }\)

Z warunków ciągłości pierwszej i drugiej pochodnej otrzymujemy odpowiednio równania:

\(\displaystyle{ S^{'}_{i}(x_{i}) = S^{'}_{i-1}(x_{i}) }\)

\(\displaystyle{ b_{i} = b_{i-1} +2c_{i-1}\cdot (x_{i} -x_{i-1}) +3d_{i-1}\cdot (x_{i} - x_{i-1})^2= }\)

\(\displaystyle{ = b_{i-1} +2c_{i-1}\cdot h_{i-1} +3d_{i-1}\cdot h^2_{i-1}, \ \ i = 2, 3, ..., n \ \ (3) }\)

\(\displaystyle{ S^{''}_{i}(x_{i}) = S^{''}_{i-1}(x_{i}) }\)

\(\displaystyle{ 2c_{i} =2c_{i-1} +6d_{i-1}\cdot (x_{i}-x_{i-1}), }\)

\(\displaystyle{ c_{i} = c_{i-1} +3d_{i-1}\cdot h_{i-1} \ \ i =2, 3, ..., n \ \ (4) }\)

Zwiększając index \(\displaystyle{ (4) }\) o jeden - równanie \(\displaystyle{ (4) }\) możemy zapisać w postaci

\(\displaystyle{ c_{i+1} = c_{i} +3d_{i}\cdot h_{i}, \ \ i =1,2, ..., n-1 \ \ (5) }\)

Wyznaczamy z równania \(\displaystyle{ (5) \ \ d_{i} }\) i podstawiamy do \(\displaystyle{ (2) }\),

otrzymując po zmianie kolejności

\(\displaystyle{ b_{i} = \frac{1}{h_{i-1}}\cdot (a_{i+1}-a_{i}) - \frac{h_{1}}{3}\cdot (c_{i+1}+2c_{i}), \ \ i=1,2,..., n-1 \ \ (5^{'}) }\)

lub

\(\displaystyle{ b_{i} = \frac{1}{h_{i}}\cdot (a_{i+1}-a_{i}) - \frac{h_{i}}{3}\cdot (c_{i+1}+2c_{i}), \ \ i=1,2,..., n-1 \ \ (6)}\)

\(\displaystyle{ b_{i-1} =\frac{1}{h_{i}}\cdot (a_{i}-a_{i-1}) - \frac{h_{i-1}}{3}\cdot (c_{i}+2c_{i-1}), \ \ i=1,2,..., n-1 }\)

Podobnie, wyznaczając z równania \(\displaystyle{ (4) \ \ d_{i-1} }\) i podstawiając do równania \(\displaystyle{ (3) }\) dostajemy

\(\displaystyle{ b_{i} = b_{i-1} + h_{i-1}\cdot (c_{i} + c_{i-1}), \ \ i=2,3,..., n-1 \ \ (7) }\)

Ostatecznie na podstawie równań \(\displaystyle{ (5), (5'), (6), (7) }\) otrzymujemy

\(\displaystyle{ h_{i-1}\cdot c_{i-1} +u_{i}\cdot c_{i} +h_{i}\cdot c_{i+1} = v_{i}, \ \ i=2,3,..., n-1 \ \ (8) }\)

gdzie

\(\displaystyle{ u_{i} = 2(h_{i-1} +h_{i}), \ \ v_{i} = 3w_{i} - 3w_{i-1} , \ \ w_{i} = \frac{1}{h_{i}}\cdot (a_{i+1} -a_{i}).}\)

Równanie \(\displaystyle{ (8) }\) jest równaniem liniowym z \(\displaystyle{ n }\) niewiadomymi i \(\displaystyle{ n-2 }\) równaniami.

Warunek konieczny aby liczba równań była równa liczbie niewiadomych uzyskujemy z naturalnych warunków brzegowych:

\(\displaystyle{ S^{''}(x_{1} = S^{''}(x_{n}) \ \ (9) }\)

stąd splajny kubiczne nazywamy naturalnymi splajnami kubicznymi.

Na podstawie \(\displaystyle{ (9) }\) otrzymujemy dwa dodatkowe równania:

\(\displaystyle{ c_{1} = 0 }\) i \(\displaystyle{ c_{n} = 0 \ \ (10) }\)

Równania \(\displaystyle{ (8) }\) tworzą układ trójdiagonalny równań liniowych

\(\displaystyle{ A\cdot \vec{x} = \vec{b}, }\)

\(\displaystyle{ A = \left [\begin{matrix} 1 & 0 \\ h_{1} & 2(h_{1}+h_{2}) & h_{2} & 0 \\ & & \ddots & &\\ & 0 & h_{n-2} & 2(h_{n-2}+h_{n-1})& h_{n-1}\\ & & & 0 &1 \end{matrix}\right], }\)

\(\displaystyle{ \vec{x} = \left[ \begin{matrix} c_{1}\\c_{2}\\ \vdots \\ c_{n-1}\\ c_{n}\end{matrix} \right],}\)

\(\displaystyle{ \vec{b} =\left[\begin{matrix} 0 & \\ 3(a_{3}-a_{2})/ h_{2} -3(a_{2}-a_{1})/h_{1} & \\ \vdots & \\ 3(a_{n} - a_{n-1})/h_{n-1}-3(a_{n-1}-a_{n-2})/h_{n-2} \\ 0 & \end{matrix} \right] }\)

Otrzymaliśmy układ ściśle diagonalnie dominujący. Możemy zastosować metodę eliminacji do jego rozwiązanania, otrzymując wartości

\(\displaystyle{ (c_{i}),\ \ i=1,2,.., n.}\)

Pozostałe warości współczynników dla splajnów:

\(\displaystyle{ S_{i} (x) = a_{i} +b_{i}\cdot (x-x_{i}) + c_{i}\cdot (x-x_{i})^2 + d_{i}\cdot (x=x_{i})^3, \ \ i=1,2,...,n-1 }\)

obliczamy z układu równań:

\(\displaystyle{ a_{i}= f(x_{i}), }\)

\(\displaystyle{ b_{i} = \frac{1}{h_{i}} \cdot (a_{i+1}-a_{i}) - \frac{h_{i}}{3}\cdot (c_{i+1}+2c_{i}), }\)

\(\displaystyle{ d_{i} = \frac{c_{i+1}-c_{i}}{3h_{i}},}\)

dla \(\displaystyle{ i = 1,2,..., n-1, n. }\)

Program w OCTAVE

Kod: Zaznacz cały

function spl3(x,y,m)
% splajn kubiczny
n=length(x);
for i=1:n-1
    h(i)= x(i+1)-x(i);
    b(i)=(y(i+1)-y(i))/h(i);
end
u(2)=2*(h(1)+h(2));
v(2)=6*(b(2)-b(1));
for i=3:n-1
    u(i)=2*(h(i-1)+h(i))-h(i-1)^2/u(i-1);
    v(i)=6*(b(i)-b(i-1))-h(i-1)*v(i-1)/u(i-1);
end
z(n)=0;
for i=n-1:-1:2
    z(i)=(v(i)-h(i)*z(i+1))/u(i);
end
z(1)=0;
disp('          splajn kubiczny')
disp('___________________________________________')
disp('       x                  S(x)           ')
disp('___________________________________________')
for i=1:n-1
    for j=1:m+1
        r=(x(i+1)-x(i))/(m+1);
        t=x(i)+(j-1)*r;
        dis=(j-1)*r;
        hh=x(i+1)-x(i);
        bb=(y(i+1)-y(i))/hh-hh*(z(i+1)+2*z(i))/6;
        q=0.5*z(i)+dis*(z(i+1)-z(i))/(6*hh);
        q=bb+dis*q;
        spl3=y(i)+dis*q;
        fprintf('%12.5f %17.5f \n',t,spl3)
    end
end
fprintf('%12.5f %17.5f \n',x(n),y(n))
fprintf('\n')
Przykład

Kod: Zaznacz cały

                                                                                                                                                                                             
>> x=[0:0.2:1]';
>> y=x.^3+2;
>>  spl3(x,y,2)
          splajn kubiczny
______________________________
         x              S(x)
______________________________
     0.00000           2.00000
     0.06667           2.00035
     0.13333           2.00244
     0.20000           2.00800
     0.26667           2.01881
     0.33333           2.03681
     0.40000           2.06400
     0.46667           2.10220
     0.53333           2.15254
     0.60000           2.21600
     0.66667           2.29418
     0.73333           2.39125
     0.80000           2.51200
     0.86667           2.65885
     0.93333           2.82468
     1.00000           3.00000
ODPOWIEDZ