Mathematica - splajn dyskretny

SAVE OFFLINE

Dane są punkty: (-4.3, 0), (-3, 0.0111), (-1.4, 0.375), (0, 1), (1.5, 0.325), (3, 0.0111). Przybliżyć trend ich ułożenia splajnem 3 stopnia. Wczytujemy odcięte i rzędne: xa = {- 4.3, - 3, - 1.4, 0, 1.5, 3} {- 4.3, - 3, - 1.4, 0, 1.5, 3} ya = {0, .0111, 0.375, 1, 0.325, 0.0111} {0, 0.0111, 0.375, 1, 0.325, 0.0111}

Znalezienie współczynników wielomianów składowych polega na rozwiązaniu równania macierzowego A wsp = w, gdzie A - macierz kwadratowa, wsp - wektor poszukiwanych współczynników, w - wektor prawych stron. Równanie to jest podane w pliku “Interpolacja” na stronie 11 (górne - dla warunków 12a-12d) Macierz A Macierz A składa się z 9 bloków, ułożonych w trzech “superwierszach” i trzech “superkolumnach”. Policzymy po kolei każdy z nich. Numer bloku składa się z cyfry oznaczającej “superwiersz” macierzy A i cyfry oznaczającej jej “superkolumnę”. Blok 11 ma wymiar n × n, gdzie n - liczba przedziałów, na każdy z których przypada jedna funkcja składowa. Inaczej mówiąc, liczba przedziałów to numer ostatniego punktu, jeśli numerowanie punktów zaczyna się od 0. Wyrazy w tym bloku występują tylko na przekątnej i każdy z nich ma postać (xi - xi-1 )3 , i = 1 ÷ n. Ponieważ jednak Mathematica nie uznaje numeru 0, trzeba zacząć numerowanie od 1, a więc ostatni punkt będzie miał u nas numer 6, czyli ogólnie taki, ile wynosi długość wektora odciętych (lub rzędnych) - Length[xa]. Pierwszy wyraz ma więc u nas postać (x2 - x1 )3 , a ostatni (x6 - x5 )3 , czyli zakres zmienności indeksu jest od 2 do 6, tzn. od 2 do Length[xa]. Ponadto, macierz diagonalną tworzymy poleceniem DiagonalMatrix[p], gdzie p - wektor złożony z wyrazów stojących na przekątnej, zaś sam wektor p tworzymy poleceniem Table: a11 = DiagonalMatrixTable(xa[[i]] - xa[[i - 1]])3 , {i, 2, Length[xa]}; MatrixForm[a11] 2.197 0. 0. 0. 0. 0. 4.096 0. 0. 0. 0. 0. 2.744 0. 0. 0. 0. 0. 3.375 0. 0. 0. 0. 0. 3.375

Jest pięć przedziałów i blok 11 ma wymiar 5 × 5. Blok 12 ma wymiar n × (n - 1). Jest o jedną kolumnę mniej - porównując do pełnej macierzy diagonalnej, nie ma pierwszej kolumny. Jeśli wziąć pełną macierz diagonalną, to wyrazy w tym bloku występują również tylko na przekątnej i każdy z nich ma postać (xi - xi-1 )2 , i = 1 ÷ n. Pierwszy wyraz przekątnej pełnej macierzy ma u nas postać (x2 - x1 )2 , a ostatni (x6 - x5 )2 , czyli zakres zmienności indeksu jest taki sam, jak w bloku 11, tzn. od 2 do Length[xa]. Pełna macierz ma wymiar 5 × 5, czyli ogólnie (Length[xa] - 1) × (Length[xa] - 1). Usuwamy z takiej macierzy pierwszą kolumnę, czyli wypisujemy kolumny od drugiej do ostatniej, tzn. do mającej numer Length[xa] - 1:

2

splajn.nb

a12 = DiagonalMatrixTable(xa[[i]] - xa[[i - 1]])2 , {i, 2, Length[xa]}[[ All, 2 ;; Length[xa] - 1]]; MatrixForm[a12] 0. 0. 0. 0. 2.56 0. 0. 0. 0. 1.96 0. 0. 0. 0. 2.25 0. 0. 0. 0. 2.25

Blok 13 jest prawie identyczny z blokiem 11, tyle że wyrazy w tym bloku mają postać xi - xi-1 , i = 1 ÷ n: a13 = DiagonalMatrix[Table[xa[[i]] - xa[[i - 1]], {i, 2, Length[xa]}]]; MatrixForm[a13] 1.3 0. 0. 0. 0. 0. 1.6 0. 0. 0. 0. 0. 1.4 0. 0. 0. 0. 0. 1.5 0. 0. 0. 0. 0. 1.5

Blok 21 ma wymiar (n - 1) × n. Jest o jeden wiersz mniej - porównując do pełnej macierzy diagonalnej, nie ma ostatniego. Jeśli wziąć pełną macierz diagonalną, to jest ona prawie taka sama jak dla bloku 12, z tym że wyrazy mają postać 3 (xi - xi-1 )2 , i = 1 ÷ n. Macierz tę tworzymy identycznie jak dla bloku 12, a następnie usuwamy z niej ostatni wiersz, czyli wypisujemy wiersze od pierwszego do przedostatniego, tzn. do mającego numer Length[xa] - 2 (bo ogólnie jest Length[xa]-1 wierszy): a21 = DiagonalMatrixTable3 (xa[[i]] - xa[[i - 1]])2 , {i, 2, Length[xa]}[[ 1 ;; Length[xa] - 2]]; MatrixForm[a21] 5.07 0. 0. 0. 0. 7.68 0. 0. 0. 0. 5.88 0. 0. 0. 0. 6.75

0. 0. 0. 0.

Blok 22 ma wymiar (n - 1) × (n - 1) - porównując do pełnej macierzy diagonalnej, nie ma ostatniego wiersza i pierwszej kolumny. Jeśli wziąć pełną macierz diagonalną, to jest ona prawie taka sama jak dla bloku 13, z tym że wyrazy mają postać 2 (xi - xi-1 ), i = 1 ÷ n. Macierz tę tworzymy identycznie jak dla bloku 13, następnie usuwamy z niej ostatni wiersz, czyli wypisujemy wiersze od pierwszego do przedostatniego, tzn. do mającego numer Length[xa] - 2, a na końcu usuwamy pierwszą kolumnę, czyli wypisujemy kolumny od drugiej do ostatniej, tzn. do mającej numer Length[xa] - 1: a22 = DiagonalMatrix[Table[2 (xa[[i]] - xa[[i - 1]]), {i, 2, Length[xa]}]][[ All, 2 ;; Length[xa] - 1]][[1 ;; Length[xa] - 2]]; MatrixForm[a22] 0. 0. 0. 3.2 0. 0. 0. 2.8 0. 0. 0. 3.

0. 0. 0. 0.

Blok 23 jest nieco bardziej złożony. Można zauważyć, że składa się on z dwóch macierzy jednostkowych pokazanych poniżej, z których wycięto po jednym wierszu: ostatni dla macierzy dodatniej, pierwszy dla macierzy ujemnej, czyli zostawiono wyrazy oznaczone na czerwono. Dodanie do siebie obu takich “okrojonych” macierzy daje właśnie blok 23. 1 0 0 0 0

0 1 0 0 0

0 0 1 0 0

0 0 0 1 0

0 0 0 , 0 1

-1 0 0 0 0 0 -1 0 0 0 0 0 -1 0 0 , 0 0 0 -1 0 0 0 0 0 -1

splajn.nb

1 0 0 0

0 1 0 0

0 0 1 0

0 0 0 1

0 0 0 0

+

0 -1 0 0 0 0 0 -1 0 0 0 0 0 -1 0 0 0 0 0 -1

=

3

1 -1 0 0 0 0 1 -1 0 0 0 0 1 -1 0 0 0 0 1 -1

Blok ten ma wymiar (n - 1) × n, czyli w naszym przypadku 4 × 5, a ogólnie (Length[xa] - 1) × (Length[xa] - 2). Macierze pierwotne mają więc wymiar (Length[xa] - 1) × (Length[xa] - 1) i ponieważ są to macierze jednostkowe, tworzymy je poleceniem IdentityMatrix[p], gdzie p = Length[xa] - 1. Z pierwszej macierzy wypisujemy wiersze od pierwszego do przedostatniego, tzn. do mającego numer Length[xa] - 2, a z drugiej macierzy - od drugiego do ostatniego, tzn. do mającego numer Length[xa] - 1: a23 = IdentityMatrix[Length[xa] - 1][[1 ;; Length[xa] - 2]] IdentityMatrix[Length[xa] - 1][[2 ;; Length[xa] - 1]]; MatrixForm[a23] 1 -1 0 0 0 0 1 -1 0 0 0 0 1 -1 0 0 0 0 1 -1

Blok 31 jest prawie identyczny z blokiem 13, tyle że wyrazy w tym bloku mają postać 6 (xi - xi-1 ), i = 1 ÷ n: a31 = DiagonalMatrix[Table[6 (xa[[i]] - xa[[i - 1]]), {i, 2, Length[xa]}]]; MatrixForm[a31] 7.8 0. 0. 0. 0. 0. 9.6 0. 0. 0. 0. 0. 8.4 0. 0. 0. 0. 0. 9. 0. 0. 0. 0. 0. 9.

Blok 32 ma budowę analogiczną do bloku 23, z tym że macierzy jednostkowych pokazanych poniżej wycina się po jednej kolumnie: pierwszą dla macierzy dodatniej, ostatnią dla macierzy ujemnej. Dodanie do siebie obu takich “okrojonych” macierzy pomnożonych przez 2 daje właśnie blok 32.

2

1 0 0 0 0

0 1 0 0 0

0 0 1 0 0

0 0 0 1 0

2

0 1 0 0 0

0 0 1 0 0

0 0 0 1 0

0 0 0 0 1

0 0 0 , 2 0 1

+ 2

-1 0 0 0 0 0 -1 0 0 0 0 0 -1 0 0 , 0 0 0 -1 0 0 0 0 0 -1

-1 0 0 0 0 -1 0 0 0 0 -1 0 0 0 0 -1 0 0 0 0

=

-2 0 0 0 2 -2 0 0 0 2 -2 0 0 0 2 -2 0 0 0 2

Blok ten ma wymiar n × (n - 1), czyli w naszym przypadku 5 × 4, a ogólnie (Length[xa] - 2) × (Length[xa] - 1). Macierze pierwotne mają wymiar (Length[xa] - 1) × (Length[xa] - 1). Z pierwszej macierzy wypisujemy kolumny od drugiej do ostatniej, tzn. do mającej numer Length[xa] - 1, a z drugiej macierzy - od pierwszej do przedostatniej, tzn. do mającej numer Length[xa] - 2:

4

splajn.nb

a32 = 2 IdentityMatrix[Length[xa] - 1][[All, 2 ;; Length[xa] - 1]] 2 IdentityMatrix[Length[xa] - 1][[All, 1 ;; Length[xa] - 2]]; MatrixForm[a32] -2 0 0 0 2 -2 0 0 0 2 -2 0 0 0 2 -2 0 0 0 2

Blok 33 to tablica zer o wymiarze takim jak liczba przedziałów, czyli (Length[xa] - 1) × (Length[xa] 1): a33 = Array[0 &, {Length[xa] - 1, Length[xa] - 1}]; MatrixForm[a33] 0 0 0 0 0

0 0 0 0 0

0 0 0 0 0

0 0 0 0 0

0 0 0 0 0

Wszystkie bloki składamy w całość poleceniem “ArrayFlatten[{superwiersz_1, superwiersz_2, ...}]”: a = ArrayFlatten[{{a11, a12, a13}, {a21, a22, a23}, {a31, a32, a33}}]; a // MatrixForm 2.197 0. 0. 0. 0. 0. 0. 0. 0. 1.3 0. 0. 0. 0. 0. 4.096 0. 0. 0. 2.56 0. 0. 0. 0. 1.6 0. 0. 0. 0. 0. 2.744 0. 0. 0. 1.96 0. 0. 0. 0. 1.4 0. 0. 0. 0. 0. 3.375 0. 0. 0. 2.25 0. 0. 0. 0. 1.5 0. 0. 0. 0. 0. 3.375 0. 0. 0. 2.25 0. 0. 0. 0. 1.5 5.07 0. 0. 0. 0. 0. 0. 0. 0. 1 -1 0 0 0 0. 7.68 0. 0. 0. 3.2 0. 0. 0. 0 1 -1 0 0 0. 0. 5.88 0. 0. 0. 2.8 0. 0. 0 0 1 -1 0 0. 0. 0. 6.75 0. 0. 0. 3. 0. 0 0 0 1 -1 7.8 0. 0. 0. 0. -2 0 0 0 0 0 0 0 0 0. 9.6 0. 0. 0. 2 -2 0 0 0 0 0 0 0 0. 0. 8.4 0. 0. 0 2 -2 0 0 0 0 0 0 0. 0. 0. 9. 0. 0 0 2 -2 0 0 0 0 0 0. 0. 0. 0. 9. 0 0 0 2 0 0 0 0 0

Wektor prawej strony w Składa się on z tylu wyrazów, ile wierszy ma macierz A; bloki pierwszego superwiersza mają po n wierszy, drugiego - po n - 1 wierszy, trzeciego - po n wierszy, więc razem jest 3n - 1, czyli 3(Length[xa] - 1) - 1 wierszy. Na początku stworzymy tablicę złożoną z tylu właśnie zer: w = Array[0 &, 3 (Length[xa] - 1) - 1] {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}

Pierwsze n wyrazów (czyli od pierwszego do mającego numer Length[xa] - 1) należy wypełnić różnicami postaci yi - yi-1 , i = 1 ÷ n. Można to zrobić komendą “Table”.

splajn.nb

5

w[[1 ;; Length[xa] - 1]] = Table[ya[[i + 1]] - ya[[i]], {i, 1, Length[xa] - 1}]; MatrixForm[w] 0.0111 0.3639 0.625 - 0.675 - 0.3139 0 0 0 0 0 0 0 0 0

Współczynniki wielomianów składowych wsp = A-1 w, czyli wsp = Inverse[a].w; MatrixForm[wsp] 0.0124767 0.0386224 - 0.196059 0.190471 - 0.0594927 0.048659 0.234047 - 0.589401 0.267717 - 0.0125471 0.0507096 0.503039 0.00554266 - 0.476984

Wielomiany składowe Każdy wielomian składowy ma postać: si = ai (x - xi-1 )3 + bi (x - xi-1 )2 + ci (x - xi-1 ) + di , gdzie i = 1 ÷ n, zaś a, b, c, d - współczynniki, będące m.in. odpowiednimi wyrazami wektora wsp, ale nie tylko. W wektorze wsp zawarte są następujące współczynniki: 0.0124767 0.0386224 -0.196059 0.190471 -0.0594927 0.048659 0.234047 -0.589401 0.267717 -0.0125471 0.0507096 0.503039 0.00554266 -0.476984

a1 a2 a3 a4 a5 b2 b3 b4 b5 c1 c2 c3 c4 c5

6

splajn.nb

Z zasady b1 = 0, zaś di = yi , i = 1 ÷ n. Najpierw tworzymy tablicę n zer, czyli tylu, ile jest wielomianów składowych (przedziałów): ff = Table[0, {i, 1, Length[xa] - 1}] {0, 0, 0, 0, 0}

Pierwszy wielomian jest nieco inny niż pozostałe, bo b1 = 0, dlatego stworzymy go oddzielnie. Niezależnie od tego, z którym przedziałem (wielomianem) mamy do czynienia, można zauważyć, że postać si jest iloczynem skalarnym dwóch wektorów: pierwszy to [ai , bi , ci , di ], drugi zaś to [(x - xi-1 )3 , (x - xi-1 )2 , (x - xi-1 )1 , (x - xi-1 )0 ]. Drugi jest nieco łatwiej stworzyć - jest to tablica potęg o malejącym wykładniku i niezmiennej podstawie: Table(x - xa[[1]])3-k , {k, 0, 3}

(nie wykonywać!) Zapis 3 - k w wykładniku jest wymuszony tym, że znacznik wyrazów tablicy (tu - k) musi rosnąć. Dla pierwszego wektora trzeba wybrać odpowiednie wyrazy z wektora wsp: - a1 to wyraz pierwszy, czyli wsp[[1]], - b1 = 0, - c1 to wyraz odpowiadający pierwszemu wierszowi trzeciego “superwiersza”, a więc ma numer n (liczba wyrazów pierwszego “superwiersza”) + n - 1 (liczba wyrazów drugiego “superwiersza”) + 1 (pierwszy wiersz trzeciego “superwiersza”) = 2n = 2(Length[xa] - 1), - d1 = y1 , czyli ya[[1]]. Cały wektor jest więc taki: [wsp[[1]], 0, wsp[[2 (Length[xa] - 1)]], ya[[1]]]

a pierwszy wielomian składowy - taki: ff[[1]] = Table(x - xa[[1]])3-k , {k, 0, 3}. {wsp[[1]], 0, wsp[[2 (Length[xa] - 1)]], ya[[1]]} // Chop - 0.0125471 (4.3 + x) + 0.0124767 (4.3 + x)3

Komenda "Chop" ucina wszelkie niepotrzebne cyfry, np. bardzo niskie potęgi 10 zastępuje zerem, upraszcza zapis typu 1.x do postaci x, czy też pomija wyrazy typu 0.x. Następne wielomiany tworzymy już poleceniem “Do”. Wektor Table(x - xa[[1]])3-k , {k, 0, 3} pozostaje bez zmian. Dla wektora grupującego współczynniki trzeba wybrać odpowiednie wyrazy z wektora wsp: - ai to wyraz o takim numerze, jaki ma dany wielomian, czyli jeśli jest to numer i, to wsp[[i]], - bi to wyraz wybierany następująco: dla drugiego wielomianu jest to wyraz o numerze odpowiadającym numerowi pierwszego wiersza drugiego “superwiersza”, dla trzeciego wielomianu - o numerze odpowiadającym numerowi drugiego wiersza drugiego “superwiersza”, itd., ogólnie dla i-tego wielomianu jest to wyraz o numerze odpowiadającym numerowi (i-1)-szego wiersza drugiego “superwiersza”, a więc o numerze n (liczba wyrazów pierwszego “superwiersza”) + (i - 1) = Length[xa] - 1 + i - 1 = i + Length[xa] - 2, - ci to wyraz wybierany analogicznie jak bi , z tym że wybór dotyczy trzeciego “superwiersza”: ogólnie dla i-tego wielomianu jest to wyraz o numerze odpowiadającym numerowi i-tego wiersza trzeciego “superwiersza”, a więc o numerze n (liczba wyrazów pierwszego “superwiersza”) + n - 1 (liczba wyrazów drugiego “superwiersza”) + i = i + 2(Length[xa] - 1) - 1, - di = yi , czyli ya[[i]]. Doff[[i]] = Table(x - xa[[i]])3-k , {k, 0, 3}.{wsp[[i]], wsp[[i + Length[xa] - 2]], wsp[[i + 2 (Length[xa] - 1) - 1]], ya[[i]]}, {i, 2, Length[xa] - 1}

splajn.nb

7

Wielomiany składowe: ff // MatrixForm - 0.0125471 (4.3 + x) + 0.0124767 (4.3 + x)3 0.0111 + 0.0507096 (3 + x) + 0.048659 (3 + x)2 + 0.0386224 (3 + x)3 0.375 + 0.503039 (1.4 + x) + 0.234047 (1.4 + x)2 - 0.196059 (1.4 + x)3 1 + 0.00554266 x - 0.589401 x2 + 0.190471 x3 0.325 - 0.476984 (- 1.5 + x) + 0.267717 (- 1.5 + x)2 - 0.0594927 (- 1.5 + x)3

Połączenie w splajn: Każdy z wielomianów składowych obowiązuje tylko w swoim przedziale, a i-ty przedział zawiera się między xi-1 a xi . Nasuwa się pomysł taki, jak przy rysowaniu funkcji złożonych z kawałków o różnych dziedzinach (z zagnieżdżanymi komendami “If”), ale wyklucza to automatyzację. Można jednak zastosować polecenie “Piecewise”, które ma następującą składnię: Piecewise[{{wzór_ 1, warunek_ 1}, {wzór_ 2}, {warunek_ 2}, ...}]

gdzie “wzór” określa w naszym przypadku wielomian składowy, a “warunek” to u nas nakaz, by argument tego wielomianu nie przekraczał górnej granicy przypisanego mu przedziału. Uwaga: Nie wprowadzałem tej komendy wcześniej, bo we wcześniejszych wersjach Mathematiki jej nie ma (jest od 9 lub 10 wzwyż) i chciałem pokazać też bardziej ogólną możliwość rysowania funkcji kawałkami ciągłych. Pojedynczy element komendy “Piecewise” będzie miał postać - np. dla wielomianu 3: 0.375 + 0.503039 (1.4 + x) + 0.234047 (1.4 + x)2 - 0.196059 (1.4 + x)3 , x ≤ 0

Ogólnie jest tak, że i-temu wielomianowi jest przypisany warunek x ≤ xi+1 (dla powyższego przykładu 0 jest czwartym wyrazem wektora xa). Wszystkie takie elementy można opisać komendą “Table”, w której znacznik zmienia się od 1 do liczby przedziałów, czyli od 1 do Length[xa] - 1: Table[{ff[[i]], x