/* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/
/* [ Created with wxMaxima version 0.8.7 ] */

/* [wxMaxima: title   start ]
Práctica 7. Integrales. Integración numérica.
   [wxMaxima: title   end   ] */

/* [wxMaxima: comment start ]
En esta práctica estudiaremos como calcular integrales definidas e indefinidas. Además
veremos algunas formas de integración numérica mediante fórmulas elementales y también 
mediante comandos específicos de Maxima.
   [wxMaxima: comment end   ] */

/* [wxMaxima: section start ]
Cálculo de primitivas (integral indefinida).
   [wxMaxima: section end   ] */

/* [wxMaxima: comment start ]
Para calcular una primitiva de una función usamos el comando integrate con la siguiente
sintaxis 
integrate(función, variable) .
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
integrate(cos(2*x),x);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Comprobamos que al derivar la primitiva se obtiene la función original.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
diff(sin(2*x)/2,x);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Ejercicio. Halle una primitiva de la función tangente.
(Sol. log(sec(x)) .)
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
El comando integrate no incluye la contante de integración, con lo que nos devuelve una 
integral indefinida, salvo una constante.
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
Las funciones obtenidas son  primitivas de las funciones introducidas, pero no se dice en qué 
intervalo. Como ejemplo, sabemos que la función log(-x) es primitiva de la función 1/x en el 
intervalo de los números estrictamente negativos. Sin embargo la primitiva de 1/x que 
proporciona integrate sólo es valida en el intevalo de los números estrictamente positivos.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
f(x):=log(-x);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
diff(f(x),x);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
integrate(1/x,x);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
A continuación vemos un ejemplo en el que dos primitivas de una misma función se expresan de 
forma diferente, son funciones diferentes y se aprecia la relevancia de la constante de 
integración. 
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
integrate(sin(2*x),x);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Así la función -cos(2*x)/2 es una primitiva de sin(2*x). También la función (sin(x))^2 es una 
primitiva de sin(2*x), lo comprobamos a continuación con el comando diff.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
diff( (sin(x))^2, x);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
trigreduce(  2*cos(x)*sin(x) );
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Las dos primitivas se diferencian en una constante.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
trigreduce(  -cos(2*x)/2-(sin(x))^2  );
/* [wxMaxima: input   end   ] */

/* [wxMaxima: section start ]
Integral definida.
   [wxMaxima: section end   ] */

/* [wxMaxima: comment start ]
El comando integrate tambien nos permite calcular integrales indefinidas si utilizamos
la sintaxis 
integrate(función, variable, a, b) donde a y b son los límites de integración.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
integrate (1-x^2,x,-1,1);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Podemos poner los límites de integración al revés, con lo que la integral cambia de signo.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
integrate(1-x^2,x,1,-1);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Ejercicio. Halle la integral definida de la función sqrt(1-x^2) entre x=-1 y x=1 .
(Esta integral halla el área de media circunferencia de radio 1)
(Sol. pi/2 ).
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
Vemos ahora un ejemplo en el que se halla el área de un recinto plano, concretamente el área 
que delimitan la función seno entre 0 y 2pi, y el eje horizontal. 
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
wxplot2d(sin(x) ,[x,0,2*%pi]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Este área no puede calcularse con la siguiente instrucción ya que una parte de la función es
negativa.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
integrate(sin(x),x,0,2*%pi);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Una manera de plantear el cálculo del área pedida consiste en considerar la función en valor 
absoluto.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
wxplot2d(abs(sin(x)),[x,0,2*%pi]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Aunque el planteamiento es correcto, Maxima devuelve un resultado incorrecto de la integral.
(El valor correcto de la integral definida es 4).
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
integrate(abs(sin(x)),x,0,2*%pi);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Para obtener el resultado correcto debemos conocer en qué intervalos la función es positiva y 
en qué intervalos es negativa.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
integrate(sin(x),x,0,%pi)+integrate(-sin(x),x,%pi,2*%pi);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Ejercicio. Halle el área que delimitan la función x^2+x-1, el eje horizontal y las rectas
x=-2 y x=2.
(La siguiente instrucción, que no es necesario que estudie, dibuja el área pedida)
(Sol. (5^(3/2)+4)/3 = 5.060113295832984 )
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
 wxdraw2d(filled_func=0,xaxis=true,explicit(x^2+x-1,x,-2,2));
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
El comando integrate tampoco nos permite hallar la integral de una función a trozos.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
h(x):=if x<0 then x+2 else -x+2;
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
wxplot2d(h(x),[x,-2,2]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
integrate(h(x),x,-2,2);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Debemos hacer por separado la integral de cada trozo.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
integrate(x+2,x,-2,0)+integrate(-x+2,x,0,2);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Ejercicio. El arco gótico 'quinto acuto' está formado por dos arcos de circunferencia cuyo 
radio es cuatro quintos de la luz del arco. La siguiente función se ha construido siguiendo 
esta descripción. Halle el área que delimitan la función y el eje horizontal, esto es la 
integral definida de f(x) entre 0 y 5.
(Sol. exacta -16*asin(3/8)+8*%pi-(3*sqrt(55))/4 ,
aproximada 13.42024397146637 )
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
f(x):=if x<5/2 then sqrt(16-(x-4)^2) else sqrt(16-(x-1)^2);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
wxplot2d(f(x),[x,0,5], [gnuplot_preamble, "set size ratio -1"]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: section start ]
Fórmulas elementales de integración numérica.
   [wxMaxima: section end   ] */

/* [wxMaxima: input   start ] */
integrate(sqrt(1+x^3),x,1,3);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
integrate(exp(-x^2),x,0,4);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Hay ocasiones en las que una integral definida no puede hallarse de forma exacta, por ejemplo
porque no se puede (o no sabemos) hallar una primitiva, o porque sólo se conoce el valor de la
función en algunos puntos. En estas ocasiones resulta conveniente usar algún método para 
aproximar el valor de la integral definida.
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
Las fórmulas de integración numérica que vamos a ver en esta sección nos permiten aproximar el
valor de una integral definida sumando áreas de figuras geométricas sencillas (rectángulos o 
trapecios).

Para efectuar estas sumas emplearemos el comando sum de Maxima. La sintaxis es
sum(expresión, índice, indice_inicial, indice_final ). Donde la expresión depende de la
variable índice y toma los valores desde indice_inicial hasta indice_final. Vemos a 
continuación un par de ejemplos.
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
Si queremos sumar 1+2+3+4+5+6+7+8+9+10, pediremos a Maxima que sume los valores de una 
variable i, de forma que i=1,2,...,10 como en la siguiente instrucción.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
sum( i , i,0,10);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Para sumar los cuadrados de los 10 primeros números naturales, tendremos que sumar los valores
de la expresión i^2, donde i=1,2,...,10 .
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
sum( i^2 , i,1,10);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Ejercicio: Sume los valores de los cubos de los 20 primeros números naturales.
(Sol. 44100 )
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
Las fórmulas de integración numérica que vamos a describir en esta práctica dividen el 
intervalo [a,b] en n subintervalos iguales.

Por ejemplo si queremos hallar la integral definida en el intervalo [1,3] y dividimos este 
intervalo en 50 subintervalos, definimos las siguientes variables.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
a:1.0;
n:50;
b:3.0;
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Con estos datos, cada uno de los 50 subintervalos tiene longitud (b-a)/n .
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
h:(b-a)/n;
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Los extremos de estos subintervalos (que llamamos nodos) pueden describirse de la siguiente 
forma:
a, a+h, a+2h, a+3h, ..., a+(n-1)h, a+nh=b .

En el ejemplo concreto que estamos considerando, los valores son:
1, 1.04, 1.08, 1.12, ..., 2.96, 3 .
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
Con estos nodos, el primer subintervalo es [a, a+1*h], el segundo es [a+1*h, a+2*h], el 
tercero es [a+2*h, a+3*h] y así el intervalo número i, el i-ésimo es [a+(i-1)*h, a+i*h].
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
Consideramos el rectángulo que tiene por base el intervalo [a+(i-1)*h, a+i*h] y por altura
la que indica la función evaluada en el nodo izquierdo, esto es, f(a+(i-1)*h) .

El área de este rectángulo es base * altura = h * f(a+(i-1)*h) .
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
Si sumamos el área de todos estos rectángulos tenemos una aproximación de la integral definida
que llamamos fórmula de los extremos izquierdos.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
f(x):=sqrt(1+x^3);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
sum( h*f(a+(i-1)*h)  , i , 1 , n );
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Ejercicio. Aproxime la integral definida de la función sqrt(1+x^3) entre x=0.0 y x=4.0
mediante la fórmula de los extremos izquierdos, considerando 10, 25, 50 y 100 subintervalos.
(Sol. 12.61024469701097, 13.42426912125635, 13.70199388697196, 13.84204814533432 )
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
Si en la fórmula que hemos visto, hubiéramos considerado como altura del rectángulo el valor 
de la función en el nodo derecho, esto es, f(a+i*h), el área de cada rectángulo sería 
base*altura = h*f(a+i*h), y tendríamos la fórmula de los extremos derechos.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
a:1.0;
n:50;
b:3.0;
h:(b-a)/n;
f(x):=sqrt(1+x^3);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
sum( h*f(a+i*h) , i , 1 , n );
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Ejercicio. Aproxime la integral definida de la función exp(-x^2) entre x=0.0 y x=4.0
mediante la fórmula de los extremos derechos, considerando 10, 25, 50 y 100 subintervalos.
(Sol. 0.68622692384901, 0.80622691891791, 0.84622691581377, 0.86622691392042 )
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
Si en las fórmulas hubiéramos considerado como altura de los rectángulos el valor de la 
función en el punto medio del subintervalo, esto es, f(a+(i-1/2)*h) tendríamos la fórmula del 
punto medio o regla del punto medio.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
a:1.0;
n:50;
b:3.0;
h:(b-a)/n;
f(x):=sqrt(1+x^3);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
sum( h*f(a+(i-1/2)*h) , i , 1 , n );
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
La cuarta y última fórmula de la sección considera en cada subintervalo el área del trapecio 
formado por el nodo a+(i-1)*h, el nodo a+i*h y el valor de la función en estos nodos, esto es 
f(a+(i-1)*h) y f(a+i*h). El área de este trapecio puede calcularse como 
base*(semisuma de las alturas), esto es h*( f(a+(i-1)*h)+f(a+i*h) )/2 .
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
Puede probarse que, si sumamos todas estas áreas, resulta la siguiente fórmula.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
(h/2)*( f(a)+f(b)+2*sum( f(a+i*h),i,1,n-1 ) );
/* [wxMaxima: input   end   ] */

/* [wxMaxima: section start ]
Comandos para integración numérica.
   [wxMaxima: section end   ] */

/* [wxMaxima: comment start ]
En esta sección veremos dos comandos de Maxima que permiten hallar integrales definidas
usando métodos de integración numérica.
   [wxMaxima: comment end   ] */

/* [wxMaxima: subsect start ]
El comando quad_qags.
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
El comando quad_qags emplea un algoritmo adaptativo para hallar una integral definida.
Su sintaxis es 
quad_qags (función, variable, a, b) donde a y b son los límites de integración.
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
Por ejemplo, para hallar la integral definida de la función exp(-x^2) entre x=1 y x=3
usamos la siguiente instrucción.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
quad_qags( exp(-x^2) , x, 1, 3 );
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
La salida de la instrucción contiene cuatro números. El primero es el valor que buscábamos, 
que es la aproximación a la integral definida. El segundo es una estimación del error absoluto
que se comete (vemos que es un número pequeño). El tercero es el número de evaluaciones que se
han hecho del integrando. El cuarto es un código de error, que si vale cero indica que no ha 
habido ningún error (otros códigos pueden consultarse en la ayuda de Maxima).  
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
Ejercicio. Halle la integral definida de la función sin(x)/x entre x=%pi y x=2*%pi.
a) usando el comando integrate .
b) usando el comando quad_qags .
(Sol. a) ..., b) -0.43378547584984 )
   [wxMaxima: comment end   ] */

/* [wxMaxima: subsect start ]
El comando romberg.
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
El método de Romberg permite también hallar una aproximación a la integral definida de una
función en un intervalo. La sintaxis es similar al comando anterior,
romberg (función, variable, a, b) donde a y b son los límites de integración.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
romberg ( sin(x)/x , x, %pi, 2*%pi );
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Este comando sólo devuelve un valor, que es una aproximación al valor de la integral definida.
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
Ejercicio. Usando el comando romberg aproxime el valor de la integral definida de
la función 1/sqrt(1-x^2) entre -0.99 y 0.99 .
(Sol. 2.85851947267857 )
   [wxMaxima: comment end   ] */

/* [wxMaxima: section start ]
Ejercicios
   [wxMaxima: section end   ] */

/* [wxMaxima: comment start ]
Ejercicio. Defina a1, a2 y g(x) como sigue. Compruebe que g(x) es una primitiva de 
1/(x^4+3*x^2+1). Intente hallar una primitiva de 1/(x^4+3*x^2+1) usando el comando integrate.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
a1:sqrt(2)/sqrt(3-sqrt(5));
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
a2:sqrt(2)/sqrt(3+sqrt(5));
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
g(x):=a1/sqrt(5)*atan(a1*x)-a2/sqrt(5)*atan(a2*x);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Matemáticas II. Grado en Edificación. (A. Palomares.)
   [wxMaxima: comment end   ] */

/* Maxima can't load/batch files which end with a comment! */
"Created with wxMaxima"$