Antes de analizar cómo resolver un sistema de ecuaciones multivariables, es útil revisar el desarrollo en serie de Taylor de una función de N dimensiones.
El método de Newton-Raphson para una variable, \(f(x)=0\), puede generalizarse de forma similar.
\(\left\{\begin{matrix}
f_1\left ( x_1,\cdots ,x_n \right )=f_1\left ( \mathbf{x} \right )
\\ f_2\left ( x_1,\cdots ,x_n \right )=f_2\left ( \mathbf{x} \right )
\\ \cdots \\
f_n\left ( x_1,\cdots ,x_n \right )=f_n\left ( \mathbf{x} \right )
\end{matrix}\right.\)
donde,
\(\mathbf{x}=\left [ x_1,\cdots ,x_n \right ]^{T}\)
\(\mathbf{f\left (x\right )}=\left [ f_1\left (\mathbf{x}\right ),\cdots ,f_n\left (\mathbf{x}\right ) \right ]^{T}\)
\(\mathbf{f\left (x\right )}=0\)
Mediante el desarrollo en serie de Taylor,
\(\begin{matrix}
f_i\left ( \mathbf{x}+\delta \mathbf{x} \right ) = f_i\left ( \mathbf{x} \right )+\displaystyle\sum_{j=1}^{N}\frac{\partial f_i\left ( \mathbf{x} \right )}{\partial x_j}\delta x_j+O\left ( \partial x^{2} \right ) \\
\approx f_i\left ( \mathbf{x} \right )+\displaystyle\sum_{j=1}^{N}\frac{\partial f_i\left ( \mathbf{x} \right )}{\partial x_j}\delta x_j
\end{matrix}\)
Escribiendo las \(N\) ecuaciones en forma vectorial,
\(\begin{matrix}
f\left ( \mathbf{x}+\delta \mathbf{x} \right ) =\displaystyle\begin{bmatrix}
f_1\left ( \mathbf{x}+\delta \mathbf{x} \right ) \\ \vdots \\ f_N\left ( \mathbf{x}+\delta \mathbf{x} \right )
\end{bmatrix} \\
\approx \displaystyle\begin{bmatrix}
f_1\left ( \mathbf{x} \right ) \\
\cdots \\
f_n\left ( \mathbf{x} \right )
\end{bmatrix}+
\displaystyle\begin{bmatrix}
\frac{\partial f_1}{\partial x_1} & \cdots & \frac{\partial f_1}{\partial x_N} \\
\vdots & \ddots &\vdots \\
\frac{\partial f_N}{\partial x_1} & \cdots & \frac{\partial f_N}{\partial x_N}
\end{bmatrix}
\end{matrix} \)
\( f\left ( \mathbf{x}+\delta \mathbf{x} \right ) =f\left ( \mathbf{x} \right )+\mathbf{J_f}\left ( \mathbf{x} \right )\delta \mathbf{x} \)
Aquí, \( \mathbf{J_f}\left ( \mathbf{x} \right ) \) es la matriz Jacobiana de dimensión \(N \times N\) de la función \(f\left ( \mathbf{x} \right )\).
\(\mathbf{J_f}\left ( \mathbf{x} \right ) =
\begin{bmatrix}
\frac{\partial f_1}{\partial x_1} & \cdots & \frac{\partial f_1}{\partial x_N} \\
\vdots & \ddots &\vdots \\
\frac{\partial f_N}{\partial x_1} & \cdots & \frac{\partial f_N}{\partial x_N}
\end{bmatrix}\)
El elemento \(i,j\)-ésimo de la matriz es \(\frac{\partial f_i}{\partial x_j}\).
A partir de la expresión
\(f\left ( \mathbf{x}+\delta \mathbf{x} \right ) =f\left ( \mathbf{x} \right )+\mathbf{J}\left ( \mathbf{x} \right )\delta \mathbf{x}\)
obtenemos \(\delta \mathbf{x}\). Aquí \(f\left ( \mathbf{x}+\delta \mathbf{x} \right )=0\).
\(\delta \mathbf{x}=-\mathbf{J}\left ( \mathbf{x} \right )^{-1} f\left ( \mathbf{x} \right )\)
Continuando,
\( \mathbf{x}+\delta \mathbf{x}=\mathbf{x}-\mathbf{J}\left ( \mathbf{x} \right )^{-1} f\left ( \mathbf{x} \right )\)
y
\( \mathbf{x}_{n+1}=\mathbf{x}_n-\mathbf{J}\left ( \mathbf{x}_n \right )^{-1} f\left ( \mathbf{x}_n \right )\)
Aquí, \(\mathbf{x}_n\) es el vector solución en la iteración \(n\), mientras que \(\mathbf{x}_0=\left [ x_1,\cdots ,x_N \right ]^{T}\) es el vector inicial.
Forma general:
\(\begin{bmatrix}x_1\\ \vdots \\ x_N\end{bmatrix}_{n+1}=
\begin{bmatrix}x_1\\\vdots\\x_N\end{bmatrix}_{n}-
\mathbf{J\left ( x_n \right )}^{-1}
\begin{bmatrix}f_1\left ( x_n \right )\\\vdots \\f_N\left ( x_n \right )\end{bmatrix}\)