Examen de prácticas de Análisis Matemático 
Primero de Informática - septiembre de 2005

Francisco Javier Pérez González
Departamento de Análisis Matemático
Universidad de Granada

Ejercicio 1

Haz un programa para calcular raíces de una función cerca de un determinado punto usando el método de Newton-Raphson. El programa debe admitir como datos de entrada la función, el valor inicial y el error permitido y debe proporcionar la solución encontrada y el número de pasos que han sido necesarios para calcularla. Debes poner una condición para que el programa acabe en un número razonable de pasos. Aplica el programa que has escrito a la función f(x) = 2  x - sen x  - (e^(-x)) con valor inicial (Subscript[x, 0]) = 0 y cota de error (10^(-4)).

Solución. Después de probar diferentes posibilidades, el siguiente programa me parece sencillo y eficaz.
f = función que hay que definir previamente.
x0 = valor inicial.
delta = cota de error de la diferencia entre dos aproximaciones sucesivas.
maxiter = máximo número de iteraciones permitidas.
digitos = número de dígitos con los que trabaja el programa. Mathematica trabaja por defecto con 16. Es conveniente incrementar este número para compensar pérdidas de precisión por redondeos. El valor por defecto de la variable "digitos" es 20. El número de dígitos se usa también para controlar la magnitud del valor de la función en la solución hallada.

NewtonRaphson[f_, x0_, delta_, maxiter_, digitos_:  20] := 
      Module[{a, b}, 
a = 
          SetPrecision[Rationalize[x0], digitos];   
For[k = 0, 
          k <  maxiter, k = k + 1, 
          If[f'[a] ==  0,  
            Print["La derivada se anula. El método no puede \
aplicarse."]; Break[], 
            b = SetPrecision[a - f[a]/(f')[a], 
                digitos]]; 
If[
            SetPrecision[Abs[f[b]], digitos] < (10^(-digitos)),  
            Print["Solución computacionalmente exacta: p = ", 
              b,  "\nCalculada en ",  
              k + 1, " iteraciones.\nf(p)=", 
              SetPrecision[f[b], digitos]]; Break[]]; 
If[
            Abs[b - a] ≤  delta,  
            Print["Solución: p =", b,  "\nCalculada en ",  
              k + 1, " iteraciones.\nDiferencia entre las dos \
últimas aproximaciones menor que ", delta, ".\nf(p)=", 
              SetPrecision[f[b], digitos]]; Break[]]; 
a = 
            b]; k; If[k == maxiter, 
          Print["Se alcanzó el límite de iteraciones sin \
encontrar una solución\nen las condiciones pedidas."]]]

Los siguientes ejemplos ilustran el funcionamiento del programa comparado con el comando FindRoot[ ].

f[x_] = 2  x - Sin[x] - Exp[-x];
NewtonRaphson[f, 0, (10^(-4)), 5]

f[x_] = 2  x - Sin[x] - Exp[-x];
NewtonRaphson[f, 0, (10^(-14)), 5]

Comparemos con el comando FindRoot[ ]. 
Observación. La siguiente celda funciona bien en las versiones de Mathematica anteriores a la 5.0. Si usas la versión 5.0 o posterior debes suprimir la primera y cuarta entrada de la celda y en la segunda entrada debes suprimir "n++;".

n = 0;
sol = 
        FindRoot[n++; f[x], {x, 0}, WorkingPrecision -> 20, 
          AccuracyGoal -> 20];
Print["Solución hallada p = ",  x /. sol]
Print["Número de iteraciones realizadas ",  
      n - 2]
Print["f(p) = ", SetPrecision[f[x /. sol], 20]]

Como puedes ver Mathematica proporciona la misma solución. Probemos ahora con una función muy llana y muy pequeña en el entorno de un cero.

Clear[f]
f[x_] = (((x - 2))^4) Exp[-(x^2)];
NewtonRaphson[f, 1, (10^(-14)), 60, 30]

n = 0;
sol = 
        FindRoot[n++; f[x], {x, 1}, WorkingPrecision -> 30, 
          AccuracyGoal -> 30, MaxIterations -> 100];
Print["Solución hallada p = ",  x /. sol]
Print["Número de iteraciones realizadas ",  
      n - 2]
Print["f(p) = ", f[x /. sol]]

No hay diferencia apreciable entre la solución de nuestro programa y la proporcionada por Mathematica.

Clear[f]
f[x_] = 2  x - 34  Cos[x] - Exp[-x];
NewtonRaphson[f, 0, (10^(-13)), 7]

n = 0;
sol = 
        FindRoot[n++; f[x], {x, 0}, WorkingPrecision -> 20, 
          AccuracyGoal -> 20];
Print["Solución hallada p = ",  x /. sol]
Print["Número de iteraciones realizadas ",  
      n - 2]
Print["f(p) = ", SetPrecision[f[x /. sol], 20]]

En este caso las soluciones son distintas porque Mathematica ha encontrado una solución más próxima al punto inicial x=0  mientras que nuestro programa se ha limitado a aplicar el algoritmo y obtiene una solución distinta.

Plot[f[x], {x, 0, 21}];

Podemos partir en nuestro programa de un punto más cercano a 5.

NewtonRaphson[f, 4, (10^(-14)), 7]

Volvemos a obtener la misma solución proporcionada por Mathematica.

Ejercicio 2

Dado el polinomio

\(TraditionalForm\`p(x) = 
      0.487857\ x\  - \ 0.417177  x\^2 - 
        4.13604  x\^3 + \ \ \ \ 3.0184  x\^4 + 13.0446  x\^5 - 
        13.0802  x\^6 - 
        0.983014  x\^7 + 
		\ \ \ \ \ \ \ \ \ 1.07625  x\^8 - 
        9.44444  x\^9 + 10  x\^10\)

Calcular y clasificar sus puntos críticos y representar su gráfica en un intervalo donde se puedan apreciar dichos puntos.

Solución. Este ejercicio puede hacerse de muchas formas. La que sigue me parece sencilla y elegante.

p[x_] = 
        0.487857 x  -  0.417177  (x^2) - 
          4.13604  (x^3) +   3.01854  (x^4) + 13.0446  (x^5) - 13.0802  (x^6) - 
          0.983014  (x^7) + 1.07625  (x^8) - 9.44444  (x^9) + 
          10  (x^10);
cerosderivada = x /. NSolve[p'[x] == 0, x];
Print["Los ceros de la derivada son ",  cerosderivada]
eligereal[x_] := Im[x] == 0
puntoscriticos = Select[cerosderivada, eligereal];
Print["Los puntos críticos son ",  
      puntoscriticos]
eligeminimo[x_] := p''[x] > 0
minimosrelativos = Select[puntoscriticos, eligeminimo];
eligemaximo[x_] := p''[x] < 0
maximosrelativos = Select[puntoscriticos, eligemaximo];
Print["La función tiene mínimos relativos en los \
puntos ",  minimosrelativos]
Print["La función tiene máximos relativos en los \
puntos ",  maximosrelativos]
cerosrealesderivadasegunda = 
        Select[x /. NSolve[p''[x] == 0, x], 
          eligereal];
eligeinflexion[x_] := Abs[p'''[x]] > 0
puntosdeinflexion = 
        Select[cerosrealesderivadasegunda, eligeinflexion];
Print["La función tiene puntos de inflexión en los \
puntos ",  puntosdeinflexion]
Plot[p[x], {x, Min[puntoscriticos] -  .1, Max[puntoscriticos] +  .1}, 
        AspectRatio -> Automatic, PlotRange -> All];

Llama la atención que haya un mínimo relativo en el punto 0.568732 y un máximo relativo en 0.571274. Veamos la gráfica de la función cerca de estos puntos.

Plot[p[x], {x, 0.565, 0.575}];

Plot[p[x], {x, -0.86, 0.2}];

Plot[p[x], {x, -0.4, 0.2}, PlotRange -> All, 
        AspectRatio ->  .75];