Informática de Sistemas
Cálculo Numérico con Mathematica

Tema 2
Resolución numérica de ecuaciones y sistemas no lineales

Nota: Se escribe al principio de la práctica el comando que permite borrar todos los símbolos almacenados (en concreto, valores de variables o funciones) desde que se inició la sesión con Mathematica. Así se facilita su uso siempre que sea necesario.

Clear["Global`*"]

Órdenes directas para aproximar soluciones de ecuaciones

Mathematica contempla dos posibilidades de resolver ecuaciones o sistemas de ecuaciones: de forma exacta (mediante cálculo simbólico) o aproximada (mediante cálculo numérico). Debido a que las ecuaciones y sistemas que pueden resolverse de manera exacta son excepcionales, es importante el estudio de su resolución numérica.

Solve

Solve es el comando de Mathematica para la resolución exacta de ecuaciones o sistemas de ecuaciones. Su sintaxis es Solve[ecuación/es,variable/s]. Está especialmente indicado para la resolución exacta de ecuaciones polinómicas de grado menor que cinco, pero también puede utilizarse en otros tipos de ecuaciones: en general, aquéllas que puedan resolverse mediante operaciones algebraicas. Cuando Mathematica no es capaz de proporcionar todas las soluciones de una ecuación mediante el comando Solve, avisa mediante un mensaje como el siguiente: ``the equations appear to involve the variables to be solved for in an essentially non-algebraic way''.
Ilustramos todo lo dicho con los siguientes ejemplos:

Ejemplo de ecuación de segundo grado respecto de la incógnita z:

Solve[3 (z^2) - a z + (a^3) == 0, z]

Una ecuación cúbica, pero sin parámetros:

Solve[4 (x^3) - 10 (x^2) - 6 x + 3 == 0, x]

N[%]

Este resultado lleva a sospechar que todas las raíces son reales. Resolvemos esta duda utilizando el comando Plot:

Plot[4 (x^3) - 10 (x^2) - 6 x + 3, {x, -1, 3.3}]

Una manera de eliminar los números muy pequeños es utilizando el comando Chop, que sustituye todos los números reales menores que (10^(-10)) por el número 0.

Chop[N[Solve[4 (x^3) - 10 (x^2) - 6 x + 3 == 0, x]]]

Una ecuación de grado 4, el grado más alto en que toda ecuación puede resolverse de manera exacta por radicales:

Solve[(x^4) - 2 (x^3) + 3 (x^2) - 4 x + 5 == 0, x]

N[%]

Como recordará y se observa en el ejemplo anterior, las raíces imaginarias de polinomios con coeficientes reales siempre se tienen a pares (una y su conjugada).

Las ecuaciones de grado 5 o superior no son resolubles con exactitud por lo general, aunque Matemática pueda dar un nombre a las raíces, como ocurre a continuación. Lo que sí es posible es aproximar esas raíces con Mathematica con el comando N. Vemos ambas cosas en el próximo ejemplo:

Solve[(x^5) + 3 (x^2) + 7 == 0, x]

N[%]

Lo que no quita que cuando algunas raíces puedan obtenerse de modo exacto Mathematica las obtenga (y no proporcione explícitamente las otras):

Solve[-70 - 11 x - 67 (x^2) - 11 (x^3) + 3 (x^4) - 7 (x^6) + 
          (x^7) - 7 (x^8) + (x^9) == 0, x]

N[%]

Observe cómo se han aproximado las soluciones ±i de la ecuación anterior.

Además de para ecuaciones polinómicas, Solve puede aplicarse a otras ecuaciones resolubles por métodos algebraicos. Por ejemplo:

Solve[(x^2) + (Sqrt[x  + 2]) - 6 == 0, x]

N[%]

Solve puede aplicarse también a sistemas, como se ve en los siguientes casos:

Solve[{(x^2) + y + z - 3 == 0, 3 (z^2) - x - 2  z - y == 0, 
        x - 2  z - y - 1 == 0}, {x, y, z}]

N[%]

Solve[{(x^2) + y + z - 3 == 0, 
        3 (z^2) - x - 2  z - y == 0}, {x, y}]

Como se ha dicho, para la mayor parte de ecuaciones y sistemas Solve no puede hacer nada. Por ejemplo:

Solve[x - Sin[x] == 0, x]

Observe que ni siquiera ha encontrado la solución x = 0, que es trivial.

NSolve

El comando NSolve es útil para encontrar aproximaciones de todas las soluciones de una ecuación polinómica (y de algunas otras ecuaciones). Equivale a aplicar primero el comando Solve y luego el comando N. Como ocurre con éste último comando, se puede aumentar la precisión. La sintaxis de este comando es pues NSolve[ecuación,variable] o NSolve[ecuación,variable,n]. A continuación lo aplicamos a varias ecuaciones que ya se resolvieron mediante el comando Solve.

NSolve[x - Sin[x] == 0, x]

NSolve[4 (x^3) - 10 (x^2) - 6 x + 3 == 0, x]

NSolve[4 (x^3) - 10 (x^2) - 6 x + 3 == 0, x, 22]

NSolve[-70 - 11 x - 67 (x^2) - 11 (x^3) + 3 (x^4) - 7 (x^6) + 
          (x^7) - 7 (x^8) + (x^9) == 0, x]

NSolve[(x^2) + (Sqrt[x  + 2]) - 6 == 0, x]

NSolve[(x^2) + (Sqrt[x  + 2]) - 6 == 0, x, 10]

NSolve puede aplicarse también a sistemas, como se ve en el siguiente caso:

NSolve[{(x^2) + y + z - 3 == 0, 3 (z^2) - x - 2  z - y == 0, 
        x - 2  z - y - 1 == 0}, {x, y, z}]

NSolve[{(x^2) + y + z - 3 == 0, 3 (z^2) - x - 2  z - y == 0, 
        x - 2  z - y - 1 == 0}, {x, y, z}, 20]

FindRoot

El comando de resolución numérica propiamente dicho es FindRoot. Como utiliza métodos numéricos, FindRoot encuentra las soluciones una a una. Es aplicable a todo tipo de ecuaciones en general. Su sintaxis más simple es FindRoot[ecuación,{variable,valorinicial}] que aplica un proceso iterativo de aproximación partiendo de valorinicial como primera aproximación de una solución. Si se desea que el ordenador trabaje con mayor precisión, podemos escribir FindRoot[ecuación,{variable,valorinicial},WorkingPrecision->n]. Apliquemos FindRoot a algunas ecuaciones:

FindRoot[4 (x^3) - 10 (x^2) - 6 x + 3 == 0, {x, 1}]
FindRoot[4 (x^3) - 10 (x^2) - 6 x + 3 == 0, {x, 2}]
FindRoot[4 (x^3) - 10 (x^2) - 6 x + 3 == 0, {x, -1}]

FindRoot puede ser convergente incluso si partimos de un punto muy lejano a una solución (pero no siempre sucederá así):

FindRoot[4 (x^3) - 10 (x^2) - 6 x + 3 == 0, {x, -1000000}]

En la siguiente ejemplo, aun probando con varios números reales de partida, FindRoot no proporciona nada porque las cuatro raíces de la ecuación son complejas. Sin embargo, si partimos de un valor complejo, sí que llega a una solución:

FindRoot[(x^4) - 2 (x^3) + 3 (x^2) - 4 x + 5 == 0, {x, 3}]

FindRoot[(x^4) - 2 (x^3) + 3 (x^2) - 4 x + 5 == 
        0, {x, I}]
FindRoot[(x^4) - 2 (x^3) + 3 (x^2) - 4 x + 5 == 
        0, {x, -I}]
FindRoot[(x^4) - 2 (x^3) + 3 (x^2) - 4 x + 5 == 0, {x, 
        1 + I}]

Otro ejemplo más:

FindRoot[3 + 3  x - 7  (x^2) - (x^3) + 2  (x^4) + 3  (x^7) - 3  (x^8) - (x^9) + 
          2  (x^10) == 0, {x, 3}]

El problema esta vez no es que no haya soluciones reales, sino que hemos partido desde un punto donde FindRoot no converge a alguna de las dos soluciones reales que vimos que tenía esta ecuación. Se puede actuar de varias formas. Una es observar previamente con un Plot dónde se sitúan las raíces, y partir desde puntos relativamente cercanos a ellas. Es lo que hacemos a continuación.

Plot[
      3 + 3 x - 7 (x^2) - (x^3) + 2 (x^4) + 3 (x^7) - 3 (x^8) - (x^9) + 
        2 (x^10), {x, -3, 3}]

Plot[
      3 + 3 x - 7 (x^2) - (x^3) + 2 (x^4) + 3 (x^7) - 3 (x^8) - (x^9) + 
        2 (x^10), {x, -1.6, 1.5}]

FindRoot[3 + 3  x - 7  (x^2) - (x^3) + 2  (x^4) + 3  (x^7) - 3  (x^8) - (x^9) + 
          2  (x^10) == 0, {x, -1.5}]
FindRoot[3 + 3  x - 7  (x^2) - (x^3) + 2  (x^4) + 3  (x^7) - 3  (x^8) - (x^9) + 
          2  (x^10) == 0, {x, -0.5}]

En el siguiente ejemplo nos encontramos con un nuevo problema:

Plot[x - Sin[x], {x, -2 π, 2 π}]

FindRoot[x - Sin[x] == 0, {x, 1}]

Aquí FindRoot no converge muy bien a la única solución x = 0: más adelante veremos que el problema está en que x = 0 es una raíz múltiple (triple en este caso) de la ecuación x  - sen x = 0. Pero podemos conseguir más exactitud con la instrucción WorkingPrecision, opción aplicable a FindRoot y a otras operaciones numéricas que especifica cuántos dígitos de precisión deben mantenerse en los cálculos (la precisión normal está en torno a los 16 dígitos).

FindRoot[x - Sin[x] == 0,{x, 1},WorkingPrecision->20]
FindRoot[x - Sin[x] == 0,{x, 1},WorkingPrecision->30]

FindRoot puede aplicarse también a sistemas, como se en el siguiente caso:

FindRoot[{(x^2) + y + z - 3 == 0, 
        3 (z^2) - x - 2 z - y == 0, 
        x - 2 z - y - 1 == 0}, {{x, 0}, {y, 0}, {z, 
          0}}]
FindRoot[{(x^2) + y + z - 3 == 0, 
        3 (z^2) - x - 2 z - y == 0, 
        x - 2 z - y - 1 == 0}, {{x, 1}, {y, 1}, {z, 
          1}}]
FindRoot[{(x^2) + y + z - 3 == 0, 
        3 (z^2) - x - 2 z - y == 0, 
        x - 2 z - y - 1 == 0}, {{x, -2}, {y, 
          3  I}, {z, -I}}]
FindRoot[{(x^2) + y + z - 3 == 0, 
        3 (z^2) - x - 2 z - y == 0, 
        x - 2 z - y - 1 == 
          0}, {{x, -2}, {y, -3 I}, {z, I}}]

Una cuestión previa sobre los algoritmos de aproximación de soluciones de ecuaciones

En todos los programas que siguen para aproximar la solución de una ecuación f(x) = 0 (o de un sistema de ecuaciones) parecería lógico incluir un criterio de parada en el algoritmo que dijera más o menos lo siguiente:

If[f[(Subscript[x, i])] == 0, Print["La solucion exacta es ", (Subscript[x, i])]; 
      Break[]]

Pero los redondeos pueden hacer que f((Subscript[x, i])) = 0 y que (Subscript[x, i]) no sea en realidad la solución. Veamos un ejemplo: consideramos la función

f[x_] := 
      100 - 75 x - 15 (x^2) + 30 (x^3) - 6 (x^4) - 3 (x^5) + (x^6)

Las soluciones de la ecuación f(x) = 0 son:

Solve[f[x] == 0, x]

Nos fijamos en la primera de sus soluciones:

InputForm[N[-(Sqrt[5])]]

Por ser solución la función se anula en este valor:

f[-2.23606797749979]

Pero, curiosamente, también se anula en valores próximos que no son solución, como por ejemplo el siguiente:

f[-2.2360679780560986]

Lo más sorprendente es que la función no se anule en valores situados entre los dos "ceros" anteriores, como por ejemplo:

f[-2.236067978]

Por ese motivo no hemos incluido en los algoritmos el criterio de parada mencionado, pues a veces es fallido (más adelante, en el penúltimo apartado de esta práctica, lo mostraremos con un ejemplo).

Método de bisección

Programa que implementa el método de bisección

Lo mostramos con un ejemplo. En este primer ejemplo la solución de la ecuación es conocida, lo que permite comparar la aproximación del método con la solución exacta.

La raíz exacta de la ecuación (x^2) = 2 en el intervalo [0, 2] es (Sqrt[2]). Una implementación del método de bisección para calcularla con un error menor que una milésima es la siguiente:

Se sabe que el número de iteraciones suficiente para asegurar el error máximo deseado es el siguiente:

a = 0. ; b = 2. ; cota = (10^(-3)); 
    maxiter = Log[((b - a)/cota)]/Log[2]

Luego se necesitarán 11 iteraciones. De todos modos, el método de bisección se puede programar sin calcular este número (como ocurrirá con los restantes métodos), si se especifica un alto número máximo de iteraciones en la instruccion For:

\(f[x_] := x\^2 - 2; \ a = 0. ; b = 2. ; cota = 10\^\(-3\); 
    For[i = 1, i ≤ 1000, \(i++\), c = \(a + b\)\/2; 
      errmax = \(b - a\)\/2; 
      Print["i=", i, \*""  (Subscript[a, i]) = "", 
        a, \*""  (Subscript[b, i]) = "", b, \*""  (Subscript[x, i]) = "", 
        c, \*""  f((Subscript[x, i])) = "", N[f[c]], "  Cota error = ", 
        errmax]; 
      If[c - a ≤ cota, Print["La solucion aproximada es ", c]; 
        Print["La cota del error que se comete es ", errmax]; 
        Print["El valor que toma la función en la solución aproximada es ", f[c]]; 
        Print["El error que verdaderamente se comete es ", c - \@2]; 
        Break[]]; If[f[c]\ f[a] < 0, b = c, a = c]]\)

Ejemplo

Aproxime la solución de la ecuación (x^3) = 3  x - 1 que se encuentra en el intervalo (0, 1) con un error inferior a una milésima (utilice el método de bisección).

Solución

En primer lugar expresamos la solución tal y como aparece en los apuntes.

\(f[x_] := x\^3 - 3\ x + 1; a = 0; b = 1; cota = 10\^\(-3\); 
    For[i = 1, i ≤ 1000, \(i++\), c = \(a + b\)\/2; 
      errmax = \(b - a\)\/2; 
      Print["i=", i, \*""  (Subscript[a, i]) = "", 
        a, \*""  (Subscript[b, i]) = "", b, \*""  (Subscript[x, i]) = "", 
        c, \*""  f((Subscript[x, i])) = "", N[f[c]], "  Cota error = ", 
        errmax]; 
      If[c - a ≤ cota, Print["La solucion aproximada es ", c]; 
        Print["La cota del error que se comete es ", errmax]; 
        Print["El valor que toma la función en la solución aproximada es ", f[c]]; Break[]]; 
      If[f[c]\ f[a] < 0, b = c, a = c]]\)

Ahora expresamos la solución con valores decimales.

\(f[x_] := x\^3 - 3\ x + 1; a = 0. ; b = 1. ; cota = 10\^\(-3\); 
    For[i = 1, i ≤ 1000, \(i++\), c = \(a + b\)\/2; 
      errmax = \(b - a\)\/2; 
      Print["i=", i, \*""  (Subscript[a, i]) = "", 
        a, \*""  (Subscript[b, i]) = "", b, \*""  (Subscript[x, i]) = "", 
        c, \*""  f((Subscript[x, i])) = "", N[f[c]], "  Cota error = ", 
        errmax]; 
      If[c - a ≤ cota, Print["La solucion aproximada es ", c]; 
        Print["La cota del error que se comete es ", errmax]; 
        Print["El valor que toma la función en la solución aproximada es ", f[c]]; Break[]]; 
      If[f[c]\ f[a] < 0, b = c, a = c]]\)

Comentario

Esta vez no conocemos la solución exacta de la ecuación, como ocurría con el ejemplo anterior. Pero, al ser una ecuación polinómica, podemos obtenerla con los dígitos exactos que deseemos utilizando la instrucción NSolve:

NSolve[f[x]==0,x,18]

Luego la diferencia entre la última aproximación c y la solución exacta es del orden de

c - 0.34729635533386069770343325353862909904`18.

Ejercicio

Dibuje la función f(x)=x sen x. ¿Qué tiene de especial el punto x=0 si lo comparamos con los demás ceros de f? ¿Podría utilizarse el método de bisección para encontrar una raíz con esas características?

Solución

Plot[x Sin[x], {x, -4 π, 4 π}]

El punto x=0 es un cero de la función donde ésta no cambia de signo, con lo es imposible aplicar el método de bisección. Si esto no se percibe con claridad, intente aplicar bisección (por ejemplo, al intervalo [-2, 5]).

Resuelva los ejercicios siguientes

Ejercicio

Utilice el método de bisección para encontrar una aproximación del punto x donde se cortan las gráficas de las funciones f(x) = (Sqrt[(x^2) + 1]) y g(x) = tg x, definidas en el intervalo (0, π/2), de manera que x se aproxime con error menor que (10^(-4)).

Solución

Plot[{(Sqrt[(x^2) + 1]), Tan[x]}, {x, 0, (π/2)}]

f[x_] := (Sqrt[(x^2) + 1]) - Tan[x];

Ejercicio

Utilice el método de bisección para aproximar una solución de la ecuación (Sqrt[x - 3]) = 3  ln(x - 10) con un error menor que (10^(-5)). Programe el algoritmo de modo que: (a) se visualicen las sucesivas aproximaciones; (b) la aproximación final se represente con al menos 8 dígitos de mantisa; (c) se muestre en cada paso una cota del error de aproximación, de manera que en el último paso ésta sea por primera vez menor que (10^(-5)).

Solución

f[x_] := (Sqrt[x + 3]) - 3  Log[x - 10]

Métodos regula falsi y de la secante

Ejemplo

Aproxime la solución de la ecuación (x^3) = 3  x - 1 que se encuentra en el intervalo (0, 1) de modo que las dos últimas aproximaciones difieran en menos de una millonésima (utilice el método de la regula falsi).

Solución

f[x_] := (x^3) - 3 x + 1; 
a = 0; b = 1. ; 
    difmax = (10^(-6)); 
    c  = a - ((f[a] (b - a))/(f[b] - f[a])); 
\(Print[\ "i=1", \*""  \!\(a\_1\) = "", 
      a, \*""  \!\(b\_1\) = "", b, \*""  \!\(x\_1\) = "", 
      c, \*""  f(\!\(x\_1\)) = "", f[c]]; \)
For[i = 2,i ≤ 100,i++,If[f[c] f[a] < 0, b = c, a = c];dif = Abs[a - ((f[a] (b - a))/(f[b] - f[a])) - c];c  = a - ((f[a] (b - a))/(f[b] - f[a]));Print[ "i=",i,"  (Subscript[a, i]) = ",a,"  (Subscript[b, i]) = ",b,"  (Subscript[x, i]) = ",InputForm[c],"  f((Subscript[x, i])) = ",f[c],"  |(Subscript[x, i])-(Subscript[x, i - 1])| = ",dif];If[dif < difmax, 
            Print["La solucion aproximada es ", InputForm[c]]; 
            Print["La diferencia entre las dos últimas \
aproximaciones es ", dif]; 
            Print["El valor que toma la función en la soluci\
ón aproximada es ", f[c]]; Break[]]]

Comentario

Comparamos la aproximación con la que obtuvimos antes mediante el comando FindRoot:

c - 0.34729635533386069770343325353862909904`18.

Ejemplo

Aproxime la solución de la ecuación (x^3) = 3  x - 1 que se encuentra en el intervalo (0, 1) de modo que las dos últimas aproximaciones difieran en menos de una millonésima (utilice el método de la secante).

Solución

f[x_] := (x^3) - 3  x + 1; a = 0. ; b = 1. ; 
    difmax = (10^(-6));  
Print["a = ",  a, "  f(a) = ",   f[a], "  b = ",  
      b, "  f(b) = ",  f[b]]
\(For[i = 1, i ≤ 100, \(i++\), 
      c = b - \(f[b] \((b - a)\)\)\/\(f[b] - f[a]\); 
      Print[\ "i=", i, \*""  (Subscript[x, i]) = "", 
        InputForm[c], \ \*""  f((Subscript[x, i])) = "", 
        f[c], \*""  |(Subscript[x, i])-(Subscript[x, i - 1])| = "", Abs[c - b]]; 
      If[Abs[c - b] < difmax, 
        Print[\*""La solución aproximada es (Subscript[x, i]) = "", \ 
          InputForm[c]]; 
        Print["La diferencia entre las dos últimas aproximaciones es "\ \ , Abs[c - b]]; 
        Print["El valor que toma la función en la solución aproximada es "\ \ , f[c]]; Break[]]; \ a = b; b = c]\)

Comentario

Comparamos la aproximación con la que obtuvimos antes mediante el comando FindRoot:

c - 0.34729635533386069770343325353862909904`18.

Resuelva los ejercicios que hay después del ejemplo:

Ejemplo

Considere la ecuación f(x) = (x^3) - x  - 3 = 0. Implemente el método de la secante con la condición de parada |f(x)| < (10^(-6)) para obtener una aproximación de la única solución real de dicha ecuación. ¿Cuál es la solución aproximada que se obtiene si si el intervalo de partida es [-10, 20]? ¿Y si es [1, 2]? ¿Cuál es la solución obtenida al utilizar los tres comandos directos de Mathematica? Comente los resultados.

Solución

f[x_] := (x^3) - x - 3;a = -10. ;b = 20. ;fmax = (10^(-6));
Print["a = ",  
      a, "  f(a) = ",   f[a], "  b = ",  b, "  f(b) = ",  
      f[b]]
For[i = 1,i ≤ 20,i++,c = b - ((f[b] (b - a))/(f[b] - f[a]));\(Print[\ "i=", i, \*""  \!\(x\_i\) = "", 
            InputForm[c], \ \*""  f(\!\(x\_i\)) = "", 
            f[c] \*"",  |\!\(x\_i\)-\!\(x\_\(i - 1\)\)| = "", 
            Abs[c - b]]\);If[Abs[f[c]] < fmax,\(Print[\*""La solución aproximada es (Subscript[x, i]) = "", \ InputForm[c]]; 
              Print["El valor que toma la función en la solución aproximada es "\ \ , f[c]]; 
              Print["La diferencia entre las dos últimas aproximaciones es "\ \ , Abs[c - b]]; Break[]\)]; a = b;b = c]

En el siguiente caso, al partir de un intervalo de menor amplitud, lo normal será que se necesiten menos iteraciones:

f[x_] := (x^3) - x - 3;a = 1. ;b = 2. ;fmax = (10^(-6));
Print["a = ",  
      a, "  f(a) = ",   f[a], "  b = ",  b, "  f(b) = ",  
      f[b]]
For[i = 1,i ≤ 20,i++,c = b - ((f[b] (b - a))/(f[b] - f[a]));\(Print[\ "i=", i, \*""  \!\(x\_i\) = "", 
            InputForm[c], \ \*""  f(\!\(x\_i\)) = "", 
            f[c] \*"",  |\!\(x\_i\)-\!\(x\_\(i - 1\)\)| = "", 
            Abs[c - b]]\);If[Abs[f[c]] < fmax,\(Print[\*""La solución aproximada es (Subscript[x, i]) = "", \ InputForm[c]]; 
              Print["El valor que toma la función en la solución aproximada es "\ \ , f[c]]; 
              Print["La diferencia entre las dos últimas aproximaciones es "\ \ , Abs[c - b]]; Break[]\)]; a = b;b = c]

Como es una ecuación polinómica, podemos utilizar los comandos Solve y NSolve (este último, con mayor o menor aproximación):

Solve[f[x] == 0, x]
NSolve[f[x] == 0, x]
NSolve[f[x] == 0, x, 20]

Y también es posible usar la instrucción FindRoot, con mayor o menor aproximación:

FindRoot[f[x] == 0, {x, 2}]
FindRoot[f[x] == 0, {x, 2}, 
      WorkingPrecision -> 20]

Las diferencias entre esta última aproximación y las dos que obtuvimos por el método de la secante son:

1.6716998816571609697481661411876159019893736795546668993969`20.-1.6716998876530076
1.6716998816571609697481661411876159019893736795546668993969`20.-1.6716998816832451

Ejercicio

Considere la ecuación x = cos x + 3. Utilice el método de la secante con la condición de parada |f(x)| < (10^(-6)) para encontrar la única solución de esa ecuación.

Solución

f[x_] := x - Cos[x] - 3

Ejercicio

Considere la ecuación (x^3) = x cos x + 20 (Sqrt[ (x^4) + (x^2)]) y utilice el método de la secante con la condición de parada |(Subscript[x, n + 1]) - (Subscript[x, n])|< (10^(-6)) para obtener la única solución no nula de dicha ecuación.

Solución

f[x_] := (x^3) - x Cos[x] - 20 (Sqrt[ (x^4) + (x^2)])

Método de Newton-Raphson

Ejemplo

Aproxime la solución de la ecuación (x^3) = 3  x - 1 que se encuentra en el intervalo (0, 1) de modo que las dos últimas aproximaciones difieran en menos de una millonésima (utilice el método de Newton-Raphson).

Solución

f[x_] := (x^3) - 3  x + 1; a =  .5; 
    difmax = (10^(-6)); 
\(Print[\ "i=", 1, \*""  \!\(x\_1\) = "", 
      a, \ \*""  f(\!\(x\_1\)) = "", f[a]]; \)
\(For[i = 2, i ≤ 100, \(i++\), b = a - f[a]/\(f'\)[a]; 
      Print[\ "i=", i, \*""  (Subscript[x, i]) = "", 
        InputForm[b], \ \*""  f((Subscript[x, i])) = "", 
        f[b], \ \*""  |(Subscript[x, i])-(Subscript[x, i - 1])| = "", 
        Abs[b - a]]; 
If[Abs[b - a] < difmax, 
        Print[\*""La solución aproximada es (Subscript[x, i]) = "", \ 
          InputForm[b]]; 
        Print["La diferencia entre las dos últimas aproximaciones es "\ \ , Abs[b - a]]; 
        Print["El valor que toma la función en la solución aproximada es "\ \ , f[b]]; Break[]]; \ a = b]\)

Comentario

Comparamos la aproximación con la que obtuvimos antes mediante el comando FindRoot:

b - 0.34729635533386069770343325353862909904`18.

Sale 0 porque, con la precisión utilizada para b, todos los dígitos de esta aproximación son correctos.

Ejemplo

Utilizando el método de Newton-Raphson y planteando la función adecuada se puede aproximar la raíz cúbica real de cualquier número real. En particular, obtenga un valor aproximado s de la raíz cúbica de 26 asegurando que |(s^3) - 
          26|<(10^(-4)).

Solución

Como s = (26^(1/3)) es el único cero real de la función f(x) = (x^3) - 26, se trata de aproximar s utilizando como condición de parada |f(x)| < (10^(-4)):

f[x_] := (x^3) - 26; a = 3. ; 
    fmax = (10^(-4)); 
\(Print[\ "i=", 1, \*""  \!\(x\_1\) = "", 
      a, \ \*""  |f(\!\(x\_1\))| =|\!\(x\_1\^3\)-26|= "", 
      f[a]]; \)
\(For[i = 2, i ≤ 10, \(i++\), b = a - f[a]/\(f'\)[a]; 
      Print[\ "i=", i, \*""  (Subscript[x, i]) = "", 
        InputForm[b], \*""  |(Subscript[x, i])-(Subscript[x, i - 1])| = "", 
        Abs[b - a], \*""  |f((Subscript[x, i]))| =|((Subscript[x, i])^3)-26| = "", 
        Abs[f[b]]]; 
If[Abs[f[b]] < fmax, 
        Print[\*""La aproximación pedida de (26^(1/3)) es (Subscript[x, i]) = "", \ InputForm[b]]; 
        Print["La diferencia entre el cubo de la última aproximación y 26 es "\ \ , Abs[f[b]]]; \ 
        Print["La diferencia entre las dos últimas aproximaciones es "\ \ , Abs[b - a]]; Break[]]; a = b]\)

Comentario

La diferencia entre la aproximación y el valor exacto de s es:

b - (26^(1/3))

Resuelva los ejercicios siguientes:

Ejercicio

La función (x^3) + (e^x) = 10 tiene una única raíz real. Determine un intervalo en el cuál esté la raíz y utilice los métodos de bisección, regula falsi, secante y Newton-Raphson hasta que la diferencia entre dos aproximaciones consecutivas sea menor que (10^(-4)).

Solución

f[x_] := (x^3) + (E^x) - 10

Ejercicio

Utilice los métodos de bisección, regula falsi, secante y Newton-Raphson para obtener aproximadamente las soluciones de la ecuación (x^3) - 10 = (e^x), de manera que la diferencia entre dos aproximaciones consecutivas sea menor que (10^(-4)).

Solución

f[x_] := (x^3) - (E^x) - 10

Ejercicio

Utilizando el método de Newton-Raphson y planteando la función adecuada se puede aproximar la raíz cuadrada positiva de cualquier número real positivo. En particular, obtega un valor aproximado s de (Sqrt[26]) asegurando que |(s^2) - 
          26|<(10^(-4)).

Solución

Ejercicio

Considere la función f(x) = x (e^(-x)) + 1. Utilice el método de Newton-Raphson para aproximar el único cero de dicha función, con la condición de parada |(Subscript[x, n + 1])-(Subscript[x, n])|<(10^(-4)). Compare este resultado con el obtenido utilizando comandos directos de Mathematica. ¿Puede explicar lo que ocurre si se toma como aproximación inicial (Subscript[x, 1]) = 4?

Solución

f[x_] := x (E^(-x)) + 1

Método de Newton-Raphson modificado para la obtención de raíces múltiples

Ejemplo

Se sabe que el siguiente polinomio tiene sus raíces reales en el intervalo [-101, 101]:

p[x_] := 
      100 - 75 x - 15 (x^2) + 30 (x^3) - 6 (x^4) - 3 (x^5) + 
        (x^6);

Localizamos las soluciones de la ecuación p(x) = 0 representando la función p(x):

Plot[p[x], {x, -101, 101}]

Después de hacer algunas representaciones más, podemos encontrar que hay dos soluciones reales, \(TraditionalForm\`s\_1 ∈ \((\(-3\), \(-2\))\)\ \ y\ \ s\_2 ∈ \((2, 3)\)\), que aparentan ser múltiples:

Plot[p[x], {x, -3, 3}]

Aproximación de la raíz positiva mediante los métodos de Newton Raphson

Aproximamos la raíz positiva mediante el método de Newton-Raphson:

a = 2. ; difmax = (10^(-4));  
\(Print[\ "i=", 1, \*""  \!\(x\_1\) = "", 
      a, \ \*""  p(\!\(x\_1\)) = "", p[a]]; \)
\(For[i = 2, i ≤ 20, \(i++\), b = a - p[a]/\(p'\)[a]; 
      Print[\ "i=", i, \*""  (Subscript[x, i]) = "", 
        InputForm[b], \ \*""  p((Subscript[x, i])) = "", 
        p[b], \ \*""  |(Subscript[x, i])-(Subscript[x, i - 1])| = "", 
        Abs[b - a]]; 
      If[Abs[b - a] < difmax, 
        Print[\*""La solución aproximada es (Subscript[x, i]) = "", \ 
          InputForm[b]]; 
        Print["La diferencia entre las dos últimas aproximaciones es "\ \ , Abs[b - a]]; 
        Print["El valor que toma la función en la solución aproximada es "\ \ , p[b]]; Break[]]; \ a = b]\)

Se observa que la convergencia es muy lenta, por lo que probamos con el método modificado para raíces dobles:

a = 2. ; 
    difmax = (10^(-4)); 
\(Print[\ "i=", 
      1, \*""  \!\(x\_1\) = "", a, \ \*""  p(\!\(x\_1\)) = "", 
      p[a]]; \)
For[i = 2,i ≤ 20,i++,b=a-2p[a]/p'[a];\(Print[\ "i=", i, \*""  \!\(x\_i\) = "", 
            InputForm[b], \ \*""  p(\!\(x\_i\)) = "", 
            p[b], \ \*""  |\!\(x\_i\)-\!\(x\_\(i - 1\)\)| = "", 
            Abs[b - a]]\);\(If[Abs[b - a] < difmax, 
            Print[\*""La solución aproximada es (Subscript[x, i]) = "", \ InputForm[b]]; 
            Print["La diferencia entre las dos últimas aproximaciones es "\ \ , Abs[b - a]]; 
            Print["El valor que toma la función en la solución aproximada es "\ \ , p[b]]; Break[]]\); a = b]

Probamos también con el método modificado para raíces triples:

a = 2. ; 
    difmax = (10^(-4)); 
\(Print[\ "i=", 
      1, \*""  \!\(x\_1\) = "", a, \ \*""  p(\!\(x\_1\)) = "", 
      p[a]]; \)
For[i = 2,i ≤ 20,i++,b=a-3p[a]/p'[a];\(Print[\ "i=", i, \*""  \!\(x\_i\) = "", 
            InputForm[b], \ \*""  p(\!\(x\_i\)) = "", 
            p[b], \ \*""  |\!\(x\_i\)-\!\(x\_\(i - 1\)\)| = "", 
            Abs[b - a]]\);\(If[Abs[b - a] < difmax, 
            Print[\*""La solución aproximada es (Subscript[x, i]) = "", \ InputForm[b]]; 
            Print["La diferencia entre las dos últimas aproximaciones es "\ \ , Abs[b - a]]; 
            Print["El valor que toma la función en la solución aproximada es "\ \ , p[b]]; Break[]]\); a = b]

Como la convergencia es más lenta que en el caso anterior, podemos sospechar que la raíz es doble.

Puede probarse que esa raíz es exactamente igual a (Sqrt[5]) y que es, efectivamente, doble. Comparamos la solución aproximada obtenida en el segundo caso con la verdadera solución:

InputForm[(Sqrt[5])-2.2360679618974384]

Aproximación de la raíz negativa mediante los métodos de Newton Raphson

Mediante el método de Newton-Raphson normal:

a = -2. ; difmax = (10^(-4)); 
\(Print[\ "i=", 1, \*""  \!\(x\_1\) = "", 
      a, \ \*""  p(\!\(x\_1\)) = "", p[a]]; \)
\(For[i = 2, i ≤ 20, \(i++\), b = a - p[a]/\(p'\)[a]; 
      Print[\ "i=", i, \*""  (Subscript[x, i]) = "", 
        InputForm[b], \ \*""  p((Subscript[x, i])) = "", 
        p[b], \ \*""  |(Subscript[x, i])-(Subscript[x, i - 1])| = "", 
        Abs[b - a]]; 
      If[Abs[b - a] < difmax, 
        Print[\*""La solución aproximada es (Subscript[x, i]) = "", \ 
          InputForm[b]]; 
        Print["La diferencia entre las dos últimas aproximaciones es "\ \ , Abs[b - a]]; 
        Print["El valor que toma la función en la solución aproximada es "\ \ , p[b]]; Break[]]; \ a = b]\)

Utilizando el método de Newton-Raphson modificado para raíces dobles:

a = -2. ; 
    difmax = (10^(-4)); 
\(Print[\ "i=", 
      1, \*""  \!\(x\_1\) = "", a, \ \*""  p(\!\(x\_1\)) = "", 
      p[a]]; \)
For[i = 2,i ≤ 20,i++,b=a-2p[a]/p'[a];\(Print[\ "i=", i, \*""  \!\(x\_i\) = "", 
            InputForm[b], \ \*""  p(\!\(x\_i\)) = "", 
            p[b], \ \*""  |\!\(x\_i\)-\!\(x\_\(i - 1\)\)| = "", 
            Abs[b - a]]\);\(If[Abs[b - a] < difmax, 
            Print[\*""La solución aproximada es (Subscript[x, i]) = "", \ InputForm[b]]; 
            Print["La diferencia entre las dos últimas aproximaciones es "\ \ , Abs[b - a]]; 
            Print["El valor que toma la función en la solución aproximada es "\ \ , p[b]]; Break[]]\); a = b]

Si en el programa anterior añado una condición de parada como la mencionada en la segunda parte de esta sección, sorprendentemente la actual solución aproximada "pasa" a ser exacta (pero no lo es). Veámoslo:

a = -2. ; 
    difmax = (10^(-4)); 
\(Print[\ "i=", 
      1, \*""  \!\(x\_1\) = "", a, \ \*""  p(\!\(x\_1\)) = "", 
      p[a]]; \)
For[i = 2,i ≤ 20,i++,b=a-2p[a]/p'[a];If[p[b] == 0, 
              Print["La solucion exacta es ", InputForm[b]]; Break[]];\(Print[\ "i=", 
            i, \*""  \!\(x\_i\) = "", 
            InputForm[b], \ \*""  p(\!\(x\_i\)) = "", 
            p[b], \ \*""  |\!\(x\_i\)-\!\(x\_\(i - 1\)\)| = "", 
            Abs[b - a]]\);\(If[Abs[b - a] < difmax, 
            Print[\*""La solución aproximada es (Subscript[x, i]) = "", \ InputForm[b]]; 
            Print["La diferencia entre las dos últimas aproximaciones es "\ \ , Abs[b - a]]; 
            Print["El valor que toma la función en la solución aproximada es "\ \ , p[b]]; Break[]]\); a = b]

Utilizando el método de Newton-Raphson modificado para raíces triples:

a = -2. ; 
    difmax = (10^(-4)); 
\(Print[\ "i=", 
      1, \*""  \!\(x\_1\) = "", a, \ \*""  p(\!\(x\_1\)) = "", 
      p[a]]; \)
For[i = 2,i ≤ 20,i++,b=a-3p[a]/p'[a];\(Print[\ "i=", i, \*""  \!\(x\_i\) = "", 
            InputForm[b], \ \*""  p(\!\(x\_i\)) = "", 
            p[b], \ \*""  |\!\(x\_i\)-\!\(x\_\(i - 1\)\)| = "", 
            Abs[b - a]]\);\(If[Abs[b - a] < difmax, 
            Print[\*""La solución aproximada es (Subscript[x, i]) = "", \ InputForm[b]]; 
            Print["La diferencia entre las dos últimas aproximaciones es "\ \ , Abs[b - a]]; 
            Print["El valor que toma la función en la solución aproximada es "\ \ , p[b]]; Break[]]\); a = b]

Como la convergencia es más lenta que en el caso anterior, podemos sospechar que la raíz es doble.

Como se vio en la segunda parte de esta práctica, esa raíz es exactamente igual a -(Sqrt[5]) y es, efectivamente, doble. Comparamos la solución aproximada obtenida en el segundo caso con la verdadera solución:

InputForm[-(Sqrt[5])-(-2.2360679780560986)]

Resuelva el ejercicio siguiente

Ejercicio

Dada la ecuación ln x + (2/((x^2) + 1)) - 1 = 0, compruebe que el método de Newton-Raphson no converge para el valor inicial (Subscript[x, 1]) = 40. Encuentre otro valor inicial para el que dicho método converja. Proporcione un método que acelere la convergencia a la raíz. Determine la multiplicidad de la raíz.

Solución

f[x_] := Log[x] + (2/((x^2) + 1)) - 1

Método de Newton-Raphson para sistemas

Ejemplo

Aplicando el método de Newton-Raphson para sistemas, calcule de forma aproximada la solución del sistema de ecuaciones
   				 12 ((x^2) - (y^2)) - 24  x + 7 = 0
9 (x - 1) y + 1 = 0
contenida en el intervalo [0, 2/5]×[0, 
        2/5], de modo que la distancia de las dos últimas aproximaciones sea menor que (10^(-8)). Obtenga del mismo modo la aproximación de otra solución real del sistema.

Solución

Se introduce el primer miembro de las ecuaciones (el segundo es 0):

f[x_, y_] := 
      12 ((x^2) - (y^2)) - 24 x + 7; 
g[x_, y_] := 9 (x - 1) y + 1;

Se obtiene el jacobiano de f y g:

J[x, y]=(({{(∂_x) f[x, y], (∂_y) f[x, y]}, {(∂_x) g[x, y], (∂_y) g[x, y]}}));MatrixForm[J[x, y]]

J[x,y] no es una función:

J[2, 3]

Para que lo sea utilizaremos la instrucción
	expresión /. asignaciones
que a cierta parte o partes de la expresión se les asigna cierto valor o variable mediante el símbolo  -> . Por ejemplo:

x + y /. x -> 2       

    u v /. {u -> 2, v -> z}

Así, J se transforma en una función del siguiente modo:

J[a_, b_] := J[x, y] /. {x -> a, y -> b}

Ahora sí podemos calcular valores de la función J(x, y):

MatrixForm[J[2, 3]]

En la iteración del método de Newton-Raphson para sistemas se requiere resolver sistemas del siguiente tipo:

Clear[a, b]; 
    Solve[J[a, b] . {u, v} == {-f[a, b], -g[a, b]}, {u, 
        v}]

Para asignar a u y v los nuevos valores, necesitamos que desaparezca una llave por cada lado, lo que se consigue eligiendo el primer elemento de la lista (que sólo esta formada por un elemento):

%[[1]]

Luego, resumiendo, el problema propuesto puede resolverse con las instrucciones anteriores, añadiendo el punto de partida, la condición de parada, el número máximo de iteraciones y el programa que muestra cada una de las iteraciones del algoritmo:

f[x_, y_] := 12 ((x^2) - (y^2)) - 24 x + 7; 
    g[x_, y_] := 9 (x - 1) y + 1; 
J[x, y]=(({{(∂_x) f[x, y], (∂_y) f[x, y]}, {(∂_x) g[x, y], (∂_y) g[x, y]}}));J[a_, b_] := J[x, y] /. {x -> a, y -> b};{a, b} = {0. , 0. };difmax = (10^(-8));
\(Print[\ "i=", 
      1, \*""  \!\(x\_1\) = "", a, \ \*""  \!\(y\_1\) = "", 
      b, \*""  f(\!\(x\_1\),\!\(y\_1\)) = "", 
      f[a, b], \*""  g(\!\(x\_1\),\!\(y\_1\)) = "", 
      g[a, b]]; \)
\(For[i = 2, 
      i ≤ 
        20, \(i++\), {c, 
          d} = {a, b} + {u, v} /. \(Solve[
              J[a, b] . {u, v} == {\(-f[a, b]\), \(-g[a, b]\)}, {u, 
                v}]\)[\([\)\(1\)\(]\)]; 
      Print[\ "i=", i, \*""  (Subscript[x, i]) = "", 
        c, \ \*""  (Subscript[y, i]) = "", 
        d, \*""  f((Subscript[x, i]),(Subscript[y, i])) = "", 
        f[c, d], \*""  g((Subscript[x, i]),(Subscript[y, i])) = "", 
        g[c, d], \*""  ||((Subscript[x, i]),(Subscript[y, i]))-((Subscript[x, i - 1]),(Subscript[y, i - 1]))|| = "", Norm[{c - a, d - b}]]; 
If[
        Norm[{c - a, d - b}] < difmax, 
        Print[\*""La solución aproximada es ((Subscript[x, i]),(Subscript[y, i])) = "", \ InputForm[{c, d}]]; 
        Print["La distancia entre las dos últimas aproximaciones es "\ \ , Norm[{c - a, d - b}]]; 
        Print[\*""El valor que toma la función en la solución aproximada es (f((Subscript[x, i]),(Subscript[y, i])),g((Subscript[x, i]),(Subscript[y, i]))) = ""\ \ , {f[c, d], g[c, d]}]; Break[]]; \ {a, b} = {c, d}]\)

Luego la solución aproximada obtenida está en el conjunto pedido. Parece que la solución exacta correspondiente es (1/3, 1/6). Efectivamente:

{f[(1/3), (1/6)], g[(1/3), (1/6)]}

Busquemos una segunda solución real, partiendo de algún punto tal que el algoritmo anterior converja a otro punto distinto de (1/3, 1/6):

f[x_, y_] := 12 ((x^2) - (y^2)) - 24 x + 7; 
    g[x_, y_] := 9 (x - 1) y + 1; 
J[x, y]=(({{(∂_x) f[x, y], (∂_y) f[x, y]}, {(∂_x) g[x, y], (∂_y) g[x, y]}}));J[a_, b_] := J[x, y] /. {x -> a, y -> b};{a, b} = {5. , -5. };difmax = (10^(-8));
\(Print[\ "i=", 
      1, \*""  \!\(x\_1\) = "", a, \ \*""  \!\(y\_1\) = "", 
      b, \*""  f(\!\(x\_1\),\!\(y\_1\)) = "", 
      f[a, b], \*""  g(\!\(x\_1\),\!\(y\_1\)) = "", 
      g[a, b]]; \)
\(For[i = 2, 
      i ≤ 
        20, \(i++\), {c, 
          d} = {a, b} + {u, v} /. \(Solve[
              J[a, b] . {u, v} == {\(-f[a, b]\), \(-g[a, b]\)}, {u, 
                v}]\)[\([\)\(1\)\(]\)]; 
      Print[\ "i=", i, \*""  (Subscript[x, i]) = "", 
        c, \ \*""  (Subscript[y, i]) = "", 
        d, \*""  f((Subscript[x, i]),(Subscript[y, i])) = "", 
        f[c, d], \*""  g((Subscript[x, i]),(Subscript[y, i])) = "", 
        g[c, d], \*""  ||((Subscript[x, i]),(Subscript[y, i]))-((Subscript[x, i - 1]),(Subscript[y, i - 1]))|| = "", Norm[{c - a, d - b}]]; 
If[
        Norm[{c - a, d - b}] < difmax, 
        Print[\*""La solución aproximada es ((Subscript[x, i]),(Subscript[y, i])) = "", \ InputForm[{c, d}]]; 
        Print["La distancia entre las dos últimas aproximaciones es "\ \ , Norm[{c - a, d - b}]]; 
        Print[\*""El valor que toma la función en la solución aproximada es (f((Subscript[x, i]),(Subscript[y, i])),g((Subscript[x, i]),(Subscript[y, i]))) = ""\ \ , {f[c, d], g[c, d]}]; Break[]]; \ {a, b} = {c, d}]\)

Parece que la solución exacta correspondiente es (5/3, -1/6). Efectivamente:

{f[(5/3), -(1/6)], g[(5/3), -(1/6)]}

Comentario

La distancia entre las aproximaciones obtenidas en el ejercicio y las correspondientes soluciones exactas son prácticamente nulas:

Norm[{0.33333333333333326,  0.16666666666666666} - {1/3, 
            1/6}]
Norm[{1.6666666666666665,  -0.16666666666666666} - {5/
              3, -1/6}]

Podemos aproximar también la solución mediante el comando directo FindRoot:

FindRoot[{12 ((x^2) - (y^2)) - 24 x + 7, 
        9 (x - 1) y + 1}, {x, 0. }, {y, 0. }, 
      WorkingPrecision -> 20]

FindRoot[{12 ((x^2) - (y^2)) - 24 x + 7, 
        9 (x - 1) y + 1}, {x, 5. }, {y, -5. }, 
      WorkingPrecision -> 20]

Ejemplo

Aplicando el método de Newton-Raphson para sistemas, calcule de forma aproximada dos soluciones reales del sistema de ecuaciones
   				({{(x^2) z - 2  y tg(x z) = 0,}, {2  (e^(z - 1))(x + 3  y) + z = 
                0, 
x - 3  y - 6  z + 4 = 0, }})
de modo que la distancia de las dos últimas aproximaciones de cada solución sea menor que (10^(-6)).

Solución

f[x_, y_, z_] := (x^2) z - 2  y Tan[x z]; 
    g[x_, y_, z_] := 2  (E^(z - 1)) (x + 3  y) + z; 
    h[x_, y_, z_] := x - 3  y - 6  z + 4; 
J[x, y, z]=(({{(∂_x) f[x, y, z], (∂_y) f[x, y, 
                    z], (∂_z) f[x, y, z]}, {(∂_x) g[x, y, z], (∂_y) g[x, y, 
                    z], (∂_z) g[x, y, z]}, {(∂_x) h[x, y, z], (∂_y) h[x, y, 
                    z], (∂_z) h[x, y, z]}}));J[a_, b_, c_] := 
        J[x, y, z] /. {x -> a, y -> b, z -> c};{a, b, c} = {1. , 0. , 0. };difmax = (10^(-6));
\(Print[\ "i=", 
      1, \*""  \!\(x\_1\) = "", a, \ \*""  \!\(y\_1\) = "", 
      b, \ \*""  \!\(z\_1\) = "", 
      c, \*""  f(\!\(x\_1\),\!\(y\_1\),\!\(z\_1\)) = "", 
      f[a, b, c], \*""  g(\!\(x\_1\),\!\(y\_1\),\!\(z\_1\)) = "", 
      g[a, b, c], \*""  h(\!\(x\_1\),\!\(y\_1\),\!\(z\_1\)) = "", 
      h[a, b, c]]; \)
\(For[i = 2, 
      i ≤ 
        20, \(i++\), {a1, b1, 
          c1} = {a, b, c} + {u, v, w} /. \(Solve[
              J[a, b, c] . {u, v, 
                    w} == {\(-f[a, b, c]\), \(-g[a, b, c]\), \(-h[a, b, 
                      c]\)}, {u, v, w}]\)[\([\)\(1\)\(]\)]; 
      Print[\ "i=", i, \*""  (Subscript[x, i]) = "", 
        a1, \ \*""  (Subscript[y, i]) = "", b1, \ \*""  (Subscript[z, i]) = "", 
        c1, \*""  f((Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i])) = "", 
        f[a1, b1, c1], \*""  g((Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i])) = "", 
        g[a1, b1, 
          c1], \ \*""  h((Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i])) = "", 
        h[a1, b1, 
          c1], \*""  ||((Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i]))-((Subscript[x, i - 1]),(Subscript[y, i - 1]),(Subscript[z, i - 1]))|| = "", 
        Norm[{a1 - a, b1 - b, c1 - c}]]; 
If[
        Norm[{a1 - a, b1 - b, c1 - c}] < difmax, 
        Print[\*""La solución aproximada es ((Subscript[x, i]),(Subscript[y, i])) = "", \ InputForm[{a1, b1, c1}]]; 
        Print["La distancia entre las dos últimas aproximaciones es "\ \ , Norm[{a1 - a, b1 - b, c1 - c}]]; 
        Print[\*""El valor que toma la función en la solución aproximada es (f((Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i])),g((Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i])),h((Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i]))) = ""\ \ , {f[a1, 
              b1, c1], g[a1, b1, c1], h[a1, b1, c1]}]; 
        Break[]]; \ {a, b, c} = {a1, b1, c1}]\)

Parece que la solución exacta correspondiente es (-2, 2/3, 0). Efectivamente:

{f[-2, (2/3), 0], g[-2, (2/3), 0], h[-2, (2/3), 0]}

Busquemos una segunda solución real:

f[x_, y_, z_] := (x^2) z - 2  y Tan[x z]; 
    g[x_, y_, z_] := 2  (E^(z - 1)) (x + 3  y) + z; 
    h[x_, y_, z_] := x - 3  y - 6  z + 4; 
J[x, y, z]=(({{(∂_x) f[x, y, z], (∂_y) f[x, y, 
                    z], (∂_z) f[x, y, z]}, {(∂_x) g[x, y, z], (∂_y) g[x, y, 
                    z], (∂_z) g[x, y, z]}, {(∂_x) h[x, y, z], (∂_y) h[x, y, 
                    z], (∂_z) h[x, y, z]}}));J[a_, b_, c_] := 
        J[x, y, z] /. {x -> a, y -> b, z -> c};{a, b, c} = {1. , -1. , -1. };difmax = (10^(-6));
\(Print[\ "i=", 
      1, \*""  \!\(x\_1\) = "", a, \ \*""  \!\(y\_1\) = "", 
      b, \ \*""  \!\(z\_1\) = "", 
      c, \*""  f(\!\(x\_1\),\!\(y\_1\),\!\(z\_1\)) = "", 
      f[a, b, c], \*""  g(\!\(x\_1\),\!\(y\_1\),\!\(z\_1\)) = "", 
      g[a, b, c], \*""  h(\!\(x\_1\),\!\(y\_1\),\!\(z\_1\)) = "", 
      h[a, b, c]]; \)
\(For[i = 2, 
      i ≤ 
        20, \(i++\), {a1, b1, 
          c1} = {a, b, c} + {u, v, w} /. \(Solve[
              J[a, b, c] . {u, v, 
                    w} == {\(-f[a, b, c]\), \(-g[a, b, c]\), \(-h[a, b, 
                      c]\)}, {u, v, w}]\)[\([\)\(1\)\(]\)]; 
      Print[\ "i=", i, \*""  (Subscript[x, i]) = "", 
        a1, \ \*""  (Subscript[y, i]) = "", b1, \ \*""  (Subscript[z, i]) = "", 
        c1, \*""  f((Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i])) = "", 
        f[a1, b1, c1], \*""  g((Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i])) = "", 
        g[a1, b1, 
          c1], \ \*""  h((Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i])) = "", 
        h[a1, b1, 
          c1], \*""  ||((Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i]))-((Subscript[x, i - 1]),(Subscript[y, i - 1]),(Subscript[z, i - 1]))|| = "", 
        Norm[{a1 - a, b1 - b, c1 - c}]]; 
If[
        Norm[{a1 - a, b1 - b, c1 - c}] < difmax, 
        Print[\*""La solución aproximada es ((Subscript[x, i]),(Subscript[y, i])) = "", \ InputForm[{a1, b1, c1}]]; 
        Print["La distancia entre las dos últimas aproximaciones es "\ \ , Norm[{a1 - a, b1 - b, c1 - c}]]; 
        Print[\*""El valor que toma la función en la solución aproximada es (f((Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i])),g((Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i])),h((Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i]))) = ""\ \ , {f[a1, 
              b1, c1], g[a1, b1, c1], h[a1, b1, c1]}]; 
        Break[]]; \ {a, b, c} = {a1, b1, c1}]\)

He encontrado otras tres soluciones reales (y probablemente se puedan encontrar más, quizá muchas más):

f[x_, y_, z_] := (x^2) z - 2  y Tan[x z]; 
    g[x_, y_, z_] := 2  (E^(z - 1)) (x + 3  y) + z; 
    h[x_, y_, z_] := x - 3  y - 6  z + 4; 
J[x, y, z]=(({{(∂_x) f[x, y, z], (∂_y) f[x, y, 
                    z], (∂_z) f[x, y, z]}, {(∂_x) g[x, y, z], (∂_y) g[x, y, 
                    z], (∂_z) g[x, y, z]}, {(∂_x) h[x, y, z], (∂_y) h[x, y, 
                    z], (∂_z) h[x, y, z]}}));J[a_, b_, c_] := 
        J[x, y, z] /. {x -> a, y -> b, z -> c};{a, b, c} = {-1. , -1. , 1. };difmax = (10^(-6));
\(Print[\ "i=", 
      1, \*""  \!\(x\_1\) = "", a, \ \*""  \!\(y\_1\) = "", 
      b, \ \*""  \!\(z\_1\) = "", 
      c, \*""  f(\!\(x\_1\),\!\(y\_1\),\!\(z\_1\)) = "", 
      f[a, b, c], \*""  g(\!\(x\_1\),\!\(y\_1\),\!\(z\_1\)) = "", 
      g[a, b, c], \*""  h(\!\(x\_1\),\!\(y\_1\),\!\(z\_1\)) = "", 
      h[a, b, c]]; \)
\(For[i = 2, 
      i ≤ 
        20, \(i++\), {a1, b1, 
          c1} = {a, b, c} + {u, v, w} /. \(Solve[
              J[a, b, c] . {u, v, 
                    w} == {\(-f[a, b, c]\), \(-g[a, b, c]\), \(-h[a, b, 
                      c]\)}, {u, v, w}]\)[\([\)\(1\)\(]\)]; 
      Print[\ "i=", i, \*""  (Subscript[x, i]) = "", 
        a1, \ \*""  (Subscript[y, i]) = "", b1, \ \*""  (Subscript[z, i]) = "", 
        c1, \*""  f((Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i])) = "", 
        f[a1, b1, c1], \*""  g((Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i])) = "", 
        g[a1, b1, 
          c1], \ \*""  h((Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i])) = "", 
        h[a1, b1, 
          c1], \*""  ||((Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i]))-((Subscript[x, i - 1]),(Subscript[y, i - 1]),(Subscript[z, i - 1]))|| = "", 
        Norm[{a1 - a, b1 - b, c1 - c}]]; 
If[
        Norm[{a1 - a, b1 - b, c1 - c}] < difmax, 
        Print[\*""La solución aproximada es ((Subscript[x, i]),(Subscript[y, i])) = "", \ InputForm[{a1, b1, c1}]]; 
        Print["La distancia entre las dos últimas aproximaciones es "\ \ , Norm[{a1 - a, b1 - b, c1 - c}]]; 
        Print[\*""El valor que toma la función en la solución aproximada es (f((Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i])),g((Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i])),h((Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i]))) = ""\ \ , {f[a1, 
              b1, c1], g[a1, b1, c1], h[a1, b1, c1]}]; 
        Break[]]; \ {a, b, c} = {a1, b1, c1}]\)

f[x_, y_, z_] := (x^2) z - 2  y Tan[x z]; 
    g[x_, y_, z_] := 2  (E^(z - 1)) (x + 3  y) + z; 
    h[x_, y_, z_] := x - 3  y - 6  z + 4; 
J[x, y, z]=(({{(∂_x) f[x, y, z], (∂_y) f[x, y, 
                    z], (∂_z) f[x, y, z]}, {(∂_x) g[x, y, z], (∂_y) g[x, y, 
                    z], (∂_z) g[x, y, z]}, {(∂_x) h[x, y, z], (∂_y) h[x, y, 
                    z], (∂_z) h[x, y, z]}}));J[a_, b_, c_] := 
        J[x, y, z] /. {x -> a, y -> b, z -> c};{a, b, c} = {1. , 1. , -1. };difmax = (10^(-6));
\(Print[\ "i=", 
      1, \*""  \!\(x\_1\) = "", a, \ \*""  \!\(y\_1\) = "", 
      b, \ \*""  \!\(z\_1\) = "", 
      c, \*""  f(\!\(x\_1\),\!\(y\_1\),\!\(z\_1\)) = "", 
      f[a, b, c], \*""  g(\!\(x\_1\),\!\(y\_1\),\!\(z\_1\)) = "", 
      g[a, b, c], \*""  h(\!\(x\_1\),\!\(y\_1\),\!\(z\_1\)) = "", 
      h[a, b, c]]; \)
\(For[i = 2, 
      i ≤ 
        20, \(i++\), {a1, b1, 
          c1} = {a, b, c} + {u, v, w} /. \(Solve[
              J[a, b, c] . {u, v, 
                    w} == {\(-f[a, b, c]\), \(-g[a, b, c]\), \(-h[a, b, 
                      c]\)}, {u, v, w}]\)[\([\)\(1\)\(]\)]; 
      Print[\ "i=", i, \*""  (Subscript[x, i]) = "", 
        a1, \ \*""  (Subscript[y, i]) = "", b1, \ \*""  (Subscript[z, i]) = "", 
        c1, \*""  f((Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i])) = "", 
        f[a1, b1, c1], \*""  g((Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i])) = "", 
        g[a1, b1, 
          c1], \ \*""  h((Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i])) = "", 
        h[a1, b1, 
          c1], \*""  ||((Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i]))-((Subscript[x, i - 1]),(Subscript[y, i - 1]),(Subscript[z, i - 1]))|| = "", 
        Norm[{a1 - a, b1 - b, c1 - c}]]; 
If[
        Norm[{a1 - a, b1 - b, c1 - c}] < difmax, 
        Print[\*""La solución aproximada es ((Subscript[x, i]),(Subscript[y, i])) = "", \ InputForm[{a1, b1, c1}]]; 
        Print["La distancia entre las dos últimas aproximaciones es "\ \ , Norm[{a1 - a, b1 - b, c1 - c}]]; 
        Print[\*""El valor que toma la función en la solución aproximada es (f((Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i])),g((Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i])),h((Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i]))) = ""\ \ , {f[a1, 
              b1, c1], g[a1, b1, c1], h[a1, b1, c1]}]; 
        Break[]]; \ {a, b, c} = {a1, b1, c1}]\)

f[x_, y_, z_] := (x^2) z - 2  y Tan[x z]; 
    g[x_, y_, z_] := 2  (E^(z - 1)) (x + 3  y) + z; 
    h[x_, y_, z_] := x - 3  y - 6  z + 4; 
J[x, y, z]=(({{(∂_x) f[x, y, z], (∂_y) f[x, y, 
                    z], (∂_z) f[x, y, z]}, {(∂_x) g[x, y, z], (∂_y) g[x, y, 
                    z], (∂_z) g[x, y, z]}, {(∂_x) h[x, y, z], (∂_y) h[x, y, 
                    z], (∂_z) h[x, y, z]}}));J[a_, b_, c_] := 
        J[x, y, z] /. {x -> a, y -> b, z -> c};{a, b, c} = {0. , -1. , -1. };difmax = (10^(-6));
\(Print[\ "i=", 
      1, \*""  \!\(x\_1\) = "", a, \ \*""  \!\(y\_1\) = "", 
      b, \ \*""  \!\(z\_1\) = "", 
      c, \*""  f(\!\(x\_1\),\!\(y\_1\),\!\(z\_1\)) = "", 
      f[a, b, c], \*""  g(\!\(x\_1\),\!\(y\_1\),\!\(z\_1\)) = "", 
      g[a, b, c], \*""  h(\!\(x\_1\),\!\(y\_1\),\!\(z\_1\)) = "", 
      h[a, b, c]]; \)
\(For[i = 2, 
      i ≤ 
        20, \(i++\), {a1, b1, 
          c1} = {a, b, c} + {u, v, w} /. \(Solve[
              J[a, b, c] . {u, v, 
                    w} == {\(-f[a, b, c]\), \(-g[a, b, c]\), \(-h[a, b, 
                      c]\)}, {u, v, w}]\)[\([\)\(1\)\(]\)]; 
      Print[\ "i=", i, \*""  (Subscript[x, i]) = "", 
        a1, \ \*""  (Subscript[y, i]) = "", b1, \ \*""  (Subscript[z, i]) = "", 
        c1, \*""  f((Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i])) = "", 
        f[a1, b1, c1], \*""  g((Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i])) = "", 
        g[a1, b1, 
          c1], \ \*""  h((Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i])) = "", 
        h[a1, b1, 
          c1], \*""  ||((Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i]))-((Subscript[x, i - 1]),(Subscript[y, i - 1]),(Subscript[z, i - 1]))|| = "", 
        Norm[{a1 - a, b1 - b, c1 - c}]]; 
If[
        Norm[{a1 - a, b1 - b, c1 - c}] < difmax, 
        Print[\*""La solución aproximada es ((Subscript[x, i]),(Subscript[y, i])) = "", \ InputForm[{a1, b1, c1}]]; 
        Print["La distancia entre las dos últimas aproximaciones es "\ \ , Norm[{a1 - a, b1 - b, c1 - c}]]; 
        Print[\*""El valor que toma la función en la solución aproximada es (f((Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i])),g((Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i])),h((Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i]))) = ""\ \ , {f[a1, 
              b1, c1], g[a1, b1, c1], h[a1, b1, c1]}]; 
        Break[]]; \ {a, b, c} = {a1, b1, c1}]\)

Comentario

Podemos aproximar también la solución mediante el comando FindRoot:

FindRoot[{(x^2) z - 2  y Tan[x z] == 0, 
        2 (E^(z - 1)) (x + 3 y) + z == 0, 
        x - 3 y - 6  z + 4 == 0}, {x, 1. }, {y, 0. }, {z, 0. }, 
      WorkingPrecision -> 20]

FindRoot[{(x^2) z - 2  y Tan[x z] == 0, 
        2 (E^(z - 1)) (x + 3 y) + z == 0, 
        x - 3 y - 6  z + 4 == 0}, {x, 
        1. }, {y, -1. }, {z, -1. }, 
      WorkingPrecision -> 20]

Luego, en este último caso, el mismo punto de partida nos ha llevado a dos soluciones distintas: una con el método de Newton-Raphson programado y otra con la instrucción FindRoot.

FindRoot[{(x^2) z - 2  y Tan[x z] == 0, 
        2 (E^(z - 1)) (x + 3 y) + z == 0, 
        x - 3 y - 6  z + 4 == 0}, {x, -1. }, {y, -1. }, {z, 
        1. }, WorkingPrecision -> 20]

FindRoot[{(x^2) z - 2  y Tan[x z] == 0, 
        2 (E^(z - 1)) (x + 3 y) + z == 0, 
        x - 3 y - 6  z + 4 == 0}, {x, 1. }, {y, 1. }, {z, -1. }, 
      WorkingPrecision -> 20]

De nuevo, el mismo punto de partida nos ha llevado a dos soluciones distintas: una con el método de Newton-Raphson programado y otra con la instrucción FindRoot.

FindRoot[{(x^2) z - 2  y Tan[x z] == 0, 
        2 (E^(z - 1)) (x + 3 y) + z == 0, 
        x - 3 y - 6  z + 4 == 0}, {x, 
        0. }, {y, -1. }, {z, -1. }, 
      WorkingPrecision -> 20]

Ejercicio

Aplicando el método de Newton-Raphson para sistemas, calcule de forma aproximada tres soluciones reales del sistema de ecuaciones
   				 {({{(x^2) sen(x y) - y = 0,}, {2  (y^2) ln(x + 3  y) + x = 0,}})
de modo que la distancia de las dos últimas aproximaciones de cada solución sea menor que (10^(-7)).

Solución

f[x_, y_] := (x^2) Sin[x y] - y; 
    g[x_, y_] := 2  (y^2) Log[x + 3  y] + x

Ejercicio

Obtenga de forma aproximada al menos tres soluciones reales del sistema de ecuaciones
   				 {({{(x^2) + (e^y) = 5,}, {sen x + cos y = 1,}})
de modo que la distancia de las dos últimas aproximaciones de cada solución sea menor que (10^(-6)).

Solución