.. -*- coding: utf-8 -*- Cálculo vectorial ================= En esta hoja vamos a estudiar problemas de cálculo diferencial e integral en varias dimensiones espaciales. En general, manejaremos funciones de varias variables simbólicas. Comencemos por representar gráficamente funciones de dos variables reales. :: sage: var('x,y') (x, y) :: sage: f1 = sin((x+y)^2) sage: f1.show() .. MATH:: \sin\left({\left(x + y\right)}^{2}\right) .. index:: plot3d Representación de funciones ::::::::::::::::::::::::::: El comando ``plot3d`` genera una superficie cuya altura sobre el punto (x,y) es proporcional al valor de una función f de dos variables simbólicas. La sintaxis es: :: plot3d(f,(x,x1,x2),(y,y1,y2)) .. - x e y son variables simbólicas - f es función de x y de y - (x1,x2) es el rango que se mostrará de la variable x - (y1,y2) es el rango que se mostrará de la variable y Al mostrar la gráfica llamando a ``show()`` , podemos pasar el argumento opcional viewer que nos permite elegir entre varias formas de mostrar la gráfica, algunas interactivas, que permiten girar la gráfica con el ratón. Como de costumbre, ``plot3d?`` y ``p.show?`` (p es una gráfica generada con ``plot3d`` ) muestran la ayuda con ejemplos de uso. :: sage: p = plot3d(f1,(x,-1,1),(y,-1,1)) sage: p.show(viewer='tachyon') .. image:: b5s3_media/cell_88_sage0.png :align: center :: sage: f2 = x*y*exp(x-y) sage: f2.show() sage: p = plot3d(f2,(x,-2,0.5),(y,-0.5,2)) sage: p.show(viewer='tachyon') .. MATH:: x y e^{\left(x - y\right)} .. image:: b5s3_media/cell_2_sage0.png :align: center .. index:: implicit_plot, contour_plot Curvas de nivel ~~~~~~~~~~~~~~~ Otra opción es dibujar curvas de nivel de la función en el plano (x,y). Usando :: implicit_plot(f,(x,x1,x2),(y,y1,y2)) Dibujamos la curva dada por f=0. Si pasamos a ``implicit_plot`` el argumento adicional ``contours``, el dibujo contiene tantas curvas de nivel como le indiquemos: :: implicit_plot(f,(x,x1,x2),(y,y1,y2), contours= 10) ó exactamente las curvas de nivel que le indiquemos: :: implicit_plot(f,(x,x1,x2),(y,y1,y2), contours = lista) :: sage: #Pasando a show el parametro adicional aspect_ratio=1 sage: #garantizamos que se mantendrá el ratio correcto entre sage: #el eje x y el eje y sage: implicit_plot(x^2 + 2*y^2 - 0.8,(x,-2,2),(y,-2,2)).show(aspect_ratio=1) .. image:: b5s3_media/cell_96_sage0.png :align: center :: sage: implicit_plot(x^2 + 2*y^2,(x,-2,2),(y,-2,2),contours=10).show(aspect_ratio=1) .. image:: b5s3_media/cell_97_sage0.png :align: center :: sage: implicit_plot(x^2 + 2*y^2,(x,-2,2),(y,-2,2),contours=srange(0,4,0.1)).show(aspect_ratio=1) .. image:: b5s3_media/cell_95_sage0.png :align: center ``contour_plot`` es muy similar al anterior, pues muestra de distinto color la región entre dos curvas de nivel consecutivas (es decir, la preimagen de un intervalo). :: sage: contour_plot(x^2+2*y^2,(x,-2,2),(y,-2,2)).show(aspect_ratio=1) .. image:: b5s3_media/cell_92_sage0.png :align: center Comparamos las distintas formas de representar una función de dos variables. :: sage: f3=log(2+sin(x*y)) sage: f3.show() sage: plot3d(f3,(x,-3,3),(y,-3,3)).show(viewer='tachyon') sage: contour_plot(f3,(x,-3,3),(y,-3,3)).show() sage: implicit_plot(f3,(x,-3,3),(y,-3,3),contours=10).show() .. MATH:: \log\left(\sin\left(x y\right) + 2\right) .. image:: b5s3_media/cell_3_sage0.png :align: center .. image:: b5s3_media/cell_3_sage1.png :align: center .. image:: b5s3_media/cell_3_sage2.png :align: center Plano tangente ~~~~~~~~~~~~~~ Una función de dos variables simbólicas se puede derivar respecto de cada una de ellas de forma independiente. De esta forma podemos aplicar la fórmula del plano tangente en un punto: .. MATH:: z = f(x_0, y_0)+\frac{\partial f}{\partial x}(x_0, y_0)(x-x_0) + \frac{\partial f}{\partial y}(x_0,y_0)(y - y_0) Las gráficas de curvas de nivel no permiten apreciar bien el plano tangente, pero podemos representar el plano junto a la función con ``plot3d`` . :: sage: f4=x*log(x^2+y^2) sage: f4.show() sage: x0 = 1 sage: y0 = 1 sage: plano = (x - x0)*f4.derivative(x,1).subs(x=x0, y=y0) + (y - y0)*f4.derivative(y,1).subs(x=x0, y=y0) + f4(x=x0, y=y0) sage: show(plano) sage: p = (plot3d(f4,(x,-1.5,1.5),(y,-1.5,1.5)) + ... plot3d(plano,(x,-1.5,1.5),(y,-1.5,1.5),color='red') + ... point3d( (x0,y0,f4(x=x0,y=y0)), color='green', pointsize='30') ) sage: p.show(viewer='tachyon') .. MATH:: x \log\left(x^{2} + y^{2}\right) .. MATH:: {\left(\log\left(2\right) + 1\right)} {\left(x - 1\right)} + y + \log\left(2\right) - 1 .. image:: b5s3_media/cell_4_sage0.png :align: center .. index:: taylor Polinomio de Taylor ~~~~~~~~~~~~~~~~~~~ La fórmula del plano tangente es el caso particular de la fórmula del polinomio de Taylor de la función en un punto: .. MATH:: f(x) =\sum_{|\alpha|=0}^k\frac{1}{\alpha!}\frac{\partial^\alpha f(a)}{\partial x^\alpha}(x-a)^\alpha+\sum_{|\alpha|=k+1}R_{\alpha}(x)(x-a)^\alpha Donde :math:`x=(x_1,\dots,x_n)` es una variable vectorial, y :math:`\alpha` es un multiíndice. En palabras, el polinomio de Taylor de f de orden k en :math:`a=(a_1,\dots,a_n)` contiene todos los monomios construídos con derivadas de orden menor o igual que k. La sintaxis para construir el polinomio de Taylor de f de orden k en las variables (x,y) y en el punto (a,b) es: :: f.taylor((x,a),(y,b),k) :: sage: #El polinomio de Taylor de orden 1 es la ecuacion del sage: #plano tangente sage: x0 = 1 sage: y0 = 1 sage: plano = (x - x0)*f4.derivative(x,1).subs(x=x0,y=y0) + (y - y0)*f4.derivative(y,1).subs(x=x0,y=y0) + f4(x=x0,y=y0) sage: show(plano) sage: f4t1 = f4.taylor((x,x0),(y,y0),1) sage: show(f4t1) .. MATH:: {\left(\log\left(2\right) + 1\right)} {\left(x - 1\right)} + y + \log\left(2\right) - 1 .. MATH:: {\left(\log\left(2\right) + 1\right)} {\left(x - 1\right)} + y + \log\left(2\right) - 1 :: sage: f5 = x*exp(x + y) sage: f5.show() sage: pol5 = f5.taylor((x,0),(y,0),2) sage: pol5.show() .. MATH:: x e^{\left(x + y\right)} .. MATH:: x^{2} + x y + x :: sage: p1=plot3d(f5,(x,-1,1),(y,-1,1)); sage: p2=plot3d(pol5,(x,-1,1),(y,-1,1),color='red'); sage: show(p1+p2,viewer='tachyon') .. image:: b5s3_media/cell_xxx_sage0.png :align: center .. index:: arrow Derivadas direccionales ~~~~~~~~~~~~~~~~~~~~~~~ Ejercicio --------- - Calcula la derivada direccional de la función :math:`f_6 = \sin\left(x y\right) + \cos\left(x y\right)` en el punto :math:`(1,\pi/2)` y la dirección :math:`(1,1)`. - Investiga la ayuda del comando ``arrow``, y dibuja el gradiente de f6 en un punto. Observa el efecto al cambiar la función, y el punto. :: sage: f6=sin(x*y)+cos(x*y) sage: f6.show() .. MATH:: \sin\left(x y\right) + \cos\left(x y\right) .. index:: plot_vector_field Campo gradiente ~~~~~~~~~~~~~~~ Usando el comando ``plot_vector_field`` , podemos dibujar el **campo gradiente** . El campo de vectores *gradiente de f* es la asignación del vector gradiente de f a cada punto del plano. Recordemos de las clases de teoría que el gradiente es perpendicular a las curvas de nivel. Es una buena ocasión para dibujar simultáneamente las curvas de nivel y el campo gradiente de una función. :: sage: f7=sin(x*y)+cos(x*y) sage: f7.show() .. MATH:: \sin\left(x y\right) + \cos\left(x y\right) :: sage: plot_vector_field(f7.gradient(),(x,-2,2),(y,-2,2)).show(aspect_ratio=1) .. image:: b5s3_media/cell_16_sage0.png :align: center :: sage: (contour_plot(f7,(x,-2,2),(y,-2,2)) + plot_vector_field(f7.gradient(),(x,-2,2),(y,-2,2))).show(aspect_ratio=1) .. image:: b5s3_media/cell_107_sage0.png :align: center .. index:: region_plot Regiones ~~~~~~~~ Otro comando interesante es el comando ``region_plot``, que sirve para dibujar regiones definidas por desigualdades. Se puede llamar de varias formas para representar regiones definidas por varias ecuaciones (ver la ayuda para más información). En esta sesión lo usaremos para visualizar regiones de integración. :: sage: #Una elipse sage: region_plot((x-1)^2 + 2*y^2 < 1,(x,-2,2),(y,-2,2)).show(aspect_ratio=1) .. image:: b5s3_media/cell_37_sage0.png :align: center :: sage: #Un algo sage: region_plot(x^4+x*y+y^2<1,(x,-2,2),(y,-2,2)) .. image:: b5s3_media/cell_45_sage0.png :align: center :: sage: #La interseccion de las dos regiones anteriores sage: region_plot([x^4+x*y+y^2<1, (x-1)^2+2*y^2<1],(x,-2,2),(y,-2,2)) .. image:: b5s3_media/cell_55_sage0.png :align: center Ejercicio --------- dibuja la región 0<=x<=1,x\*(1\-x)<=y<=sin(pi\*x) Puntos críticos y extremos :::::::::::::::::::::::::: Ceros del gradiente ~~~~~~~~~~~~~~~~~~~ La forma más naive de buscar puntos críticos es resolver el sistema de ecuaciones no lineales .. MATH:: \nabla(f)=0 El sistema es trivial de plantear usando el método gradient visto antes, y podemos intentar resolver el sistema llamando a ``solve`` : :: solve(list(f.gradient()),f.variables()) :: sage: #Esta funcion parece tener un minimo cerca de (-0.5,-0.5) sage: var('x y') sage: f = x^2 + y^4 + x + y sage: contour_plot(f,(x,-2,2),(y,-2,2)).show() .. image:: b5s3_media/cell_120_sage0.png :align: center :: sage: #Obtenemos varias soluciones "extra": sage: #?Cual corresponde al minimo? sage: solve(list(f.gradient()),f.variables()) [[x == -0.5, y == (0.314980262474 + 0.545561817986*I)], [x == -0.5, y == (0.314980262474 - 0.545561817986*I)], [x == -0.5, y == -0.629960578187]] Por supuesto, con funciones más complicadas, este tipo de sistemas de ecuaciones no lineales es casi imposible de resolver. :: sage: #Esta funcion parece tener un minimo cerca de (-0.3,-0.5) sage: var('x y') sage: f = x^2+ y^4 + sin(x+y) sage: contour_plot(f,(x,-1,1),(y,-1,1),plot_points=300).show() .. image:: b5s3_media/cell_124_sage0.png :align: center :: sage: solve(list(f.gradient()),[x,y]) [2*x + cos(x + y), 4*y^3 + cos(x + y)] .. index:: scipy.optimize, fmin Mínimos de una función de varias variables obtenidos numéricamente ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ El paquete ``scipy.optimize`` contiene varios métodos relacionados con problemas de optimización. El método ``scipy.optimize.fmin`` busca mínimos de una función de varias variables partiendo de un punto inicial. Para llamar al método, es necesario definir una función de python que acepta como argumento una lista con los valores de las variables y devuelve el valor de la función en el punto. El resultado es el mínimo propuesto, y una estimación del error cometido. Además, nos imprime un resumen de la ejecución del algoritmo numérico. :: sage: from scipy.optimize import fmin sage: def f_0(xs): ... x,y = xs ... return x^2 + y^4 + x + y sage: fmin(f_0,(0,0)) Optimization terminated successfully. Current function value: -0.722470 Iterations: 58 Function evaluations: 113 array([-0.500004 , -0.62998934]) :: sage: from scipy.optimize import fmin sage: def f_0(xs): ... x,y = xs ... return x^2 + y^4 + sin(x + y) sage: fmin(f_0,(0,0)) Optimization terminated successfully. Current function value: -0.570485 Iterations: 60 Function evaluations: 109 array([-0.32318198, -0.54469709]) Integrales :::::::::: La integración simbólica sobre regiones cuadradas no supone ningún quebradero de cabeza (siempre que el sistema sea capaz de encontrar las integrales del integrando). Para calcular la integral doble, integramos primero respecto de una variable y después respecto de la otra. Ejemplo: .. MATH:: \iint_Q x^2\,e^y\,d x\,d y, \ Q=[-1,1]\times [0,\log 2] :: sage: f = x^2*e^y sage: f.integral(x,-1,1).integral(y,0,log(2)) 2/3 Para regiones más generales, tenemos que operar igual que en clase de Cálculo, descomponiendo la región en subregiones del tipo: .. MATH:: \left\{x_1\leq x\leq x_2, g(x)\leq y \leq h(x) \right\} .. index:: scipy.integrate, dblquad Ejercicio resuelto ~~~~~~~~~~~~~~~~~~ Representa el conjunto de los valores :math:`f(x,y)` sobre :math:`Q=[0,1]\times [0,1]` y calcula el volumen del sólido así obtenido. .. MATH:: f(x,y)= \left\{ \begin{array}{ll} x+y &\quad\mbox{si}\quad x^2\le y\le 2 x^2, \\ 0 &\quad\mbox{en el resto}. \end{array}\right. :: sage: #Comenzamos por dibujar la funcion sage: #!!!Atencion!!! Esta vez el argumento de plot3d no es una sage: #funcion de dos variables simbolicas, sino una funcion sage: #normal de python que acepta como argumentos dos numeros sage: #reales y devuelve otro sage: def f(a,b): ... if a^2 <= b <= 2*a^2: ... return a+b ... else: ... return 0 sage: plot3d(f,(0,1),(0,1)).show(viewer='tachyon') .. image:: b5s3_media/cell_58_sage0.png :align: center :: sage: #Un dibujo de curvas de nivel es igualmente ilustrativo sage: contour_plot(f,(0,1),(0,1),plot_points=300).show() .. image:: b5s3_media/cell_61_sage0.png :align: center :: sage: #para calcular la integral simbolica necesitamos definir sage: #una funcion de dos variables simbolicas, pero debemos sage: #integrar solo en la region adecuada sage: var('x y') sage: g = x + y sage: print g sage: print g.integral(y,x^2,2*x^2) sage: print g.integral(y,x^2,2*x^2).integral(x,0,1) x + y 3/2*x^4 + x^3 11/20 Cuando sea imposible evaluar la integral de una función sobre una región básica del plano de forma exacta, podemos hacerlo de forma numérica con ``scipy.integrate.dblquad`` :: sage: #Es necesario importar la funcion dblquad sage: #del paquete scipy.integrate al que pertenece sage: from scipy.integrate import dblquad :: sage: dblquad? ... Calculamos la misma integral que calculamos antes de forma simbólica usando ``dblquad.`` :: sage: def g(x,y): ... return x+y sage: def gfun(x): ... return x^2 sage: def hfun(x): ... return 2*x^2 sage: dblquad(g,0,1,gfun,hfun) (0.54999999999999993, 6.1062266354383602e-15) Cambios de coordenadas ~~~~~~~~~~~~~~~~~~~~~~ Probemos a calcular una integral en coordenadas cartesianas y polares. Recordemos que: .. MATH:: \iint_R f(x,y)dx\,dy = \iint_R rf(r,\theta)dr\,d\theta :: sage: var('x y r theta') (x, y, r, theta) :: sage: f = x^2+y^2 sage: print f.integral(y,-sqrt(1-x^2),sqrt(1-x^2)).integral(x,-1,1) CODE: sage47 : integrate(sage43,sage44,sage45,sage46)$ Maxima ERROR: defint: lower limit of integration must be real; found -sqrt(1-x^2) -- an error. To debug this try: debugmode(true); Traceback (most recent call last): ... TypeError: Error executing code in Maxima El código anterior produce un error. Hurgando un poco en la traza del error, leemos: :: TypeError: Error executing code in Maxima CODE: sage9 : integrate(sage5,sage6,sage7,sage8) Maxima ERROR: defint: lower limit of integration must be real; found \-sqrt(1\-x^2) Recordemos que SAGE usa otras piezas de software con más tradición siempre que puede. En este caso, le pasa la tarea de realizar la integral a la librería **Maxima** . Esta librería arroja un error bastante ilustrativo: "el límite de integración debe ser real", pero sqrt(1\-x^2) puede ser imaginario. Entender la causa del error no siempre garantiza que se pueda superar el problema, pero en este caso la solución es informar a Sage de que puede asumir que x está entre \-1 y 1. En cualquier caso, nada nos garantiza el éxito al intentar una integral de forma simbólica, así que deberíamos estar preparadas para hacer la integral con **dblquad** si todo lo demás falla. :: sage: assume(1-x^2>0) sage: f = x^2+y^2 sage: print f.integral(y,-sqrt(1-x^2),sqrt(1-x^2)).integral(x,-1,1) 1/2*pi :: sage: #Integramos en coordenadas polares. sage: f_polar = f(x=r*cos(theta),y=r*sin(theta)) sage: print f_polar sage: #Intentamos simplificar el integrando (no es necesario) sage: f_polar = f_polar.full_simplify() sage: print f_polar sage: integrando_polar = r*f_polar() sage: integrando_polar.integral(theta,0,2*pi).integral(r,0,1) r^2*sin(theta)^2 + r^2*cos(theta)^2 r^2 1/2*pi