Álgebra Y Cálculo Básicos

Sage puede efectuar cómputos relacionados al algebra y cálculo básicos: por ejemplo, encontrar soluciones de ecuaciones, diferenciación, integración y transformadas de Laplace. Véa la documentación «Construcciones En Sage» para más ejemplos.

Resolviendo Ecuaciones

Resolviendo Ecuaciones De Manera Exacta

La función solve resuelve ecuaciones. Para usarla, primero no olvides especificar algunas variables. Los argumentos de solve son una ecuación (o un sistema de ecuaciones), junto con las variables a resolver:

sage: x = var('x')
sage: solve(x^2 + 3*x + 2, x)
[x == -2, x == -1]
>>> from sage.all import *
>>> x = var('x')
>>> solve(x**Integer(2) + Integer(3)*x + Integer(2), x)
[x == -2, x == -1]

Puedes resolver ecuaciones en una variable respecto de las demás:

sage: x, b, c = var('x b c')
sage: solve([x^2 + b*x + c == 0],x)
[x == -1/2*b - 1/2*sqrt(b^2 - 4*c), x == -1/2*b + 1/2*sqrt(b^2 - 4*c)]
>>> from sage.all import *
>>> x, b, c = var('x b c')
>>> solve([x**Integer(2) + b*x + c == Integer(0)],x)
[x == -1/2*b - 1/2*sqrt(b^2 - 4*c), x == -1/2*b + 1/2*sqrt(b^2 - 4*c)]

Puedes también resolver ecuaciones en varias variables:

sage: x, y = var('x, y')
sage: solve([x+y==6, x-y==4], x, y)
[[x == 5, y == 1]]
>>> from sage.all import *
>>> x, y = var('x, y')
>>> solve([x+y==Integer(6), x-y==Integer(4)], x, y)
[[x == 5, y == 1]]

El siguiente ejemplo del uso de Sage para resolver un sistema de ecuaciones no-lineales fue proporcionado por Jason Grout: primero, resolvemos el sistema simbólicamente:

sage: var('x y p q')
(x, y, p, q)
sage: eq1 = p+q==9
sage: eq2 = q*y+p*x==-6
sage: eq3 = q*y^2+p*x^2==24
sage: solve([eq1,eq2,eq3,p==1],p,q,x,y)
[[p == 1, q == 8, x == -4/3*sqrt(10) - 2/3, y == 1/6*sqrt(10) - 2/3], [p == 1, q == 8, x == 4/3*sqrt(10) - 2/3, y == -1/6*sqrt(10) - 2/3]]
>>> from sage.all import *
>>> var('x y p q')
(x, y, p, q)
>>> eq1 = p+q==Integer(9)
>>> eq2 = q*y+p*x==-Integer(6)
>>> eq3 = q*y**Integer(2)+p*x**Integer(2)==Integer(24)
>>> solve([eq1,eq2,eq3,p==Integer(1)],p,q,x,y)
[[p == 1, q == 8, x == -4/3*sqrt(10) - 2/3, y == 1/6*sqrt(10) - 2/3], [p == 1, q == 8, x == 4/3*sqrt(10) - 2/3, y == -1/6*sqrt(10) - 2/3]]

Si queremos aproximaciones numéricas de las soluciones, podemos usar lo siguiente:

sage: solns = solve([eq1,eq2,eq3,p==1],p,q,x,y, solution_dict=True)
sage: [[s[p].n(30), s[q].n(30), s[x].n(30), s[y].n(30)] for s in solns]
[[1.0000000, 8.0000000, -4.8830369, -0.13962039],
 [1.0000000, 8.0000000, 3.5497035, -1.1937129]]
>>> from sage.all import *
>>> solns = solve([eq1,eq2,eq3,p==Integer(1)],p,q,x,y, solution_dict=True)
>>> [[s[p].n(Integer(30)), s[q].n(Integer(30)), s[x].n(Integer(30)), s[y].n(Integer(30))] for s in solns]
[[1.0000000, 8.0000000, -4.8830369, -0.13962039],
 [1.0000000, 8.0000000, 3.5497035, -1.1937129]]

(La función n imprime una aproximación numérica, y el argumento es el número de bits de precisión.)

Resolviendo Ecuaciones Numéricamente

A menudo, solve no podrá encontrar una solución exacta para la ecuación o ecuaciones especificadas. Cuando falla, puedes usar find_root para encontrar una solución numérica. Por ejemplo, solve no devuelve nada interesante para la siguiente ecuación:

sage: theta = var('theta')
sage: solve(cos(theta)==sin(theta), theta)
[sin(theta) == cos(theta)]
>>> from sage.all import *
>>> theta = var('theta')
>>> solve(cos(theta)==sin(theta), theta)
[sin(theta) == cos(theta)]

Por otro lado, podemos usar find_root para encontrar una solución a la ecuación de arriba en el rango \(0 < \theta < \pi/2\):

sage: phi = var('phi')
sage: find_root(cos(phi)==sin(phi),0,pi/2)
0.785398163397448...
>>> from sage.all import *
>>> phi = var('phi')
>>> find_root(cos(phi)==sin(phi),Integer(0),pi/Integer(2))
0.785398163397448...

Diferenciación, Integración, etc.

Sage sabe cómo diferenciar e integrar muchas funciones. Por ejemplo, para diferenciar \(\sin(u)\) con respecto a \(u\), haz lo siguiente:

sage: u = var('u')
sage: diff(sin(u), u)
cos(u)
>>> from sage.all import *
>>> u = var('u')
>>> diff(sin(u), u)
cos(u)

Para calcular la cuarta derivada de \(\sin(x^2)\):

sage: diff(sin(x^2), x, 4)
16*x^4*sin(x^2) - 48*x^2*cos(x^2) - 12*sin(x^2)
>>> from sage.all import *
>>> diff(sin(x**Integer(2)), x, Integer(4))
16*x^4*sin(x^2) - 48*x^2*cos(x^2) - 12*sin(x^2)

Para calcular las derivadas parciales de \(x^2+17y^2\) con respecto a x e y, respectivamente:

sage: x, y = var('x,y')
sage: f = x^2 + 17*y^2
sage: f.diff(x)
2*x
sage: f.diff(y)
34*y
>>> from sage.all import *
>>> x, y = var('x,y')
>>> f = x**Integer(2) + Integer(17)*y**Integer(2)
>>> f.diff(x)
2*x
>>> f.diff(y)
34*y

También podemos calcular integrales, tanto indefinidas como definidas. Para calcular \(\int x\sin(x^2)\, dx\) y \(\int_0^1 \frac{x}{x^2+1}\, dx\)

sage: integral(x*sin(x^2), x)
-1/2*cos(x^2)
sage: integral(x/(x^2+1), x, 0, 1)
1/2*log(2)
>>> from sage.all import *
>>> integral(x*sin(x**Integer(2)), x)
-1/2*cos(x^2)
>>> integral(x/(x**Integer(2)+Integer(1)), x, Integer(0), Integer(1))
1/2*log(2)

Para calcular la descomposición en fracciones simples de \(\frac{1}{x^2-1}\):

sage: f = 1/((1+x)*(x-1))
sage: f.partial_fraction(x)
-1/2/(x + 1) + 1/2/(x - 1)
>>> from sage.all import *
>>> f = Integer(1)/((Integer(1)+x)*(x-Integer(1)))
>>> f.partial_fraction(x)
-1/2/(x + 1) + 1/2/(x - 1)

Resolviendo Ecuaciones Diferenciales

Puedes usar a Sage para investigar ecuaciones diferenciales ordinarias. Para resolver la ecuación \(x'+x-1=0\):

sage: t = var('t')    # defina una variable t
sage: x = function('x')(t)   # defina x como una función de esa variable
sage: DE = diff(x, t) + x - 1
sage: desolve(DE, [x,t])
(_C + e^t)*e^(-t)
>>> from sage.all import *
>>> t = var('t')    # defina una variable t
>>> x = function('x')(t)   # defina x como una función de esa variable
>>> DE = diff(x, t) + x - Integer(1)
>>> desolve(DE, [x,t])
(_C + e^t)*e^(-t)

Esto utiliza el interfaz a Maxima de Sage [Max], por lo que el resultado puede diferir de otros resultados de Sage. En este caso, la salida nos dice que la solución general a la ecuación diferencial es \(x(t) = e^{-t}(e^{t}+c)\).

También puedes calcular transformadas de Laplace; la transformada de Laplace de \(t^2e^t -\sin(t)\) se calcula como sigue:

sage: s = var("s")
sage: t = var("t")
sage: f = t^2*exp(t) - sin(t)
sage: f.laplace(t,s)
-1/(s^2 + 1) + 2/(s - 1)^3
>>> from sage.all import *
>>> s = var("s")
>>> t = var("t")
>>> f = t**Integer(2)*exp(t) - sin(t)
>>> f.laplace(t,s)
-1/(s^2 + 1) + 2/(s - 1)^3

Veamos un ejemplo más complicado. El desplazamiento desde el punto de equilibrio de dos resortes acoplados, sujetos a una pared a la izquierda

|------\/\/\/\/\---|masa1|----\/\/\/\/\/----|masa2|
         resorte1               resorte2

está modelado por el sistema de ecuaciones diferenciales de segundo órden

\[ \begin{align}\begin{aligned}m_1 x_1'' + (k_1+k_2) x_1 - k_2 x_2 = 0\\m_2 x_2''+ k_2 (x_2-x_1) = 0,\end{aligned}\end{align} \]

donde \(m_{i}\) es la masa del objeto i, \(x_{i}\) es el desplazamiento desde el equilibrio de la masa i, y \(k_{i}\) es la constante de elasticidad del resorte i.

Ejemplo: Utiliza Sage para resolver el problema de arriba con \(m_{1}=2\), \(m_{2}=1\), \(k_{1}=4\), \(k_{2}=2\), \(x_{1}(0)=3\), \(x_{1}'(0)=0\), \(x_{2}(0)=3\), \(x_{2}'(0)=0\).

Solución: Toma la transformada de Laplace de la primera ecuación (con la notación \(x=x_{1}\), \(y=x_{2}\)):

sage: de1 = maxima("2*diff(x(t),t, 2) + 6*x(t) - 2*y(t)")
sage: lde1 = de1.laplace("t","s"); lde1.sage()
2*s^2*laplace(x(t), t, s) - 2*s*x(0) + 6*laplace(x(t), t, s) - 2*laplace(y(t), t, s) - 2*D[0](x)(0)
>>> from sage.all import *
>>> de1 = maxima("2*diff(x(t),t, 2) + 6*x(t) - 2*y(t)")
>>> lde1 = de1.laplace("t","s"); lde1.sage()
2*s^2*laplace(x(t), t, s) - 2*s*x(0) + 6*laplace(x(t), t, s) - 2*laplace(y(t), t, s) - 2*D[0](x)(0)

El resultado puede ser difícil de leer, pero significa que

\[-2x'(0) + 2s^2*X(s) - 2sx(0) - 2Y(s) + 6X(s) = 0\]

(donde la transformada de Laplace de una función en letra minúscula como \(x(t)\) es la función en letra mayúscula \(X(s)\)). Toma la transformada de Laplace de la segunda ecuación:

sage: t,s = SR.var('t,s')
sage: x = function('x')
sage: y = function('y')
sage: f = 2*x(t).diff(t,2) + 6*x(t) - 2*y(t)
sage: f.laplace(t,s)
2*s^2*laplace(x(t), t, s) - 2*s*x(0) + 6*laplace(x(t), t, s) - 2*laplace(y(t), t, s) - 2*D[0](x)(0)
>>> from sage.all import *
>>> t,s = SR.var('t,s')
>>> x = function('x')
>>> y = function('y')
>>> f = Integer(2)*x(t).diff(t,Integer(2)) + Integer(6)*x(t) - Integer(2)*y(t)
>>> f.laplace(t,s)
2*s^2*laplace(x(t), t, s) - 2*s*x(0) + 6*laplace(x(t), t, s) - 2*laplace(y(t), t, s) - 2*D[0](x)(0)

Esto dice

\[-Y'(0) + s^2Y(s) + 2Y(s) - 2X(s) - sy(0) = 0.\]

Introduce las condiciones iniciales para \(x(0)\), \(x'(0)\), \(y(0)\) y \(y'(0)\) y resuelve las dos ecuaciones resultantes:

sage: var('s X Y')
(s, X, Y)
sage: eqns = [(2*s^2+6)*X-2*Y == 6*s, -2*X +(s^2+2)*Y == 3*s]
sage: solve(eqns, X,Y)
[[X == 3*(s^3 + 3*s)/(s^4 + 5*s^2 + 4),
  Y == 3*(s^3 + 5*s)/(s^4 + 5*s^2 + 4)]]
>>> from sage.all import *
>>> var('s X Y')
(s, X, Y)
>>> eqns = [(Integer(2)*s**Integer(2)+Integer(6))*X-Integer(2)*Y == Integer(6)*s, -Integer(2)*X +(s**Integer(2)+Integer(2))*Y == Integer(3)*s]
>>> solve(eqns, X,Y)
[[X == 3*(s^3 + 3*s)/(s^4 + 5*s^2 + 4),
  Y == 3*(s^3 + 5*s)/(s^4 + 5*s^2 + 4)]]

Ahora toma la transformada inversa de Laplace para obtener la respuesta:

sage: var('s t')
(s, t)
sage: inverse_laplace((3*s^3 + 9*s)/(s^4 + 5*s^2 + 4),s,t)
cos(2*t) + 2*cos(t)
sage: inverse_laplace((3*s^3 + 15*s)/(s^4 + 5*s^2 + 4),s,t)
-cos(2*t) + 4*cos(t)
>>> from sage.all import *
>>> var('s t')
(s, t)
>>> inverse_laplace((Integer(3)*s**Integer(3) + Integer(9)*s)/(s**Integer(4) + Integer(5)*s**Integer(2) + Integer(4)),s,t)
cos(2*t) + 2*cos(t)
>>> inverse_laplace((Integer(3)*s**Integer(3) + Integer(15)*s)/(s**Integer(4) + Integer(5)*s**Integer(2) + Integer(4)),s,t)
-cos(2*t) + 4*cos(t)

Por tanto, la solución es

\[x_1(t) = \cos(2t) + 2\cos(t), \quad x_2(t) = 4\cos(t) - \cos(2t).\]

La solución puede dibujarse paramétricamente usando

sage: t = var('t')
sage: P = parametric_plot((cos(2*t) + 2*cos(t), 4*cos(t) - cos(2*t) ),\
....: (0, 2*pi), rgbcolor=hue(0.9))
sage: show(P)
>>> from sage.all import *
>>> t = var('t')
>>> P = parametric_plot((cos(Integer(2)*t) + Integer(2)*cos(t), Integer(4)*cos(t) - cos(Integer(2)*t) ),(Integer(0), Integer(2)*pi), rgbcolor=hue(RealNumber('0.9')))
>>> show(P)

Los componentes individuales pueden dibujarse usando

sage: t = var('t')
sage: p1 = plot(cos(2*t) + 2*cos(t), 0, 2*pi, rgbcolor=hue(0.3))
sage: p2 = plot(4*cos(t) - cos(2*t), 0, 2*pi, rgbcolor=hue(0.6))
sage: show(p1 + p2)
>>> from sage.all import *
>>> t = var('t')
>>> p1 = plot(cos(Integer(2)*t) + Integer(2)*cos(t), Integer(0), Integer(2)*pi, rgbcolor=hue(RealNumber('0.3')))
>>> p2 = plot(Integer(4)*cos(t) - cos(Integer(2)*t), Integer(0), Integer(2)*pi, rgbcolor=hue(RealNumber('0.6')))
>>> show(p1 + p2)

REFERENCIAS: Nagle, Saff, Snider, Fundamentos De Ecuaciones Diferenciales, 6a ed, Addison-Wesley, 2004. (véase § 5.5).

Método De Euler Para Sistemas De Ecuaciones Diferenciales

En el siguiente ejemplo, ilustraremos el método de Euler para EDOs de primer y segundo órden. Primero, recordemos la idea básica para ecuaciones de primer órden. Dado un problema con valor inicial de la forma

\[y'=f(x,y) y(a)=c\]

queremos encontrar el valor aproximado de la solución en \(x=b\) con \(b>a\).

Recuerda de la definición de derivada que

\[y'(x) \approx \frac{y(x+h)-y(x)}{h},\]

donde \(h>0\) está dado y es pequeño. Esto, junto con la ED, dan \(f(x,y(x))\approx \frac{y(x+h)-y(x)}{h}\). Ahora resuelve para \(y(x+h)\):

\[y(x+h) \approx y(x) + h*f(x,y(x)).\]

Si llamamos a \(h f(x,y(x))\) el «término de corrección» (a falta de algo mejor), llamamos a \(y(x)\) «el valor viejo de y», y llamamos a \(y(x+h)\) el «nuevo valor de y», entonces, esta aproximación puede re-expresarse como

\[y_{nuevo} \approx y_{viejo} + h*f(x,y_{viejo}).\]

Si descomponemos el intervalo desde a a b en n pasos, de modo que \(h=\frac{b-a}{n}\), podemos guardar la información dada por este método en una tabla.

\(x\)

\(y\)

\(hf(x,y)\)

\(a\)

\(c\)

\(hf(a,c)\)

\(a+h\)

\(c+hf(a,c)\)

\(a+2h\)

\(b=a+nh\)

???

La meta es llenar todos los espacios de la tabla, una fila cada la vez, hasta que lleguemos a la casilla ???, que será la aproximación del método de Euler para \(y(b)\).

La idea para los sistemas de EDOs es similar.

Ejemplo: Aproxima numéricamente \(z(t)\) en \(t=1\) usando 4 pasos del método de Euler, donde \(z''+tz'+z=0\), \(z(0)=1\), \(z'(0)=0\).

Debemos reducir la EDO de segundo órden a un sistema de dos EDs de primer órden (usando \(x=z\), \(y=z'\)) y aplicar el método de Euler:

sage: t,x,y = PolynomialRing(RealField(10),3,"txy").gens()
sage: f = y; g = -x - y * t
sage: eulers_method_2x2(f,g, 0, 1, 0, 1/4, 1)
      t                x            h*f(t,x,y)                y       h*g(t,x,y)
      0                1                  0.00                0           -0.25
    1/4              1.0                -0.062            -0.25           -0.23
    1/2             0.94                 -0.12            -0.48           -0.17
    3/4             0.82                 -0.16            -0.66          -0.081
      1             0.65                 -0.18            -0.74           0.022
>>> from sage.all import *
>>> t,x,y = PolynomialRing(RealField(Integer(10)),Integer(3),"txy").gens()
>>> f = y; g = -x - y * t
>>> eulers_method_2x2(f,g, Integer(0), Integer(1), Integer(0), Integer(1)/Integer(4), Integer(1))
      t                x            h*f(t,x,y)                y       h*g(t,x,y)
      0                1                  0.00                0           -0.25
    1/4              1.0                -0.062            -0.25           -0.23
    1/2             0.94                 -0.12            -0.48           -0.17
    3/4             0.82                 -0.16            -0.66          -0.081
      1             0.65                 -0.18            -0.74           0.022

Por tanto, \(z(1)\approx 0.75\).

También podemos dibujar los puntos \((x,y)\) para obtener una representación aproximada de la curva. La función que hace esto es eulers_method_2x2_plot. Para poder usarla, necesitamos definir las funciones f y g que toman un argumento con tres coordenadas: (t, x,*y*).

sage: f = lambda z: z[2]        # f(t,x,y) = y
sage: g = lambda z: -sin(z[1])  # g(t,x,y) = -sin(x)
sage: P = eulers_method_2x2_plot(f,g, 0.0, 0.75, 0.0, 0.1, 1.0)
>>> from sage.all import *
>>> f = lambda z: z[Integer(2)]        # f(t,x,y) = y
>>> g = lambda z: -sin(z[Integer(1)])  # g(t,x,y) = -sin(x)
>>> P = eulers_method_2x2_plot(f,g, RealNumber('0.0'), RealNumber('0.75'), RealNumber('0.0'), RealNumber('0.1'), RealNumber('1.0'))

A estas alturas, P está guardando dos gráficas: P[0], el gráfico de x vs. t, y P[1], el gráfico de y vs. t. Podemos mostrar ámbas como sigue:

sage: show(P[0] + P[1])
>>> from sage.all import *
>>> show(P[Integer(0)] + P[Integer(1)])

Funciones Especiales

Se han implementado varios polinomios ortogonales y funciones especiales, utilizando tanto PARI [GAP] como Maxima [Max]. Estas funciones están documentadas en las secciones apropiadas («Polinomios Ortogonales» y «Funciones Especiales», respectivamente) del manual de referencia de Sage.

sage: x = polygen(QQ, 'x')
sage: chebyshev_U(2,x)
4*x^2 - 1
sage: bessel_I(1,1).n(250)
0.56515910399248502720769602760986330732889962162109200948029448947925564096
sage: bessel_I(1,1).n()
0.565159103992485
sage: bessel_I(2,1.1).n()  # los últimos digitos son al azar
0.16708949925104...
>>> from sage.all import *
>>> x = polygen(QQ, 'x')
>>> chebyshev_U(Integer(2),x)
4*x^2 - 1
>>> bessel_I(Integer(1),Integer(1)).n(Integer(250))
0.56515910399248502720769602760986330732889962162109200948029448947925564096
>>> bessel_I(Integer(1),Integer(1)).n()
0.565159103992485
>>> bessel_I(Integer(2),RealNumber('1.1')).n()  # los últimos digitos son al azar
0.16708949925104...

Hasta este punto, Sage únicamente ha encapsulado estas funciones para uso numérico. Para uso simbólico, por favor utiliza directamente la interfaz a Maxima, como en el siguiente ejemplo:

sage: maxima.eval("f:bessel_y(v, w)")
'bessel_y(v,w)'
sage: maxima.eval("diff(f,w)")
'(bessel_y(v-1,w)-bessel_y(v+1,w))/2'
>>> from sage.all import *
>>> maxima.eval("f:bessel_y(v, w)")
'bessel_y(v,w)'
>>> maxima.eval("diff(f,w)")
'(bessel_y(v-1,w)-bessel_y(v+1,w))/2'
[GAP]

El Grupo GAP, GAP - Grupos, Algorítmos y Programación, https://www.gap-system.org

[Max] (1,2)

Maxima, http://maxima.sf.net/