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

Tema 2
Resolución numérica de ecuaciones

Clear["Global`*"]

Método de Newton-Raphson para sistemas

Ejercicio

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], aplicando 5 veces la iteración del método de Newton-Raphson, partiendo del punto (0, 0) y presentando las aproximaciones con 15 cifras decimales.

Utilizando comandos directos de Mathematica

Solve[{12 ((x^2) - (y^2)) - 24 x + 7 == 0, 
        9 (x - 1) y + 1 == 0}, {x, y}]

Programando el método de Newton-Raphson para sistemas

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 define el jacobiano:

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

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

{a, b} = {0. , 0. }; m[0] = {a, b, f[a, b], g[a, b]}; 
    For[i = 1, 
      i ≤ 
        5, i++, {a, 
          b} = {a, b} + {u, 
              v} /. (Solve[
                J[a, b] . {u, v} == {-f[a, b], -g[a, b]}, {u, 
                  v}])[[1]]; m[i] = {a, b, f[a, b], g[a, b]}]; 
    TableForm[
      Table[{i, InputForm[m[i][[1]]], InputForm[m[i][[2]]], 
          InputForm[m[i][[3]]], InputForm[m[i][[4]]]}, {i, 0, 
          5}]]

Resuelva los siguientes ejercicios

Ejercicio 1

Aplicando el método de Newton-Raphson para sistemas, calcule de forma aproximada una solución del sistema de ecuaciones
   				 {({{(x^2) sen(x y) = 0,}, {2  (y^2) ln(x + 3  y) = 0,}})
partiendo del punto (1, 1) y presentando las aproximaciones con 15 cifras decimales. Obseve que en este caso la convergencia a la solución exacta (0, 1/3) es bastante lenta.

Ejercicio 2

Aplicando el método de Newton-Raphson para sistemas, calcule de forma aproximada una solución del sistema de ecuaciones
   				 {({{(x^2) sen(x y) - y = 0,}, {2  (y^2) ln(x + 3  y) + x = 0,}})
partiendo del punto (1, 1) y presentando las aproximaciones con 15 cifras decimales. Parta de otros puntos hasta encontrar aproximaciones de al menos otras dos soluciones del sistema.

Ejercicio 3

Aplicando el método de Newton-Raphson para sistemas, calcule de forma aproximada una solución 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, }})
partiendo del punto (-1, -1, 1) y presentando las aproximaciones con 15 cifras decimales. Parta de otro punto hasta encontrar una aproximación de otra solución del sistema.

Método de Sturm

Ejercicio

Aplicando el teorema de Sturm, localice y separe las raíces reales del polinomio

-500 + 575 x - 275 (x^2) + 70 (x^3) + 45 (x^4) - 72 (x^5) + 
      31 (x^6) + 2 (x^7) - 5 (x^8) + (x^9)

Solución

q[x_] := -500 + 575 x - 275 (x^2) + 70 (x^3) + 45 (x^4) - 
          72 (x^5) + 31 (x^6) + 2 (x^7) - 5 (x^8) + (x^9);

En primer lugar obtenemos el máximo común divisor ("the greatest common divisor", en inglés) del polinomio y de su derivada, para ver si el polinomio tiene raíces múltiples.

m[x] = PolynomialGCD[q[x], q'[x]]; m[x]

Como el máximo común divisor m(x) no es constante, el polinomio q(x) tiene raíces múltiples, que son las raíces de m(x), a saber:

Solve[m[x] == 0, x]

Luego (Sqrt[5]) y -(Sqrt[5]) son raíces dobles de q(x). El cociente ((q(x))/(m(x))) ya tiene todas sus raíces simples (con lo que se le podrá aplicar el teorema de Sturm), y éstas coinciden con las del polinomio de partida q(x).

Simplify[q[x]/m[x]]

Denominamos p(x) al polinomio anterior, que será el primer elemento de  la sucesión de Sturm; el segundo elemento es su derivada; y los siguientes (hasta un máximo de 8) son restos de divisiones (con el signo cambiado). El comando de Mathematica PolynomialRemainder[(Subscript[p, 1]), (Subscript[p, 2]), x] indica el resto de la división del polinomio (Subscript[p, 1])(x) entre el polinomio (Subscript[p, 2])(x). Una vez que el resto es una constante se ha completado la sucesión de Sturm.

p[x_] := 
        100 - 115 x + 75 (x^2) - 37 (x^3) + 6 (x^4) + 7 (x^5) - 5 (x^6) + 
          (x^7);
f[0] = p[x];
f[1] = p'[x];
For[i = 2, i ≤ 7, i++, 
      f[i] = -PolynomialRemainder[f[i - 2], f[i - 1], x]]

st = Table[f[i], {i, 0, 7}]; TableForm[st]

Inciso: Cuando se parte de un polinomio con raíces múltiples (como era el caso de q(x)), el último resto de la sucesión de Sturm es cero, y el último no nulo (el anterior) es el máximo común divisor del polinomio y su derivada (expresado salvo una constante). En nuestro ejemplo ese resto será el polinomio m(x) (expresado salvo una constante). Comprobémoslo:

g[0] = q[x];
g[1] = q'[x];
For[i = 2, i ≤ 8, i++, 
      g[i] = -PolynomialRemainder[g[i - 2], g[i - 1], x]]

TableForm[Table[g[i], {i, 0, 8}]]

Simplify[g[7]/m[x]]

Para localizar las raíces del polinomio p(x) (que son las mismas que las del polinomio de partida q(x)), se aplica el teorema que afirma que todas las raíces (reales e imaginarias) del polinomio con coeficientes reales

(Subscript[a, 0]) + (Subscript[a, 1]) x + (Subscript[a, 2]) (x^2) +  ...  + (Subscript[a, n - 1]) (x^(n - 1)) + 
      (Subscript[a, n]) (x^n)

se encuentran en el círculo de centro (0, 0) y radio ρ + 1, donde

ρ = (1/(|(Subscript[a, n])|)) 
          (max_(0 ≤ i ≤ n - 1)) | (Subscript[a, i]) |  .

En nuestro caso, es claro que ρ=115 y que, por tanto, las raices reales del polinomio p(x) se encuentran en el intervalo [-116, 116].

La sucesión de Sturm debe evaluarse en x = -116 y x = 116 para estudiar la diferencia de cambios de signo en uno y otro punto (nosotros sólo representamos los signos, que es lo único que interesa):

ρ = 115;
x = -ρ - 1;  Sign[st]
x = ρ + 1;   Sign[st]

O de manera más desarrollada:

ρ = 115;
x = -ρ - 1;  
    Table[ Sign[f[i]], {i, 0, 7}]
x = ρ + 1;  Table[ Sign[f[i]], {i, 0, 7}]

Luego las alternancias de signo son 5 en x = -116 y 2 en x = 116. Luego hay 5 - 2 = 3 raíces reales en p(x) (y en en q(x)). Dos de ellas son dobles (las raíces de m(x), esto es, (Sqrt[5]) y -(Sqrt[5])) y hay una tercera simple. Aplicando el teorema de Sturm podemos localizarla:

x = 0;  Sign[st]

Las alternancias de signo en x = 0 son 4. Luego hay 5 - 4 = 1 raíz real negativa en p(x) (-(Sqrt[5])) y 4 - 2 = 2 positivas ((Sqrt[5]) y otra).

x = 3;  Sign[st]

Las alternancias de signo en x = 3 son 2. Luego las 2 raíces positivas de p(x) se sitúan en el intervalo [0, 3]. Una de ellas es aproximadamente igual a

(Sqrt[5.])

x = 2.23;  Sign[st]

Hay 3 alternancias de signo en x = 2.23 < (Sqrt[5]): luego la tercera raíz buscada está en el intervalo [0, 2.23].

x = 2;  Sign[st]

Hay 4 alternancias de signo en x = 2: luego la tercera raíz buscada está en el intervalo [2, 2.23], con lo que podemos dar por localizada la tercera raíz. Ya podrá aproximarse bien por alguno de los métodos generales de obtención de raíces.

Veamos las tres raíces reales sobre la gráfica de la función polinómica:

Plot[q[x],  {x, -3, 3}];

Como no se aprecian bien dos de ellas, nos acercamos un poco:

Plot[q[x],  {x, 2, 2.5}];

Resuelva

Ejercicio

Dado el polinomio p(x) = 
            318 - 540 x - 1423 (x^2) + 2728 (x^3) - 
              863 (x^4) - 130  (x^5) + 106 (x^6) - 18 (x^7) + (x^8), realice lo siguiente:
	a) Calcule el m.c.d de p(x) y de su derivada. Como consecuencia, responda a la siguiente pregunta: ¿tiene p(x) raíces múltiples?
	b) A partir de la expresión de p(x), ¿en qué región del plano podemos asegurar que se encuentran sus raíces reales e imaginarias?
	c) Aplique el teorema de Sturm para localizar todas las raíces reales de p(x) en intervalos disjuntos cuya longitud sea como máximo igual a 1.
	d) Haga el menor número posible de representaciones de la gráfica de p(x) en la que se muestren con claridad las raíces de este polinomio.
	e) Obtenga una aproximación de las raíces reales y complejas de p(x) con el comando NSolve.

Solución

Apartado a)

Apartado b)

Apartado c)

Apartado d)

Apartado e)