Analiza numeryczna

Rozwinięcie w szereg Taylora dla x_n=2 do drugiej pochodnej: \[f(x_{n+1})=f(x_n+h)=f(x_n)+\frac{x_{n+1}-x_n}{1!}f'(x_n)+\frac{(x_{n+1}-x_n)^2}{2!}f''(x_n)\] \[f(x_{n+1})= 3 + (x_{n+1}-2)4 + \frac{(x_{n+1}-2)^2}{2}2\] \[f(x_{n+1})= 3 + 4x_{n+1}-8 + (x_{n+1}-2)^2\] \[f(x_{n+1})= 4x_{n+1}-5 + x_{n+1}^2-4x_{n+1}+4\] \[f(x_{n+1})= x_{n+1}^2-1\]

Każdą funkcję możemy rozwinąć w punkcie w szereg Taylora z dowolną precyzją, ale naszym celem jest stworzenie metody iteracyjnej – dlatego rozwijamy tylko do pierwszej pochodnej. \[f(x_{n+1})=f(x_n+h)=f(x_n)+\frac{x_{n+1}-x_n}{1!}f'(x_n)\] \[h=-\frac{f(x_n)}{f'(x_n)}=-\frac{3}{4}\]

using Plots
f(x)=x^2-1
g(x)=3+(x-2)*4 # rozwinięcie w szerg Taylora w punkcie 2
v1 = [f(i) for i in 0:0.01:3]
v2 = [g(i) for i in 0:0.01:3]
plot(0:0.01:3,v1)
plot!(0:0.01:3,v2)
vline!([1,1.25,2])
hline!([0])
savefig("newton_steps.png")

using Plots
iter = 0.975:0.001:1.035
wielomian(x::Float64) = x^4-4x^3+6x^2-4x+1
wielomian(x::Float32) = x^4-4x^3+6x^2-4x+1

sum(wielomian.(Float64.(iter)) .== 0) # 1 pierwiastek
sum(wielomian.(Float32.(iter)) .== 0) # 21 pierwiastków

plot(wielomian.(Float64.(iter));seriestype=:scatter)
plot!(wielomian.(Float32.(iter));seriestype=:scatter)
savefig("precision_float32.png")

Zmniejszenie precyzji liczby zmiennoprzecinkowej do Float32 spowodowało, że zamiast jednego pierwiastka pojawiły się 21 pierwiastki w naszym równaniu.

iter = 0.999999:0.0000001:1.000001
sum(wielomian.(Float64.(iter)) .== 0) # 11 pierwiastków

plot(wielomian.(Float64.(iter));seriestype=:scatter)
savefig("precision_float64.png")

1. Obliczanie wartości z wyrażenia jest bezpieczniejsze z punktu widzenia odporności na błędy, jeśli jest ona ewaluowana w trakcie działania algorytmu poprzez dodawanie do niej wartości korygującej, aniżeli obliczana bez niej np. wyrażenie \(c=\frac{a+b}{2}\) należałoby zmodyfikować do równoważnej postaci \(c=a+\frac{b-a}{2}\)
2. Operacje bardziej obciążające procesor i pamięć zastępujemy lżejszymi np. w metodzie równego podziału zamiast sprawdzać znak funkcji za pomocą mnożenia f(a)f(b)<0, to należałoby porównać znak dwóch wartości signbit(f(a)) == signbit(f(b)).
3. Działanie algorytmu kończy się, gdy:
   • przekroczymy dozwoloną ilość iteracji. Program wtedy powinien zwrócić informacje o błędzie zbieżności.
   • błąd pomiędzy kolejnymi dwiema wartościami zmiennej zależnej jest mniejszy od założonego błędu. Program zwraca informację o znalezieniu pierwiastka.

\[x_1 - x_2 - x_1^3 = 0\] \[x_1 + 6x_2 - x_2^3 = 0\]

\[F(X)=\begin{pmatrix}x_1 - x_2 - x_1^3\\ x_1 + 6x_2 - x_2^3\end{pmatrix}\]

\[J(X)=\begin{pmatrix}1-3x_1^2 & -1\\ 1 & 6-3x_2^2\end{pmatrix}\]

dla x₁=0 i x₂=-3

X0=[0;-3]
f₁(x1,x2)=x1-x2-x1^3
f₂(x1,x2)=x1+6*x2-x2^3
F(x1,x2)=[f₁(x1,x2);f₂(x1,x2)]
J(x1,x2)=[1-3*x1^2 -1;6 1-3*x2^2]
H = -inv(J(X0[1],X0[2]))*F(X0[1],X0[2]) 

[-3.45 -3.45]
[-2.286700925481873 -2.761577628837895]
[-1.4624857977508279 -2.4347895758642046]
[-0.7354570543294794 -2.2724461590331337]
[1.3062224730230296 -1.6089864310800697]
[1.5668617931211841 -1.9959592237888915]
[1.5582498571428895 -2.225055047084357]
[1.5665384578179202 -2.277492142300413]
[1.5693816234671674 -2.2959024556526546]
[1.5704088548735522 -2.3025033225237737]
[1.5707814198957448 -2.304892168298706]
[1.57091679006564 -2.3057594706206888]
"Pierwiastek: [1.5707814198957448, -2.304892168298706]'"

dla x₁=-1 i x₂=0.25

[-1.0642857142857143 0.1285714285714286]
[-0.9882413732658417 -0.041127807045914566]
[-0.68995037923022 -0.5987648465831944]
[-0.012635611298492222 -0.6514670268233802]
[0.6543036057064445 0.6539861752005263]
[-0.06565776713529803 0.5789011367729759]
[-0.6023343695832487 -0.5951105769940852]
[0.05962886349665608 -0.44233426993587777]
[0.4186137245525588 0.4145724869016827]
[-0.031039129251535957 0.13199212793746262]
[-0.11802725780008667 -0.11774593357068736]
[0.0007241367500920581 -0.002594467618342408]
[0.0023704439439994496 0.0023704409744387143]
[-5.489996130446373e-9 2.1149166751655707e-8]
[-1.9027973487215778e-8 -1.902797348721577e-8]
"Pierwiastek: [-5.489996130446373e-9, 2.1149166751655707e-8]'"

Wśród systemów liniowych możemy wyróżnić takie, których znalezienie rozwiązania \(x=A^{-1}b\) jest proste.

Macierz współczynników ma postać diagonalną: \[\begin{pmatrix}a_{11} & 0 & ... & 0\\ 0 & a_{22} & ... & 0 \\ ... & ... & ... & ...\\0 & 0 & ... & a_{nn}\end{pmatrix} \begin{pmatrix}x_1\\x_2\\...\\x_n\end{pmatrix}=\begin{pmatrix}b_1\\b_2\\...\\b_n\end{pmatrix}\]

To wtedy rozwiązaniem jest wektor: \[\begin{pmatrix}x_1\\x_2\\...\\x_n\end{pmatrix}=\begin{pmatrix}\frac{b_1}{a_{11}}\\\frac{b_2}{a_{22}}\\...\\\frac{b_n}{a_{nn}}\end{pmatrix}\]

Macierz współczynników ma postać dolno-trójkątną \[\begin{pmatrix}a_{11} & 0 & ... & 0\\ a_{21} & a_{22} & ... & 0 \\ ... & ... & ... & ...\\a_{n1} & a_{n2} & ... & a_{nn}\end{pmatrix} \begin{pmatrix}x_1\\x_2\\...\\x_n\end{pmatrix}=\begin{pmatrix}b_1\\b_2\\...\\b_n\end{pmatrix}\]

To wtedy rozwiązaniem jest wektor, który uzyskuje się przez zastosowanie algorytmu zastąpienia w przód (ang. forward substitution):

\[\begin{pmatrix}x_1\\x_2\\...\\x_n\end{pmatrix}=\begin{pmatrix}\frac{b_1}{a_{11}}\\\frac{b_2-a_{21}x_1}{a_{22}}\\...\\\frac{b_n-a_{n1}x1-a_{n2}x2-...}{a_{nn}}\end{pmatrix}\]

Macierz współczynników ma postać górno-trójkątną \[\begin{pmatrix}a_{11} & a_{12} & ... & a_{1n}\\ 0 & a_{22} & ... & a_{2n} \\ ... & ... & ... & ...\\0 & 0 & ... & a_{nn}\end{pmatrix} \begin{pmatrix}x_1\\x_2\\...\\x_n\end{pmatrix}=\begin{pmatrix}b_1\\b_2\\...\\b_n\end{pmatrix}\]

To wtedy rozwiązaniem jest wektor, który uzyskuje się przez zastosowanie algorytmu zastąpienia w tył (ang. back substitution): \[\begin{pmatrix}x_n\\...\\x_2\\x_1\end{pmatrix}=\begin{pmatrix}\frac{b_n}{a_{nn}}\\...\\\frac{b_2-...-a_{2n}x_n}{a_{22}}\\\frac{b_1-a_{12}x_2-...-a_{1n}x_n}{a_{11}}\end{pmatrix}\]

Macierz współczynników ma postać permutacyjnej macierzy dolno-trójkątnej, bądź górno-trójkątnej

\[\begin{pmatrix}a_{11} & 0 & ... & 0\\ ... & ... & ... & ...\\ a_{n1} & a_{n2} & ... & a_{nn}\\a_{21} & a_{22} & ... & 0 \end{pmatrix} \begin{pmatrix}x_1\\...\\x_n\\x_2\end{pmatrix}=\begin{pmatrix}b_1\\b_n\\...\\b_2\end{pmatrix}\]

Znając wektor permutacji macierzy, możemy znaleźć jej postać dolno-trójkątną, bądź górno-trójkątną. Następnie stosujemy algorytm zastąpienia w przód, bądź w tył.

Uwaga: W algorytmach zastąpienia w przód i w tył mamy dzielenie przez współczynniki znajdujące się na przekątnej macierzy. One muszą być różne od zera, aby rozwiązanie istniało. No i oczywiście macierz odwrotna musi też istnieć.

Metoda polega na rozłożeniu na czynniki macierzy tj. znalezieniu dolno-trójkątnej (L) i górno-trójkątnej (U) macierzy, których iloczyn daje macierz współczynników głównego równania: \[A=LU\]

Nasz układ równań możemy przedstawić w postaci: \[Lz=b,\text{szukamy z}\] \[Ux=z,\text{szukamy x}\] \[z=L^{-1}b\] \[Ux=L^{-1}b\] \[LUx=b\]

\[L=\begin{pmatrix}l_{11} & 0 & ... & 0 \\ l_{21} & l_{22} & ... & 0 \\ ... & ... & ... & ...\\l_{n1} & l_{n2} & ... & l_{nn}\end{pmatrix} \]

\[U=\begin{pmatrix}u_{11} & u_{12} & ... & u_{1n}\\ 0 & u_{22} & ... & u_{2n} \\ ... & ... & ... & ...\\0 & 0 & ... & u_{nn}\end{pmatrix} \]

Twierdzenie: jeśli każdy wiodący minor główny macierzy A stopnia n jest różny od zera (czyli taki minor jest nieosobliwy), to macierz A jest rozkładalna na czynniki L i U.

\[A_k=\begin{pmatrix}a_{11} & a_{12} & ... & a_{1k} \\ a_{21} & a_{22} & ... & a_{2k} \\ ... & ... & ... & ...\\a_{k1} & a_{k2} & ... & a_{kk}\end{pmatrix},1\le k \le n \]

Jeśli macierz posiada taki rozkład, to nie znaczy jeszcze, że ten rozkład jest jednoznaczny. Znajdowanie rozkładu zaczynamy od przypisania wartości niezerowych współczynnikom na przekątnej w tylko jednej macierzy np. dla \(l_{ii}=1, \text{ dla } 1\leq i\leq n\) - wtedy taką macierz nazywamy jednostkową macierzą dolno-trójkątną, dla \(u_{ii}=1, \text{ dla } 1\leq i\leq n\) - jednostkowa macierz górno-trójkątna.

Podstawowe wzory do zbudowania algorytmu:
\[a_{kk}=\sum_{s=1}^{k-1} l_{ks}u_{sk}+l_{kk}u_{kk}\] \[a_{kj}=\sum_{s=1}^{k-1} l_{ks}u_{sj}+l_{kk}u_{kj}\] \[a_{ik}=\sum_{s=1}^{k-1} l_{is}u_{sk}+l_{ik}u_{kk}\]

• jeśli \(l_{kk}\neq0\) to obliczamy u_kj z drugiego równania a l_ik z trzeciego: \[u_{kj} = a_{kj}-\sum_{s=1}^{k-1} l_{ks}u_{sj},\text{ dla } k\leq j\leq n\] \[l_{ik} = \frac{a_{ik}-\sum_{s=1}^{k-1} l_{is}u_{sk}}{u_{kk}},\text{ dla } (k+1)\leq i \leq n\]
• jeśli \(u_{kk}\neq0\) to obliczamy l_ik z trzeciego równania z u_kj z drugiego: \[l_{ik}=a_{ik}-\sum_{s=1}^{k-1} l_{is}u_{sk},\text{ dla }k\leq i \leq n\] \[u_{kj}=\frac{a_{kj}-\sum_{s=1}^{k-1} l_{ks}u_{sj}}{l_{kk}},\text{ dla } (k+1)\leq j\leq n\]

• jeśli macierz jest symetryczna i dodatnio określona to l_kk=u_kk i l_ks=u_sk: \[l_{kk}=\sqrt{a_{kk}-\sum_{s=1}^{k-1} l_{ks}^2}\]

  \[l_{ik} = \frac{a_{ik}-\sum_{s=1}^{k-1} l_{is}l_{ks}}{l_{kk}},\text{ dla } (k+1)\leq i \leq n\]

Algorytm oparty na jednostkowej macierzy dolno-trójkątnej L nazywamy metodą Doolittle‘a, a oparty na jednostkowej macierzy górno-trójkątnej U metodą Crouta. Gdy L=U^T to macierz jest symetryczna, a więc zachodzi równość l_ks=u_sk, dla dowolnych k i s - algorytm oparty na takiej własności nazywamy metodą Choleskiego. Aby zastosować metodę Choleskiego macierz oprócz tego, że musi być symetryczna, to musi być także dodatnio określona tj. \(x^TAx\geq0,\text{ dla wszystkich } x\in\mathbb R^n\smallsetminus \{0\}\)

Przykład macierzy symetrycznej, która nie jest dodatnio określona:
\[\begin{pmatrix}1 & 2 \\ 2 & 1\end{pmatrix}\] \[\begin{pmatrix}-1 & 1\end{pmatrix}\begin{pmatrix}1 & 2 \\ 2 & 1\end{pmatrix}\begin{pmatrix}-1 \\ 1\end{pmatrix}=-2<0\]

Twierdzenie: Jeśli macierz A jest rzeczywista, symetryczna i dodatnio określona, to posiada jednoznaczny rozkład na czynniki L i U \(A=LU=LL^T\), gdzie L jest dolno-trójkątna z dodatnimi współczynnikami na przekątnej

Kryterium Sylvestera: macierz symetryczna A jest dodatnio określona wtedy i tylko wtedy, gdy wszystkie jej wiodące minory główne są dodatnie

using LinearAlgebra
isposdef(A) # funkcja do sprawdzania, czy macierz A jest dodatnio określona
issymmetric(A) # funkcja do sprawdzania, czy macierz A jest symetryczna

Informacje o metodach do obliczeń rozkładu na czynniki macierzy współczynników w dokumentacji Julki: https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#man-linalg-factorizations

Układ równań liniowych 4x4 ma postać:

\[\begin{pmatrix}a_{11} & a_{12} & a_{13} & a_{14}\\ a_{21} & a_{22} & a_{23} & a_{24} \\ a_{31} & a_{32} & a_{33} & a_{34}\\a_{41} & a_{42} & a_{43} & a_{44}\end{pmatrix} \begin{pmatrix}x_1\\x_2\\x_3\\x_4\end{pmatrix}=\begin{pmatrix}b_{1}\\b_{2}\\b_{3}\\b_{4}\end{pmatrix}\]

Skrócona postać macierzowa: \(Ax=b\)

Przykład:

\[\begin{pmatrix}3 & 2 & 1 & 2\\ 6 & 8 & 1 & -1 \\ -3 & 2 & 2 & -2\\9 & 2 & 6 & 1\end{pmatrix} \begin{pmatrix}x_1\\x_2\\x_3\\x_4\end{pmatrix}=\begin{pmatrix}1\\3\\2\\-1\end{pmatrix}\]

W metodzie eliminacji Gaussa naszym celem jest, poprzez zastosowanie działań elementarnych na wierszach macierzy współczynników A, znaleźć rozkład na czynniki \(A=LU\), by rozwiązać układ równań liniowych \(LUx=b\) w takiej postaci: \(Ux=L^{-1}b\)

Działania elementarne na wierszach macierzy to (A_i i A_j oznaczają wiersze macierzy A):

1. Zamienianie miejscami wierszy \(A_i \leftrightarrow A_j\)
2. Mnożenie współczynników wiersza macierzy przez niezerową liczbę \(\lambda A_i \rightarrow A_i\)
3. Dodawanie do danego wiersza drugi przemnożony przez niezerową liczbę \(A_j + \lambda A_i \rightarrow A_j\)

W tej najprostszej postaci metody wykorzystamy tylko działanie elementarne nr 3. Działamy na wierszach. Pierwszy wiersz pozostaje niezmieniony. Pozostałe wiersze tak modyfikujemy, aby wyzerować kolumnę nr 1. Krok nr 1: \(A_2-2A_1\); \(A_3-(-1A_1)\);\(A_4-3A_1\); Te działania dotyczą także wektora współczynników b. \[\begin{pmatrix}3 & 2 & 1 & 2\\ 0 & 4 & -1 & -5 \\ 0 & 4 & 3 & 0\\0 & -4 & 3 &-5\end{pmatrix}\begin{pmatrix}x_1\\x_2\\x_3\\x_4\end{pmatrix}=\begin{pmatrix}1\\1\\3\\-4\end{pmatrix}\]

Element a₁₁=3 nazywamy współczynnikiem prowadzącym (ang. pivot element). Wiersz A₁ nazywamy wierszem prowadzących (ang. pivot row). Współczynniki, przez które przemnażane są wiersze, nazywamy mnożnikami (ang. multipliers)

Krok nr 2. \(A_3-1A_2\); \(A_4-(-1A_2)\); \[\begin{pmatrix}3 & 2 & 1 & 2\\ 0 & 4 & -1 & -5 \\ 0 & 0 & 4 & 5\\0 & 0 & 2 &-10\end{pmatrix}\begin{pmatrix}x_1\\x_2\\x_3\\x_4\end{pmatrix}=\begin{pmatrix}1\\1\\2\\-3\end{pmatrix}\]

Element a₂₂=4 jest współczynnikiem prowadzącym.

Krok nr 3. \(A_4-0.5A_3\)

\[\begin{pmatrix}3 & 2 & 1 & 2\\ 0 & 4 & -1 & -5 \\ 0 & 0 & 4 & 5\\0 & 0 & 0 &-12.5\end{pmatrix}\begin{pmatrix}x_1\\x_2\\x_3\\x_4\end{pmatrix}=\begin{pmatrix}1\\1\\2\\-4\end{pmatrix}\]

Element a₃₃=4 jest współczynnikiem prowadzącym. W kroku nr 3 uzyskaliśmy macierz górno-trójkątną.

\[U=\begin{pmatrix}3 & 2 & 1 & 2\\ 0 & 4 & -1 & -5 \\ 0 & 0 & 4 & 5\\0 & 0 & 0 &-12.5\end{pmatrix}\]

Macierz dolno-trójkąta jest złożona ze współczynników, przez które przemnażane były wiersze: \[L=\begin{pmatrix}1 & 0 & 0 & 0\\ 2 & 1 & 0 & 0 \\ -1 & 1 & 1 & 0\\3 & -1 & 0.5 &1\end{pmatrix}\]

\(A=LU\):

\[\begin{pmatrix}3 & 2 & 1 & 2\\ 6 & 8 & 1 & -1 \\ -3 & 2 & 2 & -2\\9 & 2 & 6 & 1\end{pmatrix}=\begin{pmatrix}1 & 0 & 0 & 0\\ 2 & 1 & 0 & 0 \\ -1 & 1 & 1 & 0\\3 & -1 & 0.5 &1\end{pmatrix}\begin{pmatrix}3 & 2 & 1 & 2\\ 0 & 4 & -1 & -5 \\ 0 & 0 & 4 & 5\\0 & 0 & 0 &-12.5\end{pmatrix}\]

A=[3 2 1 2;6 8 1 -1;-3 2 2 -2;9 2 6 1]
U=[3 2 1 2;0 4 -1 -5;0 0 4 5;0 0 0 -12.5]
L=[1 0 0 0;2 1 0 0;-1 1 1 0; 3 -1 0.5 1]
L*U == A

Nasz algorytm składa się z następujących kroków: \[\begin{pmatrix}a_{11} & a_{12} & a_{13} & a_{14}\\ a_{21} & a_{22} & a_{23} & a_{24} \\ a_{31} & a_{32} & a_{33} & a_{34}\\a_{41} & a_{42} & a_{43} & a_{44}\end{pmatrix}\]

function gauss(n,A)
    # dla k od 1 do n-1:
    #     dla i od k+1 do n
    #         z = u_{ik}/u_{kk} 
    #         u_{ik} = 0
    #         dla j od k+1 do n
    #               u_{ij}=u_{ij}-z*u_{kj}
    return U
end

Niezerowe współczynniki prowadzące stanowią niejedyne ograniczenie w stosowaniu wyżej opisanej metody.

Współczynnik prowadzący jest 0:
\[\begin{pmatrix}0 & 1 \\ 1 & 1 \end{pmatrix}\begin{pmatrix}x_1\\x_2\end{pmatrix}=\begin{pmatrix}1\\2\end{pmatrix}\] Współczynnik prowadzący jest bliski zeru:

\[\begin{pmatrix}\epsilon & 1 \\ 1 & 1 \end{pmatrix}\begin{pmatrix}x_1\\x_2\end{pmatrix}=\begin{pmatrix}1\\2\end{pmatrix}\]

Po zastosowaniu metody eliminacji Gaussa otrzymamy: \[\begin{pmatrix}\epsilon & 1 \\ 0 & 1-\epsilon^{-1} \end{pmatrix}\begin{pmatrix}x_1\\x_2\end{pmatrix}=\begin{pmatrix}1\\2-\epsilon^{-1}\end{pmatrix}\]

Rozwiązaniem analitycznym jest: \(x_2=\frac{2-\epsilon^{-1}}{1-\epsilon^{-1}}\approx 1\) \(x_1=\frac{1-x_2}{\epsilon}\approx 1\)

Ale rozwiązaniem komputerowym dla \(\epsilon < eps(Float64)\) będzie: \(x_2=\frac{2-\epsilon^{-1}}{1-\epsilon^{-1}}= 1\); \(x_1=\frac{1-x_2}{\epsilon}= 0\)

Problem generuje nam \(\epsilon^{-1}\), które to wyrażenie jest bardzo dużą liczbą. Problem jest związany z dokładnością reprezentacji liczb podczas ich odejmowania, tj. \(2-\epsilon^{-1}\) i \(1-\epsilon^{-1}\)

 Int64(0b10000110110)-1023 # –> 55
 bitstring(-1/0.00000000000000002) == bitstring(2-1/0.00000000000000002) # –> true
 Int64(0b10000110101)-1023 # –> 54
 bitstring(-1/0.00000000000000005) == bitstring(2-1/0.00000000000000005) # –> true

 Int64(0b10000110100)-1023 # –> 53
 bitstring(-1/0.00000000000000009) == bitstring(2-1/0.00000000000000009) # –> false

 Int64(0b10000110011)-1023 # —> 52
 bitstring(-1/0.0000000000000002) == bitstring(2-1/0.0000000000000002) # –> false


bitstring(-0.5e17) == bitstring(2-0.5e17) # –> true
bitstring(-0.5e16) == bitstring(2-0.5e16) # –> false

bitstring(-1/0.00000000000000002)
bitstring(-0.5e17)
# "1100001101100110001101000101011110000101110110001010000000000000"
bitstring(2-1/0.00000000000000002)
bitstring(2-0.5e17)
# "1100001101100110001101000101011110000101110110001010000000000000"

# Komputer odejmując dwie liczba sprowadza je do tej samej unormowanej
# reprezentacji naukowej, czyli mamy takie działanie: 
0.00000000000000002*1e17 - 0.5*1e17 # –> -0.5e17 
# Dostajemy taki wynik, gdyż liczba 2 w tej reprezentacji jest poza zakresem
# dokładności reprezentowania liczby zminnoprzecinkowej i dlatego nam się zeruje
# podczas próby odejmownia.

# W tym przypadku liczba 2 nie będzie zerowana, gdyż jest w zakresie dokładności:
# eps(Float64) –> 2.220446049250313e-16
0.0000000000000002*1e16 - 0.5*1e16 # –> -0.4999999999999998e16

bitstring(-1/0.0000000000000002)
bitstring(-0.5e16)
# "1100001100110001110000110111100100110111111000001000000000000000" 
bitstring(2-1/0.0000000000000002)
bitstring(2-0.5e16)
# "1100001100110001110000110111100100110111111000000111111111111110"

W celu rozwiązania tego problemu do naszego wcześniejszego algorytmu dodamy wyszukiwanie względnie największego współczynnika prowadzącego w danym wierszu. Na podstawie określenia dla każdego wiersza takiego współczynnika stworzymy wektor permutacji P, który będzie określał kolejność działań na wierszach. \(PAx=Pb\)

Nasz algorytm podzielmy na dwie części: znajdowania rozkładu na czynniki (ang. factorization) oraz znajdowania rozwiązania równania (ang. solution). Pierwsza część będzie oparta o algorytm eliminacji w przód (ang. forward elimination), a druga część o algorytm zastąpienia w tył (ang. back substitution)

W celu znalezienia wektora permutacji szukamy maksymalnego elementu w każdym wierszu: \(s_i = max{(|a_{i1}|,|a_{i2}|,...,|a_{in}|)},1 \leq i \leq n\)

a następnie szukamy względnie największych współczynników prowadzących: \(P_i=|a_{i1}|/s_i\), na których podstawie tworzymy wektor permutacji P tj. każdy kolejny element wektora P jest mniejszy od poprzedniego: \(P_1 \geq P_2 \geq ... \geq P_n\)

S=zeros(n);
for i in 1:n
    S[i]=findmax(abs.(A[i,:]))[1];
end
P=sortperm(abs.(A[:,1])./S,rev=true)

\[\begin{pmatrix}3 & 2 & 4\\ 6 & 8 & 1 \\ -6 & 2 & 2\end{pmatrix} \begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}=\begin{pmatrix}1\\2\\3\end{pmatrix}\]

A=[3 2 4;6 8 1;-6 2 2]
P=[1;2;3]
S=zeros(n);
  for i in 1:n
        S[i]=findmax(abs.(A[i,:]))[1]; # –> S={4,8,6}
  end
abs.(A[:,1])./S # –> 0.75,0.75,1.0
pivot_list=sortperm(abs.(A[:,1])./S,rev=true) # –> P={3,1,2}

# pierwszy wiersz prowadzący to 3; podmeniamy odpowiednio w P wiersz 1 z wierszem 3
P=[3;2;1]

Zgodnie w wektorem permutacji naszym pierwszym wierszem prowadzącym jest wiersz nr 3.

1. \(A_2-(-1A_3),A_1-(-0.5A3)\)

\[\begin{pmatrix}0 & 3 & 5\\ 0 & 10 & 3 \\ -6 & 2 & 2\end{pmatrix} \begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}=\begin{pmatrix}2.5\\5\\3\end{pmatrix}\]

Kolejnym wierszem jest wiersz nr 1

S = [4;8;6]
P=[3;2;1]
A=[0 3 5;0 10 3;-6 2 2]
abs.(A[P[2:3],2])./S[P[2:3]] # –> 1.25,0.75
pivot_list=sortperm(abs.(A[P[2:3],1])./S[P[2:3]],rev=true) # –> P={1,2}
# drugi wiersz prowadzący to 1; podmeniamy odpowiednio w P wiersz 2 z wierszem 1
P=[3;1;2]

1. \(A_2-10/3A_1)\)

\[\begin{pmatrix}0 & 3 & 5\\ 0 & 0 & -41/3 \\ -6 & 2 & 2\end{pmatrix} \begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}=\begin{pmatrix}2.5\\-10/3\\3\end{pmatrix}\]

Znaleźliśmy macierz permutacyjną górno-trójkątną U: \[U=\begin{pmatrix}0 & 3 & 5\\ 0 & 0 & -41/3 \\ -6 & 2 & 2\end{pmatrix}\]

Macierz permutacyjną dolno-trójkątną stanowią mnożniki: \[L=\begin{pmatrix}-0.5 & 1 & 0\\ -1 & 10/3 & 1 \\ 1 & 0 & 0\end{pmatrix}\]

czyli mamy: \[PA= \begin{pmatrix}-6 & 2 & 2\\ 3 & 2 & 4 \\ 6 & 8 & 1\end{pmatrix} = \begin{pmatrix}1 & 0 & 0\\ -0.5 & 1 & 0 \\ -1 & 10/3 & 1\end{pmatrix}\begin{pmatrix}-6 & 2 & 2\\ 0 & 3 & 5 \\ 0 & 0 & -41/3\end{pmatrix}\] \[Ux=L^{-1}b = \begin{pmatrix}-6 & 2 & 2\\ 0 & 3 & 5 \\ 0 & 0 & -41/3\end{pmatrix}\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}=\begin{pmatrix}1 & 0 & 0\\ -0.5 & 1 & 0 \\ -1 & 10/3 & 1\end{pmatrix}^{-1}\begin{pmatrix}3\\1\\2\end{pmatrix}\]

Sprawdzenie zgodności obliczenia rozkładu na czynniki L i U:

A=[-6.0 2.0 2.0;3.0 2.0 4.0;6.0 8.0 1.0];
B=[3,1,2];
L=[1 0 0;-0.5 1 0;-1 10/3 1];
U=[-6 2 2;0 3 5;0 0 -41/3];
L*U

3×3 Matrix{Float64}:
 -6.0  2.0  2.0
  3.0  2.0  4.0
  6.0  8.0  1.0

Sprawdzenie zgodności obliczeń współczynników b

B=[3,1,2] ;
L=[1 0 0;-0.5 1 0;-1 10/3 1];
inv(L)*B

3-element Vector{Float64}:
  3.0
  2.5
 -3.333333333333334

Pseudokod:

  function factorization(n,A)
        for i in 1:n
            S[i]=findmax(abs.(U[i,:]))[1];
        end
        P=sortperm(abs.(A[:,1])./S,rev=true)
          # dla k od 1 do n-1:
          #     dla i od k+1 do n
          #         z = u_{pi,k}/u_{pk,k} 
          #         zapisujemy do tablicy mnożniki zamiast 0
          #         u_{pi,k} = z
          #         dla j od k+1 do n
          #                 u_{pi,j}=u_{pi,j}-z*u_{pk,j}

        return U,P
  end

  function solution(n,U,P,B)
        # dla k od 1 do (n-1)
        #     dla i od (k+1) do n
        #         b_{pi}=b_{pi}-u{pi,k}*b_{pk}
      # stosujemy algorymt zliczania w tył:       
        # x_{n}=b_{pn}/u_{pn,n}
        # dla i od (n-1), o krok -1, do 1
        #     x_{i}=(b_{pi}-sum(u_{pi,j}*x_{j} dla j od (i+1) do n)/u_{pi,i}

        return X
  end

function gausspivot(n,A,B)
    U,P=factorization(n,A);
    X=solution(n,U,P,B);
    return X
end

Problem znajdowania rozwiązania równań liniowych definiowaliśmy jako: \(Ax=b\). W metodach iteracyjnych modyfikujemy naszą definicje, dodając do niej macierz dzieląca (ang. splitting matrix): \(Qx = (Q-A)x + b\), które jest równoznaczne z podstawą definicją, gdyż \(Qx = Qx-Ax + b \rightarrow Ax=b\)

W postaci iteracyjnej nasza nowa definicja problemu wygląda następująco: \[Qx^{(k)} = (Q-A)x^{(k-1)} + b\,\text{ dla }k\geq 1 \]

Po dwóch stronach równania otrzymaliśmy zmienną x. W przeciwieństwie do metody eliminacji Gaussa, która umożliwia otrzymanie rozwiązania ścisłego, to w metodach iteracyjnych uzyskujemy rozwiązanie przybliżone.

Wektor początkowy x⁰ może być dowolny. Pytanie więc, jakie warunki musi spełniać A i Q, aby dla dowolnego wektora początkowego rozwiązanie x^(k) istniało, a tym samym algorytm był zbieżny: \(\lim\limits_{k\rightarrow \infty}||x^{(k)}-x||=0\)

Po pierwsze macierze Q i A muszą być nieosobliwe, tj. istnieją dla nich macierze odwrotne. Kolejny warunek jest opisany w twierdzeniu poniżej.

Gdy przemnożone zostanie równanie przez macierz odwrotną do Q, to otrzymamy: \[x^{(k)} = (I-Q^{-1}A)x^{(k-1)} + Q^{-1}b\]

Twierdzenie o zbieżności metod iteracyjnych: jeśli \(||I-Q^{-1}A||<1\) dla
pewnej podporządkowanej normy macierzowej, to dla dowolnego wektora początkowego
rezultat x^k działania algorytmu zbiega do rozwiązania równania Ax=b.

Przykładowe podporządkowane normy macierzowe:
Norma L¹: \(||A||_{1}=\max\limits_{1\leq j\leq n}\sum_{i=1}^{n}|a_{ij}|\)

Norma L² - norma Euklidesowa: \(||A||_{2}=\max\limits_{1\leq i\leq n} |\sigma_{i}|=\sqrt{\rho(A^{T}A)}\), \(\rho(A^{T}A)\) - promień spektralny, czyli największa wartość własna \(\lambda\) macierzy \(A^{T}A\)

Norma L^∞: \(||A||_{\infty}=\max\limits_{1\leq i\leq n}\sum_{j=1}^{n}|a_{ij}|\)

# l_1 
function l1(n,A)
        findmax(sum(abs.(A[i,:]) for i in 1:n))
end
# l_2
function l2(A)
        sqrt(findmax(eigvals(A'*A))[1])
end
function linf(n,A)
        # l_inf
        findmax(sum(abs.(A[:,i]) for i in 1:n))
end

Przykład dla normy L¹

\[\begin{pmatrix}5 & 2 &1\\ 5 & 7 &1\\ 3 &2 &6 \end{pmatrix}\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}=\begin{pmatrix}1\\2\\3\end{pmatrix}\]

Wiersz 1: |5|+|2|+|1|=8
Wiersz 2: |5|+|7|+|1|=13
Wiersz 3: |3|+|2|+|6|=8

\(||A||_{1} = 13\)

Norma L¹ mówi, jak bardzo macierz A może rozciągnąć lub skurczyć wektory. W naszym przypadku macierz A maksymalnie może rozciągnąć wektor jednostkowy aż do długości 13.

Zbadajmy przypadek rozciągania wektorów \(e_{1}=[\frac{1}{3},\frac{1}{3},\frac{1}{3}]\) oraz \(e_{2}=[1,0,0]\)

\(y_{1}= A \cdot e_{1} = \begin{pmatrix}2.67\\4.33\\3.67\end{pmatrix}\)
\(y_{2}= A \cdot e_{2} = \begin{pmatrix}5\\5\\3\end{pmatrix}\)

\(||y_{1}||_{1} = 10.67\)
\(||y_{2}||_{1} = 13\)

W Julce normy macierzowe możemy obliczyć z wykorzystaniem opnorm: https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.opnorm

Przykład: \[\begin{pmatrix}5 & 2 &1\\ 5 & 7 &1\\ 3 &2 &6 \end{pmatrix}\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}=\begin{pmatrix}1\\2\\3\end{pmatrix}\]

Pytanie: czy w tym przykładzie spełnione jest twierdzenie o zbieżności metod iteracyjnych, tj. czy $||I-Q^-1A||<1 ?

W metodzie Jakobiego Q jest macierz diagonalną macierzy A tzn. Q=diag(A), więc nasze równanie możemy zapisać w postaci \(Qx = (Q-A)x + b\):

\[\begin{pmatrix}5 & 0 &0\\ 0 & 7 &0\\ 0 &0 &6 \end{pmatrix}\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}=\Biggl[\begin{pmatrix}5 & 0 &0\\ 0 & 7 &0\\ 0 &0 &6 \end{pmatrix}-\begin{pmatrix}5 & 2 &1\\ 5 & 7 &1\\ 3 &2 &6 \end{pmatrix}\Biggr]\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}+\begin{pmatrix}1\\2\\3\end{pmatrix}\]

\[\begin{pmatrix}5 & 0 &0\\ 0 & 7 &0\\ 0 &0 &6 \end{pmatrix}\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}=\Biggl[\begin{pmatrix}0 & -2 &-1\\ -5 & 0 &-1\\ -3 &-2 &0 \end{pmatrix}\Biggr]\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}+\begin{pmatrix}1\\2\\3\end{pmatrix}\]

\[x_{1}^{(k)}=\frac{1}{5}-\frac{2}{5}x_{2}^{(k-1)}-\frac{1}{5}x_{3}^{(k-1)}\] \[x_{2}^{(k)}=\frac{2}{7}-\frac{5}{7}x_{1}^{(k-1)}-\frac{1}{7}x_{3}^{(k-1)}\] \[x_{3}^{(k)}=\frac{3}{6}-\frac{3}{6}x_{1}^{(k-1)}-\frac{2}{6}x_{2}^{(k-1)}\]

A=[5 2 1;5 7 1;3 2 6];
Q=one(A) .* diag(A);

Po przekształceniach do postaci \(x = (I-Q^{-1}A)x + Q^{-1}b\) mamy: \[\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}=\Biggl[\begin{pmatrix}1 & 0 &0\\ 0 & 1 &0\\ 0 &0 &1 \end{pmatrix}-\begin{pmatrix}1/5 & 0 &0\\ 0 & 1/7 &0\\ 0 &0 &1/6 \end{pmatrix}\begin{pmatrix}5 & 2 &1\\ 5 & 7 &1\\ 3 &2 &6 \end{pmatrix}\Biggr]\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}+\begin{pmatrix}1/5 & 0 &0\\ 0 & 1/7 &0\\ 0 &0 &1/6 \end{pmatrix}\begin{pmatrix}1\\2\\3\end{pmatrix}\]

Sprawdzamy, czy spełnione jest twierdzenie o zbieżności metod iteracyjnych:

Ione = one(A);
Qinv=inv(Q);
linf(3,Ione-Qinv*A)[1]<1 # true

\[I-Q^{-1}A = \begin{pmatrix}0 & -2/5 &-1/5\\ -5/7 & 0 & -1/7\\ -3/6 &-2/6 &0\end{pmatrix}\]

Zgodnie z twierdzeniem o zbieżności metod iteracyjnych stwierdzanym, że rezultat działania algorytmu będzie zbiegał do rozwiązania równania Ax=b dla dowolnych warunków początkowych.

Działanie algorytmu możemy rozpocząć od wybrania losowych wartości początkowych \(x_{1}^{(0)}\) i \(x_{2}^{(0)}\), \(x_{3}^{(0)}\) . Rezultatem działania algorytmu jest wektor \((x_{1}^{(k)},x_{2}^{(k)},x_{3}^{(k)})^{T}\), który reprezentuje przybliżone rozwiązanie układu tj. takie, którego odległość od rozwiązania dla k-1 jest mniejsza od oczekiwanej dokładności \(\epsilon\): \(||x^{(k)}-x^{(k-1)}||<\epsilon\)

function jakobi(n,A,B,X,M,Eps)
    # n - wymiar macierzy
    # A - mecierz współczynników równania liniowego
    # B - wektor wyrazów wolnych 
    # X - wektor niewiadomych
    # M - maksymalna ilość iteracji
    # Eps - dopuszczalny błąd
    # dla k od 1 do M
    #     dla i od 1 do n
    #         U_i = (B_i-sum(A_{ij}*X_j dla j od 1 do n jeśli j różne od i))/A_{i,i}
    #     jeśli sqrt(sum((X-U).^2))<Eps) zwróc X
    #     dla i od 1 do n
    #         X_i = U_i
    # zwróć "Algorithm failed to converge"
end

Poniższe twierdzenie jest wersją twierdzenia o zbieżności metod iteracyjnych dla metody Jakobiego:

Twierdzenie o zbieżności metody Jakobiego: Jeśli macierz A jest diagonalnie
dominująca, to rezultat działania metody Jakobiego zbiega do rozwiązania Ax=b dla
dowolnej wartości początkowej.

Macierz diagonalnie dominująca to taka, dla której spełniona jest nierówność: \[|a_{ii}|>\sum_{j=1,j\neq i}^{n} |a_{ij}|,\, \text{dla } 1\leq i\leq n\]

Aby uniknąć dzielenia w algorytmie opisanym przez wzory na początku tego podrozdziału, możemy przekształcić nasze równanie, przemnażając obie strony przez macierz odwrotną D=diag(A) tj: \(D^{-1}Ax=D^{-1}b\)

A=[5 2 1;5 7 1;3 2 6];
D=one(A) .* diag(A);
Dinv=inv(D);
Dinv*A

\(D^{-1}Ax=D^{-1}b\rightarrow\begin{pmatrix}1 & 2/5 &1/5\\ 5/7 & 1 & 1/7\\ 3/6 &2/6 &1\end{pmatrix}\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}=\begin{pmatrix}1/5\\2/7\\3/6\end{pmatrix}\)

Dzięki temu, że na przekątnej są wartości 1, to do algorytmu nie musimy dodawać operacji dzielenia podczas obliczania x. W celu zrealizowania tego przemnożenia przez macierz odwrotną do Q=diag(A), to do naszego algorytmu dodajemy poniższy kod i odpowiednio modyfikujemy pozostałą część. W algorytmie nadal jest operacja dzielenia, ale w wyrażeniu o wiele mniej złożonym.

# dla i od 1 do n
#     d=1/A_{i,i}
#     B_i=d*B_i;
#     dla j od 1 do n
#         A_{i,j}=d*A_{i,j};

Ta metoda działa bardzo podobnie do metody Jakobiego. Różnica jest taka, że w metodzie Gaussa-Seidla do obliczania kolejnej niewiadomej np x₂ wykorzystujemy obliczenia poprzedniej x₁ oraz macierz Q jest macierzą dolno-trójkątną Q=L(A).

\[\begin{pmatrix}5 & 2 &1\\ 5 & 7 &1\\ 3 &2 &6 \end{pmatrix}\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}=\begin{pmatrix}1\\2\\3\end{pmatrix}\]

Twierdzenie o zbieżności metody Gaussa-Seidla: Jeśli macierz A jest diagonalnie
dominująca, to rezultat działania metody Gauss-Siedla zbiega do rozwiązania Ax=b dla
dowolnej wartości początkowej.

\[\begin{pmatrix}5 & 0 &0\\ 5 & 7 &0\\ 3 &2 &6 \end{pmatrix}\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}=\Biggl[\begin{pmatrix}5 & 0 &0\\ 5 & 7 &0\\ 3 &2 &6 \end{pmatrix}-\begin{pmatrix}5 & 2 &1\\ 5 & 7 &1\\ 3 &2 &6 \end{pmatrix}\Biggr]\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}+\begin{pmatrix}1\\2\\3\end{pmatrix}\]

\[\begin{pmatrix}5 & 0 &0\\ 5 & 7 &0\\ 3 &2 &6 \end{pmatrix}\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}=\begin{pmatrix}0 & -2 &-1\\ 0 & 0 &-1\\ 0 &0 &0 \end{pmatrix}\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}+\begin{pmatrix}1\\2\\3\end{pmatrix}\]

\[x_{1}^{(k)}=\frac{1}{5}-\frac{2}{5}x_{2}^{(k-1)}-\frac{1}{5}x_{3}^{(k-1)}\] \[x_{2}^{(k)}=\frac{2}{7}-\frac{5}{7}x_{1}^{(k)}-\frac{1}{7}x_{3}^{(k-1)}\] \[x_{3}^{(k)}=\frac{3}{6}-\frac{3}{6}x_{1}^{(k)}-\frac{2}{6}x_{2}^{(k)}\]

A=[5 2 1;5 7 1;3 2 6];
Q=[5 0 0;5 7 0; 3 2 6];

Po przekształceniach do postaci \(x = (I-Q^{-1}A)x + Q^{-1}b\) mamy: \[\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}=\Biggl[\begin{pmatrix}1 & 0 &0\\ 0 & 1 &0\\ 0 &0 &1 \end{pmatrix}-\begin{pmatrix}1/5 & 0 &0\\ -1/7 & 1/7 &0\\ -0.052381 &-0.047619 &1/6 \end{pmatrix}\begin{pmatrix}5 & 2 &1\\ 5 & 7 &1\\ 3 &2 &6 \end{pmatrix}\Biggr]\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}+\begin{pmatrix}1/5 & 0 &0\\ -1/7 & 1/7 &0\\ -0.052381 &-0.047619 &1/6 \end{pmatrix}\begin{pmatrix}1\\2\\3\end{pmatrix}\]

Ione = one(A);
Qinv=inv(Q);
linf(3,Ione-Qinv*A)[1]<1 # true
Ione-Qinv*A

\[I-Q^{-1}A = \begin{pmatrix}0 & -2/5 &-1/5\\ 0 & -2/7 & 0\\ 0 & 0.104762 &1/10\end{pmatrix}\]

Zgodnie z twierdzeniem o zbieżności metod iteracyjnych stwierdzanym, że rezultat działania algorytmu będzie zbiegał do rozwiązania równania Ax=b dla dowolnych warników początkowych.

Podobnie jak w metodzie Jakobiego, w celu uniknięcia dodawania operacji dzielenia podczas obliczania x przemnażamy przez macierz odwrotną do D=diag(A).

A=[5 2 1;5 7 1;3 2 6];
D=one(A) .* diag(A);
Dinv=inv(D);
Dinv*A

\(D^{-1}Ax=D^{-1}b\rightarrow\begin{pmatrix}1 & 2/5 &1/5\\ 5/7 & 1 & 1/7\\ 3/6 &2/6 &1\end{pmatrix}\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}=\begin{pmatrix}1/5\\2/7\\3/6\end{pmatrix}\)

Q=[1 0 0;5/7 1 0; 3/6 2/6 1];
Qinv=inv(Q)

\[\begin{pmatrix}1 & 0 &0\\ 5/7 & 1 & 0\\ 3/6 &2/6 &1\end{pmatrix}\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}=\Biggl[\begin{pmatrix}1 & 0 &0\\ 5/7 & 1 & 0\\ 3/6 &2/6 &1\end{pmatrix}-\begin{pmatrix}1 & 2/5 &1/5\\ 5/7 & 1 &1/7\\ 3/6 &2/6 &1/6 \end{pmatrix}\Biggr]\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}+\begin{pmatrix}1/5\\2/7\\3/6\end{pmatrix}\]

\[\begin{pmatrix}1 & 0 &0\\ 5/7 & 1 & 0\\ 3/6 &2/6 &1\end{pmatrix}\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}=\begin{pmatrix}0& -0.4 &-0.2\\ 0 & 0 & -0.142857\\ 0 & 0 &0\end{pmatrix}\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}+\begin{pmatrix}1/5\\2/7\\3/6\end{pmatrix}\]

function gausssiedel(n,A,B,X,M,Eps)
    # dla i od 1 do n
     #     d=1/A_{i,i}
     #     B_i=d*B_i;
     #     dla j od 1 do n
     #         A_{i,j}=d*A_{i,j};
     # dla k od 1 do M
     #     dla i od 1 do n
     #         X_i = B_i-sum(A_{ij}*X_j dla j od 1 do n jeśli j różne od i)
     #     jeśli sqrt(sum((X-U).^2))<Eps) zwróc X
     #     dla i od 1 do n
     #         U_i = X_i
     # zwróć "Algorithm failed to converge"
end

Jeśli macierz współczynników układu równań liniowych A

• jest rzeczywista
• kwadratowa nxn
• symetryczna \(A^T=A\)
• dodatnio określona \(x^TAx>0,x\neq 0\)

to problem znajdowania x spełniającego równanie Ax=b jest równoważny ze znajdowaniem x, dla którego funkcja kwadratowa \(q(x)=\langle x,Ax\rangle-2\langle x,b\rangle\) osiąga minimum. Ten związek można udowodnić obliczając gradient funkcji g(x): \(\nabla q(x) = 2Ax - 2b\). Minimum funkcji \(q(x)\) występuje gdy \(\nabla q(x) = 0\), czyli: 2Ax - 2b = 0 –> Ax=b.

Funkcja q(x) zawsze ma minimum jedno miejsce zerowe – w punkcie 0. Jeśli \(b=0\), to funkcja jest dodatnia dla każdego x, gdyż macierz A jest dodatnio określona. Gdy \(b \neq 0\), to funkcja przyjmuje wartości ujemne dla x między miejscami zerowymi. Funkcja q(x) z macierzą dodatnio określoną ma minimum globalne.

Funkcję kwadratową zdefiniowaliśmy, wykorzystując iloczyn skalarny dwóch wektorów x i y: \(\langle x,y\rangle = \sum_{i=1}^n x_iy_i\)

Własności iloczynu:

1. Przemmienność: \(\langle x,y\rangle=\langle y,x\rangle\)
2. Jednorodność: \(\langle ax,y\rangle=a\langle x,y\rangle,a\in \mathbb{R}\)
3. Rozłączność: \(\langle x+y,z\rangle=\langle x,z\rangle + \langle y,z\rangle\)
4. Dla operacji macierzowych: \(\langle x,Ay\rangle=\langle A^Tx,y\rangle\)

Przykład dla A=5 i b=2
\(5x-2=0\) -> \(x=\frac{2}{5}\)
\(q(x)=\langle x,5x\rangle-2\langle x,2\rangle=5x^2-4x\)

Funkcja kwadratowa osiąga minimum, gdy jej pochodna jest równa 0. \(g'(x)=10x-4=0\) –> \(x=\frac{2}{5}\), to rozwiązanie jest równe rozwiązaniu równania liniowego \(5x-b=0\)

Wyjaśnienie graficzne: https://www.geogebra.org/m/uqgsmvye

Przykład dla A=[5 4;4 5] i b=[2;2]
\[\begin{pmatrix}5& 4 \\ 4 & 5\end{pmatrix}\begin{pmatrix}x_1\\x_2\end{pmatrix}=\begin{pmatrix}2\\2\end{pmatrix}\] \[5x_1+4x_2=2, 4x_1+5x_2=2\] \[x_2=\frac{1}{2}-\frac{5}{4}x_1, 4x_1+\frac{5}{2}-\frac{25}{4}x_1=2\] \[\frac{9}{4}x_1=\frac{1}{2},x_1=\frac{2}{9},x_2=\frac{1}{2}-\frac{5}{18}=\frac{2}{9}\]

\[q(x)=\langle \begin{pmatrix}x_1\\x_2\end{pmatrix},\begin{pmatrix}5& 4 \\ 4 & 5\end{pmatrix}\begin{pmatrix}x_1\\x_2\end{pmatrix}\rangle-2\langle \begin{pmatrix}x_1\\x_2\end{pmatrix},\begin{pmatrix}2\\2\end{pmatrix}\rangle\] \[q(x)=\langle \begin{pmatrix}x_1\\x_2\end{pmatrix},\begin{pmatrix}5x_1+4x_2 \\ 4x_1 +5x_2\end{pmatrix}\rangle-2\langle \begin{pmatrix}x_1\\x_2\end{pmatrix},\begin{pmatrix}2\\2\end{pmatrix}\rangle\] \[q(x)=5x_1^2+4x_1x_2+5x_2^2+4x_1x_2 - 4x_1-4x_2\] \[q(x)=5x_1^2+5x_2^2+8x_1x_2- 4x_1-4x_2\] \[\frac{dq}{dx_1}=10x_1+8x_2- 4=5x_1+4x_2-2=0\] \[\frac{dq}{dx_2}=10x_1+8x_2- 4=5x_1+4x_2-2=0\]

Wyjaśnienie graficzne: https://www.geogebra.org/3d/fbfgpezb

Jak znaleźć rozwiązanie numerycznie bez liczenia pochodnej względem wektora X? Znajdźmy jednowymiarową krzywą wzdłuż badanej powierzchni funkcji. Nasza krzywa rozpoczyna się w x i zmienia się zgodnie z wektorem v, który przemnażamy z każdym krokiem przez t.

\(x^{k+1}=x^{k}+t_kv^{k}\)

Metoda najszybszego spadku bazuje na prostej intuicji: iteracyjnie poruszamy się w kierunku przeciwnym do gradientu funkcji q(x), czyli w kierunku najszybszego spadku wartości funkcji. Jest to analogiczne do schodzenia z góry po najstromszej ścieżce. Kierunek spadku v jest jest zgodny z gradientem funkcji q(x): \(v^k = −(2Ax^k - 2b) = 2(b-Ax^k)\). Dla uproszczenia przyjmujemy \(v^k = b-Ax^k = r^k\), gdzie \(r^k\) to wektor reszt (residual), reprezentujący "błąd" w przybliżeniu równania \(Ax=b\). Pozostaje wyznaczyć optymalną wielkość kroku t.

\[q(x+tv) = \langle x+tv, A(x+tv)\rangle-2\langle x+tv,b\rangle\] \[q(x+tv) = \langle x+tv, Ax+Atv)\rangle-2\langle x+tv,b\rangle\] \[q(x+tv) = \langle x,Ax+Atv\rangle+\langle tv,Ax+Atv\rangle-2\langle x+tv,b\rangle\] \[q(x+tv) = \langle Ax,x\rangle+\langle Ax,tv\rangle+\langle tv,Ax\rangle+\langle tv,Atv\rangle-2\langle x+tv,b\rangle\] \[q(x+tv) = \langle Ax,x\rangle+t\langle Ax,v\rangle+t\langle v,Ax\rangle+t^2\langle v,Av\rangle-2\langle x,b\rangle-2t\langle v,b\rangle\] \[q(x+tv) = \langle Ax,x\rangle-2\langle x,b\rangle+2t\langle Ax,v\rangle-2t\langle b,v\rangle+t^2\langle v,Av\rangle\] \[q(x+tv) = \langle Ax,x\rangle-2\langle x,b\rangle+2t\langle Ax-b,v\rangle+t^2\langle v,Av\rangle\] \[q(x+tv) = q(x)+2t\langle Ax-b,v\rangle+t^2\langle v,Av\rangle\]

Ta funkcja względem t jest funkcją kwadratową, tj. ma jedno minimum, a nie ma maksimum. Wartość parametru t, dla którego pochodna tej funkcji równa się 0 wykorzystamy do znalezienia wartości minimalnej naszej funkcji.

\[\frac{q(x+tv)}{dt}=2\langle Ax-b,v\rangle+2t\langle Av,v\rangle=0\] \[t^*=\frac{\langle b-Ax,v\rangle}{\langle Av,v\rangle}\] \[t^*=\frac{\langle v,v\rangle}{\langle Av,v\rangle}\] \[t^*=\frac{||v||^2}{\langle Av,v\rangle}\]

Dla t^* funkcja q(x+tv) osiąga minimum:
\[q(x+t^*v) = q(x)+t^*[2\langle Ax-b,v\rangle+t^*\langle Av,v\rangle]\] \[q(x+t^*v) = q(x)+t^*[2\langle Ax-b,v\rangle+\langle b-Ax,v\rangle]\] \[q(x+t^*v) = q(x)-t^*\langle b-Ax,v\rangle\] \[q(x+t^*v) = q(x)-\frac{\langle b-Ax,v\rangle^2}{\langle Av,v\rangle}\]

Wyrażenie ułamkowe zawsze będzie większe od zera, czyli q(x+t^*v) będzie mniejsze od q(x).

Sprawdźmy, czy dla naszych wcześniej omawiany przykładów znajdziemy wartość minimalną funkcji.

Dla A=5 i b=2 z warunkiem początkowym x=1:
\(q(1)=5\cdot 1^2-4\cdot 1=1\)
\(q(1+t^*v)=1-\frac{\langle 2-5,v\rangle^2}{\langle 5v,v\rangle}\)
\(q(1+t^*v)=1-\frac{9v^2}{5v^2}\)
\(q(1+t^*v)=1-\frac{9}{5}=-\frac{4}{5}\)

Podstawiając tą wartość do funkcji, którą wykorzystujmy do znalezienia miejsca zerowego układów równań liniowych, to otrzymamy: \(-\frac{4}{5}=5x^2-4x\), wiemy, że miejsce zerowe naszej funkcji liniowej jest dla x=2/5; gdy podstawimy tą wartość do naszego równania, to otrzymamy: \(-\frac{4}{5}=5\cdot \frac{4}{25}-4\frac{2}{5}=\frac{4}{5}-\frac{8}{5}=-\frac{4}{5}\), wszystko się zgadza.

Drugi przykład dla A=[5 4;4 5] i b=[2;2] z warunkiem początkowych x=[1;1] \[q(x)=5x_1^2+5x_2^2+8x_1x_2- 4x_1-4x_2\] \[q([1;1])=5+5+8-4-4=10\] \[q([1;1]+t^*v)=10-\frac{\langle \begin{pmatrix}2\\2\end{pmatrix}-\begin{pmatrix}5& 4 \\ 4 & 5\end{pmatrix}\begin{pmatrix}1\\1\end{pmatrix},\begin{pmatrix}v_1\\v_2\end{pmatrix}\rangle^2}{\langle \begin{pmatrix}5& 4 \\ 4 & 5\end{pmatrix}\begin{pmatrix}v_1\\v_2\end{pmatrix},\begin{pmatrix}v_1\\v_2\end{pmatrix}\rangle}\] \[q([1;1]+t^*v)=10-\frac{\langle \begin{pmatrix}-7\\-7\end{pmatrix},\begin{pmatrix}v_1\\v_2\end{pmatrix}\rangle^2}{\langle \begin{pmatrix}5& 4 \\ 4 & 5\end{pmatrix}\begin{pmatrix}v_1\\v_2\end{pmatrix},\begin{pmatrix}v_1\\v_2\end{pmatrix}\rangle}\]

Wektor v określa kierunek iteracji, więc powinien to być gradient spadku (gradient wzrostu $∇ q(x)), który jest równy reszcie r: \[r=b-Ax\]. To jest nasza funkcja kosztu: sprawdzamy jak daleko nasze obliczenia Ax są od oczekiwanych b. Kierunek wektor r jest zgodny z gradientem wzrostu, ale ma przeciwny zwrot. W naszym przykładzie: \[v=r=\begin{pmatrix}2\\2\end{pmatrix}-\begin{pmatrix}5& 4 \\ 4 & 5\end{pmatrix}\begin{pmatrix}1\\1\end{pmatrix}=\begin{pmatrix}-7\\-7\end{pmatrix}\rangle\] \[\frac{dq(x)}{x_1}(x_1=1,x_2=1)=14=-2(b-Ax)\] \[\frac{dq(x)}{x_2}(x_1=1,x_2=1)=14=-2(b-Ax)\]

\[q([1;1]+t^*v)=10-\frac{\langle \begin{pmatrix}-7\\-7\end{pmatrix},\begin{pmatrix}-7\\-7\end{pmatrix}\rangle^2}{\langle \begin{pmatrix}5& 4 \\ 4 & 5\end{pmatrix}\begin{pmatrix}-7\\-7\end{pmatrix},\begin{pmatrix}-7\\-7\end{pmatrix}\rangle}\] \[q([1;1]+t^*v)=10-\frac{9604}{882}=-\frac{8}{9}\]

Podstawiamy do q(x):
\(-\frac{8}{9}=5\cdot \frac{2}{9}^2+5\cdot \frac{2}{9}^2+8\cdot \frac{2}{9}\cdot \frac{2}{9}- 4\cdot \frac{2}{9}-4\cdot \frac{2}{9}\), \(-\frac{8}{9}=-\frac{8}{9}\), wszystko się zgadza.

Wzory na zbudowanie algorytmu:
\(v^{k}=b-Ax^{k}\)
\(t_k=\frac{\langle v^{k},v^{k}\rangle}{\langle Av^{k},v^{k}\rangle}\)
\(x^{k+1}=x^{k}+t_kv^{k}\)

Gdy wektor r unormujemy \(\lVert {v^{k}}\rVert = 1\), to t_k zmierzymy odległość pomiędzy wektorami x^k i x^k+1.

Szybkość zbieżności zależy od uwarunkowania macierzy A (stosunku największej do najmniejszej wartości własnej). Dla źle uwarunkowanych macierzy, metoda najszybszego spadku może wymagać wielu iteracji, gdyż tworzy charakterystyczny wzór "zig-zag"

Mamy dane uczące: \({(x_{i},y_{i})}^n_{i=1}\), gdzie:

• i – to indeks pojedynczej obserwacji (od 1 do n)
• \(x_{i}\) – to wektor cech dla i-tej obserwacji np dla dwóch cech (x_i1, x_i2); wielkość tego wektora to d+1, gdzie d to wymiar wielomianu; dla regresji liniowej d=1
• \(y_{i}\) – to wartość rzeczywista (odpowiedź) dla tej obserwacji.

chcemy dopasować model liniowy: \(\hat{y} = w^T \phi(x_{i})\), gdzie:

• \(w \in \mathbb{R}^{d+1}\) – to wektor parametrów (współczynników regresji wraz z elementem wolnym),
• \(\phi(x_{i}) \in \mathbb{R}^{d+1}\) – przekształcenie wejścia x_i na wektor cech (wartości liczbowych używanych przez model). Dla regresji liniowej \(x_{i}=\phi(x_{i})\)
  • dla regresji liniowej z wyrazem wolnym \(\phi(x_{i}) = \begin{pmatrix} 1 \\ x_{i} \end{pmatrix}\)
  • dla regresji wielomianowej stopnia 2: \(\phi(x_{i}) = \begin{pmatrix} 1 \\ x_{i} \\ x_{i}^2 \end{pmatrix}\)
• \(\hat{y}\) – to estymowana wartość dla x_i

Minimalizujemy:

• błąd średniokwadratowy (MSE): \(\text{J}(w)=\text{MSE}(w) = \frac{1}{n} \sum_{i=1}^n \left( \hat{y}_i - y_i \right)^2 = \frac{1}{n} \sum_{i=1}^n \left( w^T \varphi(x_i) - y_i \right)^2\); J(w) – to ogólne oznaczenie funkcji kosztu (ang. cost function lub loss function) w optymalizacji.

lub:

• \(min_{w} \sum_{i=1}^n \left( w^T \varphi(x_i) - y_i \right)^2\)

Dla regresji liniowej:

• \(w \in \mathbb{R}^{2}\) np. \(\alpha\) i \(\beta\)
• \(\phi(x_{i}) \in \mathbb{R}^{2}\), \(\phi(x_{i}) = \begin{pmatrix} 1 \\ x_{i1} \end{pmatrix}\)

\(min_{\alpha, \beta} \sum_{i=1}^n \left( \alpha + \beta x_i - y_i \right)^2\)

\(\sum_{i=1}^n (\alpha + \beta x_i) - y_i)^2 = \left\lVert \begin{pmatrix} \alpha + \beta x_1) - y_1 \\ . \\ . \\ . \\ \alpha + \beta x_n) - y_n \end{pmatrix} \right\rVert_2^2 = \left\lVert \begin{pmatrix} 1 & x_{1} \\ . \\ . \\ . \\ 1 & x_{n} \end{pmatrix} \begin{pmatrix} \alpha \\ \beta \end{pmatrix} - \begin{pmatrix} y_1 \\ . \\ . \\ . \\ y_n \end{pmatrix} \right\rVert_2^2 = \left\lVert Ax - b \right\rVert_2^2\)

gdzie:

\(A = \begin{pmatrix} 1 & x_{1} \\ . \\ . \\ . \\ 1 & x_{n} \end{pmatrix}\), \(x = \begin{pmatrix} \alpha \\ \beta \end{pmatrix}\); \(b = \begin{pmatrix} y_1 \\ . \\ . \\ . \\ y_n \end{pmatrix}\)

to nasz problem znajdowania krzywej regresji przyjmuje postać: \(min_{x} \left\lVert Ax - b \right\rVert_2^2\), \(A \in \mathbb{R}^{n \times (d+1)}\), \(x \in \mathbb{R}^{d+1}\), \(b \in \mathbb{R}^n\), \(n \ge (d+1)\).

Przejście do postaci kwadratowej: q(x)=⟨x,Ax⟩−2⟨x,b⟩ – związanie naszego problemu ze znajdowaniem globalnego minimum funkcji.

\(q(x) = \left\lVert Ax - b \right\rVert_2^2 = (Ax - b)^{T}(Ax - b)\)

Wzór skróconego mnożenia: \((u - v)^T (u - v) = u^T u - u^T v - v^T u + v^T v\)

Własność transpozycji: \((A x)^T = x^T A^T\)

\(q(x)=(Ax - b)^{T}(Ax - b) = (Ax)^T Ax - (Ax)^T b - b^T Ax + b^T b\) \(q(x)=x^T A^T Ax - x^T A^T b - b^T Ax + b^T b\)

Ponadto: \(b^T Ax\) to skalar, a skalar jest równy swojej transpozycji. \(b^T Ax = (b^T Ax)^T = x^TA^Tb\)

\(q(x)=x^T A^T Ax - x^T A^T b - x^T A^T b + b^T b\) \(q(x)=x^T A^T Ax - 2x^T A^T b + b^T b\)

\(\nabla q(x) = 2 A^T A x - 2 A^T b = 2A^T(Ax-b)\)

\(min_{x} \left\lVert Ax - b \right\rVert_2^2 \Leftrightarrow \nabla q(x) = 0\)

\(2A^T(Ax-b)=0\)

\(x=(A^TA)^{-1} (A^T b)\)

Wielomian możemy przedstawić w postaci liniowej kombinacji wielomianów 1,x,x² … np.
\(w(x)=a_{0}+a_{1}x+a_{2}x^2+\text{ ... }+a_{n}x^n\)

Dla warunku interpolacyjnego \(w(x_{i})=y_{i},\text{ gdzie } 0\leq i \leq n\) taki wielomian następnie możemy przedstawić w postaci macierzowej z n+1 równaniami liniowymi oraz n+1 niewiadomymi np. dla wielomianu stopnia n=2 mamy:

\(\begin{pmatrix}1 & x_{0} &x_{0}^2\\ 1 & x_{1} & x_{1}^2\\ 1 &x_{2}&x_{2}^2\end{pmatrix}\begin{pmatrix}a_0\\a_1\\a_2\end{pmatrix}=\begin{pmatrix}y_0\\y_1\\y_2\end{pmatrix}\)

\(w(x)=y_{0}l_{0}(x)+y_{1}l_{1}(x)+\text{ ... }+y_{n}l_{n}(x)=\sum\limits_{k=0}^n y_{k}l_{k}(x)\)

\(l_{i}=\prod\limits_{j=0;j\neq i}^n \frac{x-x_{j}}{x_{i}-x_{j}},\,0 \leq i \leq n\)

l_i - wielomian, który zależy od każdej wartości x_i, a nie zależy od wartości interpolowanej funkcji. Wielomian interpolacyjny Langrangea jest sumą iloczynów wartości interpolowanej funkcji oraz wielomianów l_i. Inaczej powstaje wielomian interpolacyjny Newtona, omawiany w następny rozdziale – w każdym kroku iteracji powstający tam wielomian interpolacyjny Newtona zależy od wartości funkcji interpolacyjnej oraz od kolejnych wartości x_i.

Przykład dla \(f(x)=x^3+5x^2-2x-5\) i wybranych punktów interpolacyjnych x₀=2,x₁=-1,x₂=-4,x₃=1
\(l_{0}=\frac{x-x_{1}}{x_{0}-x_{1}}\frac{x-x_{2}}{x_{0}-x_{2}}\frac{x-x_{3}}{x_{0}-x_{3}}=\frac{(x+1)(x+4)(x-1)}{(2+1)(2+4)(2-1)}=\frac{1}{18}(x+1)(x+4)(x-1)\)
\(l_{1}=\frac{x-x_{0}}{x_{1}-x_{0}}\frac{x-x_{2}}{x_{1}-x_{2}}\frac{x-x_{3}}{x_{1}-x_{3}}=\frac{(x-2)(x+4)(x-1)}{(-1-2)(-1+4)(-1-1)}=\frac{1}{18}(x-2)(x+4)(x-1)\)
\(l_{2}=\frac{x-x_{0}}{x_{2}-x_{0}}\frac{x-x_{1}}{x_{2}-x_{1}}\frac{x-x_{3}}{x_{2}-x_{3}}=\frac{(x-2)(x+1)(x-1)}{(-4-2)(-4+1)(-4-1)}=-\frac{1}{90}(x-2)(x+1)(x-1)\)
\(l_{3}=\frac{x-x_{0}}{x_{3}-x_{0}}\frac{x-x_{1}}{x_{3}-x_{1}}\frac{x-x_{2}}{x_{3}-x_{2}}=\frac{(x-2)(x+1)(x+4)}{(1-2)(1+1)(1+4)}=-\frac{1}{10}(x-2)(x+1)(x+4)\)

\(y_{0} = f(2)=19; y_{1}=f(-1) = 1; y_{2}=f(-4) = 19;y_{3}=f(1) = -1\)

Wielomian interpolacyjny Langrangea dla badanej funkcji ma postać:
\(w(x)=19l_{0}(x)+l_{1}(x)+19l_{2}(x)-l_{3}(x)\)

Jeśli w wielomianie Langrangea wykonamy podstawienie x=x_i, to otrzymamy postać macierzową:
\(\begin{pmatrix}l_{0}(x_{0}) & l_{1}(x_{0}) &l_{2}(x_{0}) &l_{3}(x_{0})\\ l_{0}(x_{1}) & l_{1}(x_{1}) &l_{2}(x_{1}) &l_{3}(x_{1})\\ l_{0}(x_{2}) & l_{1}(x_{2}) &l_{2}(x_{2}) &l_{3}(x_{2})\\ l_{0}(x_{3}) & l_{1}(x_{3}) &l_{2}(x_{3}) &l_{3}(x_{3})\end{pmatrix}\begin{pmatrix}y_0\\y_1\\y_2\\y_3 \end{pmatrix}=\begin{pmatrix}w_{0}\\w_{1}\\w_{2}\\w_{3}\end{pmatrix}\)

W tym przypadku macierz współczynników jest macierzą jednostkową, gdyż wielomiany znajdujące się poza przekątną, zerują się (sprawdź wzory na l_i).

\(w_{0} = 19l_0=\frac{19}{18}(x+1)(x+4)(x-1)\)
\(w_{1} = 1l_1=\frac{1}{18}(x-2)(x+4)(x-1)\)
\(w_{2} = 19l_2=-\frac{19}{90}(x-2)(x+1)(x-1)\)
\(w_{3} = -1l_3=\frac{1}{10}(x-2)(x+1)(x+4)\)

\(w(x)=\sum\limits_{i=0}^n w_{i}\)

\(0:\,w_{0}=a_{0}\)
\(1:\,w_{1}=a_{0}+a_{1}(x-x_{0})\)
\(2:\,w_{2}=a_{0}+a_{1}(x-x_{0})+a_{2}(x-x_{0})(x-x_{1})\)
\(3:\,w_{3}=a_{0}+a_{1}(x-x_{0})+a_{2}(x-x_{0})(x-x_{1})+a_{3}(x-x_{0})(x-x_{1})(x-x_{2})\)
\(n:\,w_{n}=a_{0}+a_{1}(x-x_{0})+a_{2}(x-x_{0})(x-x_{1})+\text{ ... }+a_{n}(x-x_{0})\text{ ... }(x-x_{n-1})\)

Postać uproszczona: \(w_{n}=\sum\limits_{i=0}^n a_{i}\prod\limits_{j=0}^{i-1}(x-x_{j})\)

Jeśli wykonamy podstawienie x=x_i, to otrzymamy postać macierzową dolno-trójkątną: \(\begin{pmatrix}1 & 0 & 0\\ 1 & (x_{1}-x_{0}) & 0\\ 1 & (x_{2}-x_{0})& (x_{2}-x_{0})(x_{2}-x_{1})\end{pmatrix}\begin{pmatrix}a_0\\a_1\\a_2\end{pmatrix}=\begin{pmatrix}f(x_0)\\f(x_1)\\f(x_2)\end{pmatrix}\)

\(a_{0} = f(x_{0})\)
\(a_{1} = \frac{f(x_{1})-a_{0}}{x_{1}-x_{0}}=\frac{f(x_{1})-f(x_{0})}{x_{1}-x_{0}}\)

\(a_{2} = \frac{f(x_{2})-a_{0}-(x_{2}-x_{0})a_{1}}{(x_{2}-x_{0})(x_{2}-x_{1})}= \frac{f(x_{2})-f(x_{0})-(x_{2}-x_{0})\frac{f(x_1)-f(x_0)}{x_1-x_0}}{(x_{2}-x_{0})(x_{2}-x_{1})}= \frac{(x_1-x_0)(f(x_{2})-f(x_{0}))-(x_{2}-x_{0})(f(x_1)-f(x_0))}{(x_{2}-x_{0})(x_{2}-x_{1})(x_1-x_0)}\)

Do obliczania współczynników a będziemy korzystać z poniższych wzorów na ilorazy różnicowe:
\(a_0=f[x_{0}]=f(x_{0})\)
\(a_1=f[x_0,x_1] = \frac{f(x_{1})-f(x_{0})}{x_{1}-x_{0}}\)
\(a_2=f[x_0,x_1,x_2]=\frac{f[x_1,x_2]-f[x_0,x_1]}{x_2-x_0}\)

\(\frac{f[x_1,x_2]-f[x_0,x_1]}{x_2-x_0}= \frac{\frac{f(x_2)-f(x_1)}{x_2-x_1}-\frac{f(x_1)-f(x_0)}{x_1-x_0}}{x_2-x_0}= \frac{(x_1-x_0)(f(x_{2})-f(x_{0}))-(x_{2}-x_{0})(f(x_1)-f(x_0))}{(x_{2}-x_{0})(x_{2}-x_{1})(x_1-x_0)}\)

Wzór ogólny na ilorazy różnicowe:
\(f[x_i,x_{i+1},...,x_{i+j}]=\frac{f[x_{i+1},x_{i+2},...,x_{i+j}]-f[x_{i},x_{i+1},...,x_{i+j-1}]}{x_{i+j}-x_i}\)

Kolejność wykonywania operacji od lewej do prawej:

Szukany wielomian jest liniową kombinacją współczynników z pierwszego wiersza.

Przykład dla \(f(x)=x^3+5x^2-2x-5\) i wybranych punktów interpolacyjnych x₀=2,x₁=-1,x₂=-4,x₃=1

Warunki początkowe:

Obliczenia:

w₀ = 19
w₁ = 19 + 6(x-2)
w₂ = 19 + 6(x-2)+2(x-2)(x+1)
w₃ = 19 + 6(x-2)+2(x-2)(x+1)+(x-2)(x+1)(x+4)

19 + 6(x-2)+2(x-2)(x+1)+(x-2)(x+1)(x+4) = \(x^3+5x^2-2x-5\) - to musi być tym samym, gdyż te wielomiany stopnia trzeciego są takie same w punktach \(w(x_{i})=y_{i}\). A dla tych punktów istnieje tylko jeden wielomian stopnia trzeciego.

Algorytm do obliczania współczynników wielomianu, gdzie A jest macierzą odpowiadającą tabeli pokazanej wyżej, tj. zawierającą wszystkie ilorazy różnicowe oraz wartości początkowe x_i.

A # macierz Nx(N+1)
A[:,1]=# wartości początkowe X
A[:,2]=# wartości badanej funkcji f(x)
dla j = 3:(n+1)
        dla i = 1:(n+2-j)
            A[{i}{j}] = (A[{i+1}{j-1}] - A[{i}{j-1}])/(A[{i+j-2}{1}] - A[{i}{1}])

\(\int\limits_{a}^b f(x)dx \approx f(\frac{a+b}{2})(b-a)\)

Nasza funkcja g jest wielomianem stopnia 0, czyli jest wartością stałą. Geometrycznie reprezentuje pole powierzchni prostokąta.

Źródło: Wikipedia

Funkcję podcałkową przybliżamy wielomianem stopnia pierwszego, czyli funkcją liniową. Poprzez takie przybliżenie, pole powierzchni, które tworzy funkcja g, jest polem trapezu:
Źródło: Wikipedia

To pole jest równe:
\(f(x_i)\Delta x + \frac{1}{2}(f(x_{i+1})-f(x_{i}))\Delta x = \frac{f(x_i) + f(x_{i+1})}{2}\Delta x\)

\(\int\limits_{a}^b f(x)dx \approx \frac{f(x_i) + f(x_{i+1})}{2}\Delta x\)

\(\int\limits_{a}^b f(x)dx \approx \frac{b-a}{6}[f(a)+4f(\frac{a+b}{2})+f(b)]\)

Funkcję podcałkową przybliżamy wielomianem stopnia drugiego, czyli funkcją kwadratową. Poprzez takie przybliżenie, pole powierzchni, które tworzy funkcja g jest polem pod parabolą:
Źródło: Wikipedia