Métodos Numéricos: Método de Newton-Raphson para sistemas de ecuaciones no lineales

Este método es un método de aproximación de soluciones de sistemas de ecuaciones no lineales, por ejemplo,

y=x^2-4x, \hskip1cm (1)

  (x-2)^2+(y+4)=16.\hskip1cm (2)

Observemos que resolver este sistema de ecuaciones consiste en hallar los puntos (x_0,y_0) de intersección de la parábola asociada a la primera ecuación, la cual puede reescribirse como

y=(x-2)^2-4,\hskip1cm (1.a)

o de manera equivalente como

y=x(x-4)\hskip1cm (1.b)

De $latex (1.a)$ observemos que la parábola abre hacia arriba y que tiene vértice en el punto (2,-4); mientras que de (1.b) se obtiene que intersecta el eje horizontal en los puntos (0,0) y (4,0). Por otra parte, asociada a la ecuación (2) tenemos la circunferencia de radio 4 centrada en el punto (2,-4). El lector podrá graficar sin ninguna dificultad ambos objetos geométricos simultáneamente y observar que deben existir dos puntos de intersección entre la parábola y la circunferencia dadas.

Recordemos que en el caso unidimensional, el método de Newton-Raphson, para aproximar raíces de ecuaciones de la forma

f(x)=0 \hskip1cm (3)

consiste en iterar la función

g(x)=x-\frac{f(x)}{f'(x)}

de modo que al aproximar  puntos fijos de la función g(x) estamos aproximando soluciones de la ecuación f(x)=0; es decir, raíces de la función f(x).

La construcción unidimensional del método de Newton-Raphson se basa en el uso de la aproximación, a primer orden o lineal, de la función f(x), mediante el desarrollo de Taylor alrededor de un punto x_0, aproximación dada por la recta tangente a la gráfica de la función f(x) en el punto (x_0,f(x_0)), y cuya ecuación corresponde a

y=f(x_0)+f'(x_0)(x-x_0) \hskip1cm (4)

Recordemos que la aproximación es para puntos x suficientemente cercanos al punto x_0. Nos es difícil que el lector verifique que la intersección de esta recta tangente (4) con el eje horizontal se da cuando x=x_0-\frac{f(x_0)}{f'(x_0)}. Así, una mejor aproximación de la raíz de f(x) estará dada por el punto

x_1=x_0-\frac{f(x_0)}{f'(x_0)}

siempre que f'(x_0) no se anule o tome un valor muy cercano a cero.

De esta forma, el método de Newton-Raphson, consiste en construir la sucesión de valores \left( x_k\right)_k mediante el algoritmo

x_{k+1}=x_k-\frac{f(x_k)}{f'(x_k)} \hskip1cm (5)

Antes de intentar obtener un algoritmo semejante al anterior para el sistema de ecuaciones (1)-(2), en caso de ser posible, reescribamos dicho sistema de ecuaciones en una forma semejante a (3); es decir,

F(X)=\overline {\bf 0},

donde X=(x,y)^t y \overline 0=(0,0)^t.

Observemos que si tomamos

f_1(x,y)=y-x^2-4x

y

f_2(x,y)=(x-2)^2+(y+4)^2-16,

entonces resolver el sistema de ecuaciones (1)-(2) es equivalente a resolver el sistema

\begin{array}{lcl} f_1(x,y) & = & 0 \\ f_2(x,y) & = & 0 \end{array}

Además, si tomamos

X=(x,y)^t,

 F(X)= \left(\begin{array}{c} f_1(x,y)  \\ f_2(x,y) \end{array}\right)=\left(\begin{array}{c}f_1(x,y)\\f_2(x,y)\end{array}\right),

ahora es claro que el sistema de ecuaciones (1)-(2) es equivalente al sistema

F(X)=\overline{\bf 0}. \hskip2cm (6).

Lo expuesto anteriormente nos muestra que es posible cualquier sistema de dos ecuaciones no lineales en la forma matricial  (6). Más aún, cualquier sistema de n ecuaciones (lineales o no lineales) con m incógnitas puede escribirse en la forma (6).

Ahora sí, está justificada la posibilidad de que exista un algoritmo semejante al dado por (5) para sistemas de ecuaciones no lineales, los cuales tienen la forma (6).

La idea unidimensional consistió en usar la aproximación de Taylo de orden 1 para f cerca de un punto x=x_0.

Para una función g: \mathbb{R}^2\rightarrow  \mathbb{R}, su aproximación lineal para puntos (x,y) cercanos al punto (x_0,y_0) está dada por

g(x,y)\approx g(x_0,y_0)+\frac{\partial g}{\partial x}(x-x_0)+\frac{\partial g}{\partial y}(y-y_0).

Así, para  funciones f_1(x,y) y f_2(x,y) tendremos

f_1(x,y)\approx f_1(x_0,y_0)+\frac{\partial f_1}{\partial x}(x-x_0)+\frac{\partial f_1}{\partial y}(y-y_0),\hskip1cm (7.1)

y

f_2(x,y)\approx f_2(x_0,y_0)+\frac{\partial f_2}{\partial x}(x-x_0)+\frac{\partial f_2}{\partial y}(y-y_0). \hskip1cm (7.2)

Ahora procedamos a escribir estas dos últimas ecuaciones en forma vectorial usando las siguientes igualdades

X=\left(\begin{array}{c} x \\ y \end{array}\right)\hskip0.5cm X_0=\left(\begin{array}{c}x_0\\y_0\end{array}\right),\hskip0.5cm X-X_0=\left(\begin{array}{c} x-x_0\\y-y_0\end{array}\right),

F(X)=\left(\begin{array}{c} f_1(X) \\ f_2(X) \end{array}\right)=\left(\begin{array}{c} f_1(x,y)\\f_2(x,y)\end{array}\right).

DF(X_0)=\left(\begin{array}{cc} \frac{\partial f_1(x_0,y_0)}{\partial x}  & \frac{\partial f_1(x_0,y_0)}{\partial y} \\ \frac{\partial f_2(x_0,y_0)}{\partial x} & \frac{\partial f_2(x_0,y_0)}{\partial y} \end{array}\right)

Ahora, usted, amable lector de este blog, puede comprobar mediante un cálculo directo que el sistema (7.1)-(7.2) puede escribirse matricialmente como

F(X)\approx F(X_0)+DF(X_0)(X-X_0) \hskip1cm (8)

Como deseamos aproximar soluciones de F(X)=\overline {\bf 0}, de esta igualdad y (8), se tiene que el vector X que satisface ambas ecuaciones debe estar dado por

X=X_0-\left( DF(X_0)\right)^{-1}F(X_0)   \hskip1cm (9)

Por lo que definimos a la nueva aproximación de la raíz de F(x) como

X_1=X_0-\left( DF(X_0)\right)^{-1}F(X_0)  \hskip1cm (10),

De esta forma, al algoritmo asociado a este esquema desarrollado está dado por

X_{k+1}=X_k-\left( DF(X_k)\right)^{-1}F(X_k)  \hskip1cm (11),

el cual es la generalización del algoritmo unidimensional (5).

Es importante notar que en el caso bidimensional, el cálculo de las iteraciones puede no presentar demasiadas dificultades, pues al saber que  una matriz

A=\left(\begin{array}{ c c} a_{11} & a_{12}\\a_{21} & a_{22}\end{array}\right)

con determinante no nulo, su matriz inversa está dada por

A^{-1}=\frac{1}{a_{11}a_{22}-a_{21}a_{12}}\left(\begin{array}{ c c} a_{22} & -a_{12}\\-a_{21} & a_{11}\end{array}\right)=\frac{1}{det(A)}\left(\begin{array}{ c c} a_{22} & -a_{12}\\-a_{21} & a_{11}\end{array}\right)

De esta forma, la expresión (11), puede reescribirse al tomar

X_k=\left(\begin{array}{c} x_k\\y_k\end{array}\right)

F(X_k)=\left(\begin{array}{c} f_1(x_k,y_k)\\f_2(x_k,y_k)\end{array}\right)

DF(X_k)^{-1}=\frac{1}{det\left( DF(X_k)\right)}\left(\begin{array}{cc} \frac{\partial f_2(x_k,y_k)}{\partial y_k} & -\frac{\partial f_1(x_k,y_k)}{\partial y_k} \\- \frac{\partial f_2(x_k,y_k)}{\partial x_k} &\frac{\partial f_1(x_k,y_k)}{\partial x_k}\end{array}\right)

donde

det(DF(X_k))=\frac{\partial f_1(x_k,y_k)}{\partial x_k}\frac{\partial f_2(x_k,y_k)}{\partial y_k}-\frac{\partial f_2(x_k,y_k)}{\partial x_k}\frac{\partial f_1(x_k,y_k)}{\partial y_k}

y obtener las siguientes expresiones para el algoritmo de Newton-Raphson en el caso bidimensional, coordenada a coordenada

x_{k+1}=x_k-\frac{1}{det(F(X_k))}\left( f_1(x_k,y_k)\frac{\partial f_2(x_k,y_k)}{\partial y_k}-f_2(x_k,y_k)\frac{\partial f_1(x_k,y_k)}{\partial y_k}\right),\hskip1cm (12.1)

y_{k+1}=y_k-\frac{1}{det(F(X_k))}\left( -f_1(x_k,y_k)\frac{\partial f_2(x_k,y_k)}{\partial x_k}+f_2(x_k,y_k)\frac{\partial f_1(x_k,y_k)}{\partial x_k}\right). \hskip1cm (12.2)

Determinar fórmulas equivalentes a las anteriores para sistemas de más ecuaciones, es claramente, bastante laborioso.

En lugar de ello, no tomaremos el camino de calcular la matriz inversa de

DF(X_k)

pues su cálculo implica dificultades numéricas. Mejor, reescribamos la ecuación (11) de manera más amigable.

X_{k+1}-X_k=-\left(DF(X_k)\right)^{-1}F(X_k)

o equivalentemente,

DF(X_k)\left(X_{k+1}-X_k\right)=-F(X_k)\hskip1cm (13).

Hasta el momento, suponemos que se conoce el elemento X_k de la sucesión de aproximaciones de la raíz de F(X) y que deseamos calcular el elemento X_{k+1} de la sucesión de aproximaciones de la raíz.

Si escribimos

\Delta_k=X_{k+1}-X_k\hskip1cm (14),

el sistema (13) tiene la forma

DF(X_k)\Delta_k=-F(X_k)\hskip1cm (15)

y que representa un sistema de ecuaciones lineales donde conocemos a la matriz DF(X_k) y al vector F(X_k), pues conocemos X_k. Resolviendo el sistema (15) para \Delta_k, podemos obtener X_{k+1}, pues de (14), tenemos que

X_{k+1}=\Delta_k+X_k.\hskip1cm (16).

 

METODOLOGÍA NEWTON-RAPHSON PARA SISTEMAS DE ECUACIONES NO LINEALES

Lo desarrollado anteriormente nos provee de una metodología para calcular aproximaciones de raíces de un sistema de ecuaciones  no lineales
F(x)=\overline{\bf 0}.
I) Considerar una primera aproximación X_0 de la raíz de F(X).
II) Calcular
F(X_0).
III) Calcular
DF(X_0).
IV) Obtener la solución \Delta_0 del sistema de ecuaciones lineales
DF(X_0)\Delta_0=-F(X_0).
V) Obtener la aproximación X1 usando X0 y \Delta_0, mediante la igualdad
X_1=\Delta_0+X_0.
Posteriormente repetir este proceso para obtener X_2 usando X_1 y el procedimiento anterior.
VI) Parar el proceso cuando se satisfaga alguna de las condiciones para una precisión \epsilon>0 deseada.
VI.1)   \|X_{k+1}-X_k\|<\epsilon, o
VI.2)  \frac{\|X_{k+1}-X_k\|}{\|X_{k+1}\|}<\epsilon, o
VI.3) \|F(X_k)\|<\epsilon.

USO DE MATLAB EN NEWTON-RAPHSON

A continuación veremos un ejemplo para aproximar una solución de un sistema de dos ecuaciones no lineales con dos incógnitas. Con el fin de revisar la metodología desarrollada y practicar el uso de dicha metodología, usaremos comandos de Matlab, que esperamos el lector interesado podrá seguir en su computadora para verificar que, efectivamente, ha comprendido el algoritmo que hemos visto.

Consideremos el sistema de ecuaciones lineales

x^2+xy=10 \hskip2cm y+3xy^2=57

Observemos que x=2, y=3 es una solución de este sistema. Observemos que el sistema de ecuaciones anterior se puede escribir como

x^2+xy-10=0 \hskip2cm y+3xy^2-57=0

En el entorno de Matlab usaremos la Caja de Herramienta de cálculo simbólico. Usaremos tres colores para distinguir entre

i) Los comentarios que hagamos para desarrollar el algoritmo de Newton-Raphson, y

ii) Los comandos en el entorno de Matlab

iii) La salida correspondiente al comando implementado

Así pues, manos a la obra.

Primero definamos las variables simbólicas con las que trabajaremos

>> syms x y

 Ahora definamos en una sola línea las funciones 

f_1(x,y)=x^2+xy-10, \hskip 1cm f_2(x,y)=y+3xy^2-57

Para ello, demos el siguiente comando en el entorno y posteriormente apretemos la tecla de ENTER para obtener lo siguiente

>> f1=x^2+x*y-10, f2=y+3x*y^2-57

f1 =

x*y^2 + y – 10

f2 =

3*x*y^2 + y – 57

Ahora calculemos la matriz jacobina de la función F=\left(\begin{array}{c} f_1(x,y)\\ f_2(x,y)\end{array}\right) con los comandos siguiente y obtenemos la salida

>> DFX=[diff(f1,x) diff(f1,y); diff(f2,x) diff(f2,y) ]

DFX =

[ y^2, 2*x*y + 1]
[ 3*y^2, 6*x*y + 1]

El lector puede comprobar mentalmente que esta matriz corresponde a la matriz jacobiana de F.

Ahora definamos la función F como matriz columna, que tiene como coordenadas a las funciones f_1 y f_2

>> F=[f1;f2]
F =
x*y^2 + y – 10
3*x*y^2 + y – 57

Ya tenemos los elementos simbólicos, variables x,y, función F y su jacobiana DF, que son los ingredientes para empezar a iterar nuestro algoritmo dada en la ecuación (15). Además, necesitamos una primera aproximación X_0 de la raíz. Para ello usaremos el punto

X_0=\left(\begin{array}{c} 1.5\\ 3.5\end{array}\right).

A continuación evaluemos F(X_0) y DF(X_0) mediante la sucesión de comandos

>> FX0=subs(F,{x,y},{1.5,3.5})

FX0 =

95/8
13/8  

>> DFX0=subs(DFX,{x,y},{1.5,3.5})

DFX0 =

[ 49/4, 23/2]
[ 147/4, 65/2]

Ya contamos con los elementos de los pasos II y III descritos en nuestra metodología y debemos resolver el sistema de ecuaciones lineales descrito en el apartado III. Para esto, usemos los siguientes comandos, con lo que obtenemos la correspondiente salida \Delta_0, símbolo que en los comandos reemplazaremos por D0.

Para  calcular la siguiente iteración X_1, de cuerdo con el inciso IV descrito en nuestra metodología, definamos nuestra primera aproximación X_0 vectorialmente y calculemos X_1 de acuerdo con IV. Para ello, debemos saber que el sistema de ecuaciones lineales AX=b se resuelve en la ventana de comandos de Matlab mediante el comando

>> X=A\b

Ahora sí, calculemos X_1 como sigue:

>> D0=DFX0\-FX0

D0 =

1469/98
-17

>> X0=[1.5;3.5]

X0 =

1.5000
3.5000

>> X1=D0+X0

X1 =

808/49
-27/2

Si deseamos conocer los valores de las entradas en su presentación decimal, debemos hacer uso del comando vpa , (variable precision arithmetic) 

>> vpa(X1)

ans =

16.489795918367346938775510204082
-13.5

A primera vista esta aproximación no pareciera afortunada, pero si continuamos iterando bajo el proceso anterior en Matlab, obtenemos los siguientes datos

CR1 =subs(F,X,X1′)

2981.7653061224489775271617730823
8945.2959183673469325814853192469

X2 =

0.12894375857338820301783264746228
-13.5

CR2 =subs(F,X,X2′)

-1.1754943508222875079687365372222e-37
-3.5264830524668625239062096116667e-37

X3 =

0.12894375857338820301783264746228
-13.5

CR3 =subs(F,X,X3′)

0
0

X4 =

0.12894375857338820301783264746228
-13.5

CR4 =subs(F,X,X4′)

0
0

X5 =

0.12894375857338820301783264746228
-13.5

CR5 =subs(F,X,X5′)

0
0

X6 =

0.12894375857338820301783264746228
-13.5

CR6 =

0
0

X7 =

0.12894375857338820301783264746228
-13.5

CR7 =

0
0

De lo anterior vemos a partir de la tercera iteración, X3, obtenemos una muy buena aproximación de una solución, pues hemos evaluado la función F en X3 y hemos obtenido el vector \overline{\bf 0}.

En el siguiente documento pdf encontrarán la descripción del método de Newton-Raphson para sistemas de ecuaciones no lineales, así como ejemplo trabajado en Matlab, así como una serie de ejercicios.

Newton-Raphson Sistemas

Anuncios

Métodos numéricos: valores, vectores propios y radio espectral de una matriz

En el siguiente documento pdf encontrarán ejemplos de cálculo de valores y vectores propios de una matriz, así como de ejercicios para trabajar y complementar la comprensión de este tema. También se da la definición de radio espectral de una matriz, ejemplos y algunos ejercicios

Notaspdf

Ecuaciones diferenciales lineales 2do orden no homogéneas: coeficientes indeterminados o ¿podemos adivinar una solución?

 

 

En el siguiente archivo pdf encontrará un ejemplo de como se resuelve una ecuación diferenciales lineal de segundo orden no homogénea con condiciones iniciales de posición y velocidad:

 

y''(x)+5y'(x)+6y(x)=(x+1)\sin(x)  \hskip2cm y(0)=y'(0)=0

solución ejemplo

Debemos leer el archivo con cuidado y haciendo uso de papel y lápiz para llevar a cabo de manera simultánea los cálculos realizados, con el fin de lograr una mejor comprensión de la metodología usada para la solución del ejemplo propuesto.

Recordemos que para hallar la solución general y_g(x) de la ecuación diferencial lineal de segundo orden y no homogénea

y''(x)+p(x)y'(x)+q(x)y(x)=f(x)    (1)

requerimos de dos ingredientes:

i) La solución general y_h(x), de la ecuación lineal de segundo orden homogénea

y''(x)+p(x)y'(x)+q(x)y(x)=0      (2)

Para calcular dicha solución requerimos de dos soluciones y_1(x)y_2(x) de (2) que sean linealmente independientes; es decir, tal que su wronskiano

W[y_1,y_2](x)=y_1(x)y_2'(x)-y_1'(x)y_2(x)

 

no se anule.

 

Así, la solución general y_h(x) de la ecuación homogénea (2) se puede escribir como combinación lineal de y_1(x) y de y_2(x); es decir, existen constantes c_1 y c_2 tales que

y_h(x)=c_1y_1(x)+c_2y_2(x).

 

ii) Una solución particular y_p(x) de la ecuación no homogénea (1).

 

En el ejemplo a que se hace mención en esta entrada se hace uso del método de coeficientes indeterminados y que algún texto da en llamar de adivinación razonada o juiciosa.

 

Una vez que contamos con la solución genera

y_h(x)=c_1y_1(x)+c_2y_2(x)

de la ecuación homogénea y con una solución particular y_p(x) de la ecuación no homogénea (1),  tenemos que la solución general y_g(x) de la ecuación no homogénea (1) tiene la forma

 

y_g(x)=y_h(x)+y_p(x)=c_1y_1(x)+c_2y_2(x)+y_p(x),

donde las constantes c_1 y c_2 dependen de las condiciones iniciales y(x_o)=u_0 y y'(x_0)=v_0.