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()

\sin\left({\left(x + y\right)}^{2}\right)

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')
_images/cell_88_sage0.png
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')

x y e^{\left(x - y\right)}

_images/cell_2_sage01.png

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)
_images/cell_96_sage0.png
sage: implicit_plot(x^2 + 2*y^2,(x,-2,2),(y,-2,2),contours=10).show(aspect_ratio=1)
_images/cell_97_sage0.png
sage: implicit_plot(x^2 + 2*y^2,(x,-2,2),(y,-2,2),contours=srange(0,4,0.1)).show(aspect_ratio=1)
_images/cell_95_sage0.png

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)
_images/cell_92_sage0.png

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()

\log\left(\sin\left(x y\right) + 2\right)

_images/cell_3_sage03.png _images/cell_3_sage1.png _images/cell_3_sage2.png

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:

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')

x \log\left(x^{2} + y^{2}\right)

{\left(\log\left(2\right) + 1\right)} {\left(x - 1\right)} + y + \log\left(2\right) - 1

_images/cell_4_sage01.png

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:

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 x=(x_1,\dots,x_n) es una variable vectorial, y \alpha es un multiíndice. En palabras, el polinomio de Taylor de f de orden k en 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)

{\left(\log\left(2\right) + 1\right)} {\left(x - 1\right)} + y + \log\left(2\right) - 1

{\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()

x e^{\left(x + y\right)}

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')
_images/cell_xxx_sage0.png

Derivadas direccionales

Ejercicio

  • Calcula la derivada direccional de la función f_6 = \sin\left(x y\right) + \cos\left(x y\right) en el punto (1,\pi/2) y la dirección (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()

\sin\left(x y\right) + \cos\left(x y\right)

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()

\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)
_images/cell_16_sage0.png
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)
_images/cell_107_sage0.png

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)
_images/cell_37_sage02.png
sage: #Un algo
sage: region_plot(x^4+x*y+y^2<1,(x,-2,2),(y,-2,2))
_images/cell_45_sage0.png
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))
_images/cell_55_sage0.png

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

\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()
_images/cell_120_sage01.png
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()
_images/cell_124_sage0.png
sage: solve(list(f.gradient()),[x,y])
[2*x + cos(x + y), 4*y^3 + cos(x + y)]

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:

\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:

\left\{x_1\leq x\leq x_2, g(x)\leq y \leq h(x)  \right\}

Ejercicio resuelto

Representa el conjunto de los valores f(x,y) sobre Q=[0,1]\times [0,1] y calcula el volumen del sólido así obtenido.

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')
_images/cell_58_sage01.png
sage: #Un dibujo de curvas de nivel es igualmente ilustrativo
sage: contour_plot(f,(0,1),(0,1),plot_points=300).show()
_images/cell_61_sage0.png
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?
<html>...</html>

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:

\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