Interpolacja splajnami kubicznymi

Przybliżanie, metoda najmniejszych kwadratów, wielomiany interpolacyjne i inne.
janusz47
Użytkownik
Użytkownik
Posty: 7923
Rejestracja: 18 mar 2009, o 16:24
Płeć: Mężczyzna
Podziękował: 30 razy
Pomógł: 1673 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
Fibik
Użytkownik
Użytkownik
Posty: 980
Rejestracja: 27 wrz 2005, o 22:56
Płeć: Mężczyzna
Lokalizacja: Wrocław
Podziękował: 12 razy
Pomógł: 75 razy

Re: Interpolacja splajnami kubicznymi

Post autor: Fibik »

Interesuje mnie operacja odwrotna,
czyli mając krzywą beziera rzędu 3, znaczy te 4-ry punkty:

\(\displaystyle{ B = [P_1,P_2,P_3,P_4]}\)

i potrzebuję to przekształcić na serię odcinków,
przy warunku minimalizacji liczby tych linii dla zadanej dokładności,
np. odchyłka (od tej krzywej beziera) nie może przekraczać 0.2 [mm].

Jest na to gotowa metoda, czy też znowu muszę to samodzielnie załatwić? :)
ODPOWIEDZ