Dany jest układ równań nieliniowych z \(\displaystyle{ n - }\) niewiadomymi:
\(\displaystyle{ \begin{cases} f_{1}(x_{1}, x_{2},...,x_{n}) = 0 \\ f_{2}(x_{1}, x_{2},...,x_{n}) = 0 \\.....................\ \ (1)\\ f_{n}(x_{1}, x_{2},...,x_{n}) = 0, \end{cases} }\)
który możemy zapisać w postaci wektorowej
\(\displaystyle{ \vec{f}(\vec{x}) = \vec{0} }\)
gdzie:
\(\displaystyle{ \vec{f}(\vec{x}) = \left[ \begin{matrix} f_{1}(\vec{x})\\ f_{2}(\vec{x}) \\ .....\\ f_{n}(\vec{x})\end{matrix} \right], \ \ \vec{x} = \left[ \begin{matrix} x_{1} \\ x_{2}\\ ...\\ x_{n} \end{matrix} \right]^{T}. }\)
O funkcjach \(\displaystyle{ f_{i}(x_{1},x_{2},...,x_{n}), \ \ i=1,2,...,n }\) zakładamy, że mają ciągłe pochodne pierwszeg rzędu w pewnym obszarze
zawierającym pierwiastek układu równań \(\displaystyle{ (1). }\)
Niech \(\displaystyle{ \vec{x}^{(k)} =(x^{(k)}_{1}, x^{(k)}_{2}, ..., x^{(k)}_{n}) }\) będzie \(\displaystyle{ k - }\) tym przybliżeniem pierwiastka \(\displaystyle{
\vec{\alpha} = (\alpha_{1}, \alpha_{2},...,\alpha_{n}). }\)
Dokładną wartość pierwiastka \(\displaystyle{ \vec{\alpha} }\) można wyrazić w postaci równania:
\(\displaystyle{ \vec{\alpha} =\vec{x}^{(k)} + \vec{\epsilon}^{(k)} .}\)
gdzie:
\(\displaystyle{ \vec{\epsilon}^{(k)} = (\varepsilon^{(k)}_{1}, \varepsilon^{(k)}_{2},...,\varepsilon^{(k)}_{n}) }\)
jest błędem przybliżenia pierwiastka \(\displaystyle{ vec{x}^{(k)} }\) w \(\displaystyle{ k - }\) tej iteracji.
Wobec istnienia ciągłej pochodnej funkcji \(\displaystyle{ \vec{f}(\vec{x}) }\) możemy napisać
\(\displaystyle{ \vec{f}(\vec{\alpha}) \approx \vec{f}^{(k)}(\vec{x}^{(k)}) + \vec{f'}^{(k)}(\vec{x}^{(k)})\cdot \vec{\epsilon}^{(k)} \ \ (1)}\)
Zastępując błąd \(\displaystyle{ \vec{\epsilon}^{(k)} }\) przyrostem \(\displaystyle{ \vec{h}^{(k)} = \vec{x}_{k+1} - \vec{x}_{k} }\) i porównując prawą stronę
równania \(\displaystyle{ (1) }\) do zera, otrzymujemy równanie liniowe
\(\displaystyle{ \vec{f}(\vec{x}^{(k)}) +\cdot (\vec{x}^{(k+1)} - \vec{x}^{(k)})= \vec{0} \ \ (2)}\)
Wzór \(\displaystyle{ (2) }\) jest wzorem rekurencyjnym dla Wielowymiarowej Metody Newtona w postaci macierzowej.
Pochodna \(\displaystyle{ \vec{f'}(\vec{x})}\) jest macierzą Jacobiego.
Zakadając jej nieosobliwość, możemy równanie \(\displaystyle{ (2) }\) napisać w postaci:
\(\displaystyle{ \vec{x}^{(k+1)} = \vec{x}^{(k)} + \vec{f^{-1}}(\vec{x}^{(k)})\cdot \vec{f}(\vec{x}^{(k)}).}\)
Przykład
Rozwiązujemy układ równań:
\(\displaystyle{ \begin{cases} f_{1}(x_{1},x_{2}) = x^3_{1} +3x^2_{2} -21 = 0 \\f_{2}(x_{1},x_{2}) = x^2_{1} +2x_{2} +2 = 0 \end{cases} }\)
Jakobian macierzy układu
\(\displaystyle{ \vec{f'}(x_{1},x_{2}) = \left[ \begin{matrix} 3x^2_{1} & 6x_{2} \\ 2x_{1} & 2 \end{matrix} \right] }\)
Przyjmujemy wektor początkowy \(\displaystyle{ \vec{x}^{(0)} = [1, -1] }\)
Obliczamy dla tego wektora wartość wektor funkcji układu i wartość jej Jakobianu.
\(\displaystyle{ \vec{f}(1,-1) = \left[ \begin{matrix} -17 \\ 1 \end{matrix} \right], \ \ \vec{f'}(1,-1)= \left[ \begin{matrix} 3 & -6 \\ 2 & 2 \end{matrix} \right].}\)
Rozwiązujemy układ równań:
\(\displaystyle{ \left[ \begin{matrix} 3 & -6 \\ 2 & 2 \end{matrix} \right]\cdot \left[ \begin{matrix} h^{(0)}_{1} \\ h^{(0)}_{2} \end{matrix} \right] = - \left[\begin{matrix} -17 \\1 \end{matrix} \right]. }\)
Otrzymujemy
\(\displaystyle{ \vec{h}^{(0)} = \left[ \begin{matrix} 1.555556\\ -2.055560 \end{matrix} \right] }\)
Znajdujemy wektor w następnej iteracji:
\(\displaystyle{ \vec{x}^{(1)} = \left [ \begin{matrix} x^{(0)}_{1} \\ x^{(0)}_{2} \end{matrix}\right] = \left[ \begin{matrix} 1 \\ -1 \end{matrix} \right] + \left[ \begin{matrix} 1.555556\\ -2.055560 \end{matrix} \right] = \left[ \begin{matrix} 2.555556 \\ -3.055556 \end{matrix} \right] }\)
Podobnie znajdujemy pozostałe wektory.
Program w OCTAVE
Kod: Zaznacz cały
function newton_sys(F,JF,xx0,tol,maxit)
x0=xx0;
iter=1;
while(iter<=maxit)
y=feval(JF,x0)\feval(F,x0);
xn=x0+y;
err=max(abs(xn-x0));
if (err<=tol)
x=xn;
fprintf('Metoda Newtona zbieżna po %3.0f iteracjach do \n',iter)
x
return;
else
x0=xn;
end
iter=iter+1;
endwhile
disp('Metoda Newtona nie jest zbieżna')
x=xn;
function f=f1(X)
x=X(1);y=X(2);
f=[x.^3+3*y.^2-21;x.^2+2*y+2];
function df=df1(X)
x=X(1);y=X(2);
df=[3*x.^2,6*y;2*x,2];
>> newton_sys('f1','df1',[1,-1]',10^(8),40)
Metoda Newtona zbieżna po 6 iteracjach do
x =
1.643038
-2.349787