Interpolacja splajnami kubicznymi
: 27 lip 2022, o 19:54
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
Przykład
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')
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
