Soluciones de los ejercicios del examen de prácticas de Análisis Matemático de primer curso (grupo A) de Ingeniería Superior de Informática

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 en un intervalo usando el método de bisección. El programa debe admitir como datos de entrada la función, los extremos del intervalo y el error admitido; y debe proporcionar la solución encontrada, el número de pasos que han sido necesarios para calcularla y la cota de error permitida. Debes evitar la posibilidad de que al ejecutar el programa se entre en un ``bucle infinito''.
Aplica a los siguientes casos el programa que has escrito.
a)  f(x)=2  x-sen x - (e^x)   en el intervalo  [0,1]  y cota de error (10^(-4)).
b) g(x) = 
        (x^6) - 3  (x^5) + 2  (x^3) + 3  (x^2) - x - 1  en el intervalo [1.55,2.55]  y cota de error (10^(-4)).

Solución

Este ejercicio es "clásico" y está propuesto en la práctica 5 cuyo título es Resolución de ecuaciones: Métodos de bisección y de Newton-Raphson. Además, dentro de su estilo, es el ejercicio más fácil que se puede poner pues ni siquiera es preciso conocer Mathematica para hacerlo correctamente. Solamente hay que hacer un sencillo bucle con dos condicionamientos y eso puede hacerlo cualquiera que tenga conocimientos básicos de algún lenguaje de programación. Por si eso fuera poco, disponíais de la ayuda de Mathematica y de cualquier información que hubiérais traído al examen (apuntes, ejercicios resueltos de prácticas, tutoriales de Mathematica,...) - dicho sea de paso, el ejercicio está resuelto en un tutorial de Mathematica que está en el directorio de la asignatura.

(* Es evidente que el error que se comete al calcular la raíz es menor que la longitud del intervalo .  Por tanto debemos hacer iteraciones hasta que la longitud del intervalo sea menor que el error permitido. Como en cada iteración la longitud del intervalo se divide por 2, el número de iteraciones necesarias será el más pequeño número natural n tal que (b - a)/(2^n) < error, donde a y b son los extremos iniciales del intervalo .  Por tanto, \(imponiendo\ la\ condición\ indicada\ en\ rojo\ y\ en\ cursiva\ como\ control\ nunca\ podremos\ entrar\ en\ un\ "bucle infinito"\) *)

For[i = 1, b - a ≥  error, i = i + 1, 
            c = (a + b)/2; 
If[f[c] f[a] < 0, b = c, 
              a = c]
];
Print["Solución aproximada: ", 
          c, " con error < ", b - a, "\ncalculada en ", 
          i - 1, " pasos."]

(* hay que poner "i-1 pasos" porque el valor de i fuera del bucle está incrementado en una unidad respecto al valor de i en la última iteración .  Pongo la orden Print[] fuera del bucle porque creo que no tiene ningún interés imprimir el resultado de cada iteración .  Eso solamente sirve para llenar la pantalla de texto cuando lo único que interesa es el resultado final *)

a = N[0];
b = N[1];
error = N[10^-4];
f[x_] = 2  x - Sin[x] - Exp[-x];
For[i = 1, b - a ≥  error, i = i + 1, 
        c = (a + b)/2; 
If[f[c] f[a] < 0, b = c, 
          a = c]];
Print["Solución aproxiamada: ", c, " con error < ", 
      b - a, "\ncalculada en ", i - 1, " pasos."]

a = N[1.55];
b = N[2.55];
error = N[10^-4];
f[x_] = 
        (x^6) - 3  (x^5) + 2  (x^3) + 3  (x^2) - x - 
          1;
For[i = 1, b - a ≥  error, i = i + 1, 
        c = (a + b)/2; 
If[f[c] f[a] < 0, b = c, 
          a = c]];
Print["Solución aproximada: ", c, " con error < ", 
      b - a, "\ncalculada en ", i - 1, " pasos."]

Es inmediato convertir este programa básico en un comando típico de Mathematica. Al hacerlo, conviene añadir una pequeña mejora para controlar si es posible aplicar el método de bisección en el intervalo elegido.

biseccion[f_, xmin_, xmax_, cotadeerror_] := 
      Module[{a = N[xmin], b = N[xmax], 
          error = N[cotadeerror]}, 
If[
          f[a]*f[b] < 0, 
For[i = 1, 
            b - a ≥  error, i = i + 1, 
            c = (a + b)/2; 
If[f[c] f[a] < 0, b = c, 
              a = c]]; 
Print["Solución \
aproximada: ", c, " con error < ", b - a, "\ncalculada en ", 
            i - 1, " pasos."], 
          Print["El método de bisección no puede aplicarse \
en el intervalo elegido."]]
]

Clear[f]
f[x_] = 2  x - Sin[x] - Exp[-x];
g[x_] = 
        (x^6) - 3  (x^5) + 2  (x^3) + 3  (x^2) - x - 
          1;
biseccion[f, 0, 1, 10^-4]
biseccion[g, 1.55, 2.55, 10^-4]
biseccion[f, 1, 2, 10^-4]

Es frecuente añadir a este programa básico algo que no me convence mucho. A saber, la condición If[f[c]==0, Print["Solución exacta: ",N[c,10],"calculada en ",i,"pasos."]. Es casi imposible que calculemos una raíz exacta por el método de bisección pues tendría que ocurrir que dicha raíz coincidiera con uno de los puntos de división del intervalo. La palabra "exacta" en este contexto sólo puede entenderse como "exacta dentro de la precisión de la aritmética de punto flotante" del ordenador que estemos usando. Con este discutible añadido el programa quedaría como sigue.

\(biseccion[f_, xmin_, xmax_, cotadeerror_] := 
      Module[{a = N[xmin], b = N[xmax], 
          error = N[cotadeerror]}, 
If[
          f[a]*f[b] < 0, 
For[i = 1, 
            b - a ≥ \ error, i = i + 1, 
            c = \((a + b)\)/2; 
If[f[c] == 0, 
              Print["Solución exacta: ", c, " calculada en ", 
                i, " pasos."]; 
Break[]]; 
If[f[c] f[a] < 0, b = c, a = c]]; 
If[
            f[c] ≠ 0, 
            Print["Solución aproximada: ", 
              c, " con error < ", b - a, "
calculada en ", 
              i - 1, " pasos."]], 
          Print[\ "El método de bisección no puede aplicarse en el intervalo elegido."]]
]\)

Cosa distinta es que estemos interesados en controlar la magnitud de f[c]. En tal caso podemos poner un control del tipo If[Abs[f[c]]<delta,...] donde "delta" es una cota de error independiente de la cota de error para la raíz aproximada. Podemos entonces interrumpir el bucle cuando se verifique alguno de dichos controles. Observa que podría ocurrir (por ejemplo, para una función muy plana y muy pequeña) que fuera Abs[f[c]]<delta pero c estuviera lejos de una raíz exacta. La condición Abs[f[c]]<delta no sirve para controlar el error que se comete en el cálculo de una raíz exacta. Ahora bien, en un ordenador los cáculos se realizan con una precisión limitada (salvo que trabajemos simbólicamente lo que hace que todo el proceso sea lento) por lo que parece razonable considerar como solución computacionalmente exacta aquella para la cual el valor de f[c] es tan pequeño que 1.0+ f[c] no puede distinguirse de 1.0. Esto ocurre cuando Abs[f[c]]<$MachineEpsilon (recuerda la primera práctica).
A partir de aquí pueden hacerse muchísimas mejoras al programa básico: controlar la precisión en los cálculos internos, incluir como variable un número máximo de iteraciones introducido por el usuario, proporcionar como salida representaciones gráficas, cambiar de intervalo si el que se introduce no está en condiciones de aplicar el método de bisección, y otras muchas que pueden ocurrírsete a ti mismo si lo piensas un poco. Quizás pueda serte útil estudiar el siguiente programa que he hecho. Posiblemente a ti se te ocurran mejoras y simplificaciones que yo no he sabido ver.

Off[General::spell1]
\(Biseccion[func_, u_, v_, error_, digitos_:  16, maxiter_:  Infinity] := 
      Module[{f, maxpasos, a, b, ll, m, M, prec = digitos, delta = error, 
          pasitos = maxiter, c, a2, b2, c2, a1, b1, a0, b0, k, 
          resp}, 
Set @@ {f[x_], func};  (*\ 
          esto\ hace\ que\ el\ programa\ funcione\ tanto\ si\ se\ pone\ la\ función\ en\ la\ forma\ "func[x]"\ o\ simplemente\ "func"\ *) 
a = N[Rationalize[u], prec];  (*\ 
          Rationalize[]\ permite\ controlar\ de\ forma\ precisa\ la\ introducción\ de\ los\ datos\ *) 
b = 
          N[Rationalize[v], prec]; 
maxpasos = 
          Ceiling[N[
              Log[2, \((b - a)\)/delta]]]; 
maxpasos = 
          Min[maxpasos, 
            pasitos]; 
If[
f[a]*f[b] > 
            0, 
m = Min[f[Range[a, b, \((b - a)\)/50]]]; 
          M = Max[f[Range[a, b, \((b - a)\)/50]]]; 
          If[
m*M > 
              0, 
Print["Es posible que la función ", \ 
              f[x], \ " no tenga raíces en el intervalo [", 
              u, ",", v, "]."]; 
            Plot[f[x], {x, a, b}, 
              PlotRange -> {{Min[a, b, 0], Max[a, b, 0]}, {m -  .5, 
                    M +  .5}}, AxesOrigin -> {0, 0}, 
              PlotLabel -> "Gráfica de " <> 
                  ToString[f[x]] <> "
 en el intervalo [" <> 
                  ToString[N[a, 2]] <> ", " <> 
                  ToString[N[b, 2]] <> "]."\ ]; 
            Return[]\ , 
ll = 
              Table[{a + \((k - 1)\) \((b - a)\)/100, 
                  a + k*\((b - a)\)/100}, {k, 1, 
                  100}]; 
{a0, b0} = 
              First[Select[ll, \((f[#[\([1]\)]]*f[#[\([2]\)]] < 0)\) &]]; 
            Pause[ .5]; 
            Plot[f[x], {x, a, b}, 
              PlotRange -> {{Min[a, b, 0], Max[a, b, 0]}, {m -  .5, 
                    M +  .5}}, AxesOrigin -> {0, 0}, 
              PlotLabel -> "Gráfica de " <> 
                  ToString[f[x]] <> "
 en el intervalo [" <> 
                  ToString[N[a, 2]] <> ", " <> 
                  ToString[N[b, 2]] <> "]."\ ]; 
            Print["Para poder aplicar el método se ha cambiado el intervalo
[", u, ",", v, "] por [", N[a0], ",", 
              N[b0], "]."]]
]; 
Which[
          f[a]*f[b] < 0, a1 = a; b1 = b, f[a]*f[b] > 0 && m*M < 0, a1 = a0; 
          b1 = b0]; 
For[
k = 1, 
          k ≤ \ 
            maxpasos, \(k++\), \ 
		\ c = \(a1 + b1\)\/2; \ \ 
		If[
Abs[f[c]] < $MachineEpsilon, 
            Print["Solución computacionalmente exacta: ", 
              NumberForm[c, prec], \ "
calculada en ", \ 
              k, " iteraciones"]; 
Return[
              NumberForm[c, prec]]
		]; 
		If[
            Sign[f[a1]] ≠ Sign[f[c]], b1 = c, 
            a1 = c]; 
		If[
b1 - a1 ≤ \ 
              N[Rationalize[delta], prec], 
            Print["Solución aproximada: ", 
              NumberForm[c, 
                prec]]; 
Print["Número de iteraciones realizadas: ", k]; 
	
            Print["Error máximo cometido: ", 
              N[delta]]; 
Return[
              NumberForm[c, 
                prec]]
		]; 
If[
k + 1 > maxpasos, 
            Print["Se llegó al número máximo de iteraciones permitido ", 
              maxpasos, "
sin encontrar una solución dentro de la cota de error permitida."]; 
Print[\ "La solución aproximada encontrada es ", \ 
              NumberForm[c, 
                prec], \ "
El error máximo cometido es ", 
              b1 - a1]; 
Pause[ .3]; 
resp = 
              Input["¿Quieres una aproximación mejor? Responde: S (si) o N (no)"]; 
            If[resp == S, 
Pause[ .5]; 
              Plot[f[x], {x, a, b}, 
                PlotRange -> {{Min[a, b, 0], Max[a, b, 0]}, All}, 
                AxesOrigin -> {0, 0}]; Pause[ .5]; 
              a2 = Input["A la vista de la gráfica introduce un nuevo valor para el extremo izquierdo del intervalo"]; Pause[ .5]; 
              b2 = Input["A la vista de la gráfica introduce un nuevo valor para el extremo derecho del intervalo"]; 
Biseccion[f[x], a2, b2, delta, prec, Infinity]]; 
            If[resp == N, 
              Return[NumberForm[c, 
                  prec]
	]
	]
	]
	]
	]\)
On[General::spell1]

Puedes comprobar el funcionamiento de este programa con los ejemplos que siguen.

f[x_] := x*Sin[x] - 1

Plot[f[x], {x, -7.5, 7.5}];

Biseccion[f[x], 0, 2, 10^-18, 
        18]//Timing

Biseccion[f[x], 3, 7, 10^-16, 18]

Biseccion[f[x], 3, 7, 10^-16, 20]

Biseccion[f[x], 3, 7, 10^-16, 18, 51]

Biseccion[f[x], 3, 7, 10^-16, 18, 51]

Biseccion[f[x], 3, 6, 10^-10, 18]

Biseccion[f[x], 0, 3, 10^-20, 20]

Ejercicio 2

Define un comando que calcule sumas de Riemann de una función f en un intervalo [a,b]. Dicho comando debe dividir el intervalo [a,b] en n partes iguales, evaluar la función f en el punto medio de cada uno de los subintervalos que resultan, sumar todos los valores así obtenidos de f y multiplicar por la longitud de los subintervalos. Define este comando de manera que trabaje numéricamente.
Aplica a las siguientes funciones el comando que has escrito y compara los valores obtenidos con los que proporciona el comando "NIntegrate[ ]" con las opciones, WorkingPrecision->30, MaxRecursion->30. Debes
conseguir que la diferencia entre ambos resultados sea menor que (10^(-2)). En cada caso indica el número de partes en que divides el intervalo para lograr lo que se pide.
a) f(x) = cos(1/x)  en el intervalo  [0.01,1].
b) g(x)=x sen((e^x)) en el intervalo [0,6].

Solución

Este ejercicio está propuesto, casi con las mismas palabras, en la práctica 6 donde se dan también algunas sugerencias para hacerlo. El ejercicio es tan fácil que lo que me extraña es que no lo hayáis hecho bien todos. Además, al igual que el ejercicio 1,  puede hacerse sin necesidad de saber mucho de Mathematica. Un elemental proceso iterativo es todo lo que se necesita.

(* programa básico *) 
 (* 
      a = N[aquí va el extremo inferior del intervalo] *) 
 (* naturalmente,  
      pongo N[] para asgurarme de que el programa trabajará numéricamente *) 
 (* 
      b = N[aquí va el extremo superior del intervalo] *) 
 (* 
      n =  número de divisiones *) 
 (* f[x_] =  
        aquí va la función *) 
SumaRiemann = 
      Module[{suma = 0, k = 0, inc = N[(b - a)/n]}, 
Do[
          suma = suma + f[a + 0.5*inc + k*inc], {k, 0, 
            n - 1}]; 
suma = 
          suma*inc]

     (* observa que debemos pararnos en n - 
        1 para no salirnos fuera del intervalo *)

a = N[0.01];
b = N[1];
n = 1000;   
f[x_] = Cos[1/x];
Module[{suma = 0.0, k = 0, 
        inc = N[(b - a)/n]}, 
Do[
        suma = suma + f[a + 0.5*inc + k*inc], {k, 0, 
          n - 1}]; 
suma = suma*inc]

NIntegrate[Cos[1/x], {x, 0.01, 1}, WorkingPrecision -> 30, 
      MaxRecursion -> 30]

Abs[% - %%] < 10^-2

a = N[0];
b = N[6];
n = 1000;    
f[x_] = 
      x*Sin[Exp[x]];
Module[{suma = 0.0, k = 0, 
        inc = N[(b - a)/n]}, 
Do[
        suma = suma + f[a + 0.5*inc + k*inc], {k, 0, 
          n - 1}]; 
suma = suma*inc]

NIntegrate[x*Sin[Exp[x]], {x, 0, 6}, WorkingPrecision -> 30, 
      MaxRecursion -> 30]

Abs[% - %%] < 10^-2

Es inmediato convertir en un comando el programita anterior. La única precaución que hay que tomar es no introducir en un proceso iterativo un parámetro formal (como ya sabes de la práctica 4).

SumaRiemann[func_, xmin_, xmax_, n_] := 
      Module[{suma = 0.0, k = 0, a = N[xmin], b = N[xmax], p = n}, 
        inc = N[(b - a)/p]; 
Do[
          suma = suma + func[a + 0.5*inc + k*inc], {k, 0, 
            p - 1}]; 
suma = suma*inc]

f[x_] = Cos[1/x];
SumaRiemann[f, 0.01, 1, 1000]

g[x_] = x*Sin[Exp[x]];
SumaRiemann[g, 0, 6, 1000]

Si usas el comando de Mathematica Sum[ ] puedes ahorrarte trabajo.

nueva  sumariemann[func_, a_, b_, 
          n_] := Sum[
            func[a + 0.5*inc + k*N[(b - a)/n]], {k, 0, 
              n - 1}]*N[(b - a)/n]

nuevasumariemann[g, 0, 6, 1000]
nuevasumariemann[f, 0.01, 1, 1000]

Criterios de calificación

Para calificar he seguido los criterios que se dieron por escrito (aquí no hay sorpresas). Primero: los programas básicos y que funcionen. Quien los haya hecho correctamente está aprobado. Cuando en dichos programas se ha introducido alguna mejora o se han convertido en comandos de Mathematica la calificación ha sido de notable o, en algunos pocos casos, de sobresaliente. He quitado puntos a quienes han cometido fallos imperdonables de sintaxis en la definicón de las funciones, como escribir "senx" en vez de Sin[x] y cosas así. También he quitado puntos a quienes han definido comandos que no se evalúan numéricamente y no han tenido la honradez de advertirlo. Me parece una falta grave entregar un examen de prácticas sabiendo que los comandos definidos no funcionan correctamente (porque eso se sabe) y no informar de ello a quien debe calificarlo.