Álgebra Elementar e Cálculo

O Sage pode realizar diversos cálculos em álgebra elementar e cálculo diferencial e integral: por exemplo, encontrar soluções de equações, diferenciar, integrar, e calcular a transformada de Laplace. Veja a documentação em Sage Constructions para mais exemplos.

Resolvendo equações

Resolvendo equações exatamente

A função solve resolve equações. Para usá-la, primeiro especifique algumas variáveis; então os argumentos de solve são uma equação (ou um sistema de equações), juntamente com as variáveis para as quais 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]

Você pode resolver equações para uma variável em termos das outras:

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

Você pode resolver para diversas variáveis:

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]]

O seguinte exemplo, que mostra como usar o Sage para resolver um sistema de equações não-lineares, foi sugerido por Jason Grout: primeiro, resolvemos o sistemas simbolicamente:

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]]

Para obter soluções numéricas aproximadas, podemos usar:

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]]

(A função n imprime uma aproximação numérica, e o argumento é o número de bits de precisão.)

Resolvendo Equações Numericamente

Frequentemente, solve não será capaz de encontrar uma solução exata para uma equação ou sistema de equações. Nesse caso, você pode usar find_root para encontrar uma solução numérica. Por exemplo, solve não encontra uma solução para a equação abaixo:

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 outro lado, podemos usar find_root para encontrar uma solução para a equação acima no intervalo \(0 < \phi < \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...

Diferenciação, Integração, etc.

O Sage é capaz de diferenciar e integrar diversas funções. Por exemplo, para diferenciar \(\sin(u)\) com respeito a \(u\), faça o seguinte:

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 a quarta 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 as derivadas parciais de \(x^2+17y^2\) com respeito 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

Passamos agora para integrais, tanto indefinidas como definidas. Para calcular \(\int x\sin(x^2)\, dx\) e \(\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 a decomposição em frações parciais 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)

Resolvendo Equações Diferenciais

Você pode usar o Sage para investigar equações diferenciais ordinárias. Para resolver a equação \(x'+x-1=0\):

sage: t = var('t')    # define a variable t
sage: x = function('x')(t)   # define x to be a function of that 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')    # define a variable t
>>> x = function('x')(t)   # define x to be a function of that variable
>>> DE = diff(x, t) + x - Integer(1)
>>> desolve(DE, [x,t])
(_C + e^t)*e^(-t)

Esse método usa a interface do Sage para o Maxima [Max]. Logo, o formato dos resultados é um pouco diferente de outros cálculos realizados no Sage. Nesse caso, o resultado diz que a solução geral da equação diferencial é \(x(t) = e^{-t}(e^{t}+c)\).

Você pode calcular a transformada de Laplace também; a transformada de Laplace de \(t^2e^t -\sin(t)\) é calculada da seguinte forma:

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

A seguir, um exemplo mais complicado. O deslocamento, com respeito à posição de equilíbrio, de duas massas presas a uma parede através de molas, conforme a figura abaixo,

|------\/\/\/\/\---|massa1|----\/\/\/\/\/----|massa2|
         mola1                    mola2

é modelado pelo sistema de equações diferenciais de segunda ordem

\[ \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} \]

onde, para \(i=1,2\), \(m_{i}\) é a massa do objeto i, \(x_{i}\) é o deslocamento com respeito à posição de equilíbrio da massa i, e \(k_{i}\) é a constante de mola para a mola i.

Exemplo: Use o Sage para resolver o problema acima com \(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\).

Solução: Primeiramente, calcule a transformada de Laplace da primeira equação (usando a notação \(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)

O resultado é um pouco difícil de ler, mas diz que

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

(onde a transformada de Laplace de uma função em letra minúscula \(x(t)\) é a função em letra maiúscula \(X(s)\)). Agora, calcule a transformada de Laplace da segunda equação:

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)

O resultado significa que

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

Em seguida, substitua a condição inicial para \(x(0)\), \(x'(0)\), \(y(0)\), e \(y'(0)\), e resolva as equações 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)]]

Agora calcule a transformada de Laplace inversa para obter a resposta:

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)

Portanto, a solução é

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

Ela pode ser representada em um gráfico parametricamente usando os comandos

sage: t = var('t')
sage: P = parametric_plot((cos(2*t) + 2*cos(t), 4*cos(t) - cos(2*t) ),
....: (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) ),
... (t, Integer(0), Integer(2)*pi), rgbcolor=hue(RealNumber('0.9')))
>>> show(P)

As componentes individuais podem ser representadas em gráfico usando

sage: t = var('t')
sage: p1 = plot(cos(2*t) + 2*cos(t), (t,0, 2*pi), rgbcolor=hue(0.3))
sage: p2 = plot(4*cos(t) - cos(2*t), (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), (t,Integer(0), Integer(2)*pi), rgbcolor=hue(RealNumber('0.3')))
>>> p2 = plot(Integer(4)*cos(t) - cos(Integer(2)*t), (t,Integer(0), Integer(2)*pi), rgbcolor=hue(RealNumber('0.6')))
>>> show(p1 + p2)

Leia mais sobre gráficos em Gráficos. Veja a seção 5.5 de [NagleEtAl2004] (em inglês) para mais informações sobre equações diferenciais.

Método de Euler para Sistemas de Equações Diferenciais

No próximo exemplo, vamos ilustrar o método de Euler para EDOs de primeira e segunda ordem. Primeiro, relembramos a ideia básica para equações de primeira ordem. Dado um problema de valor inicial da forma

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

queremos encontrar o valor aproximado da solução em \(x=b\) com \(b>a\).

Da definição de derivada segue que

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

onde \(h>0\) é um número pequeno. Isso, juntamente com a equação diferencial, implica que \(f(x,y(x))\approx \frac{y(x+h)-y(x)}{h}\). Agora resolvemos para \(y(x+h)\):

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

Se chamarmos \(h f(x,y(x))\) de “termo de correção”, \(y(x)\) de “valor antigo de y”, e \(y(x+h)\) de “novo valor de y”, então essa aproximação pode ser reescrita como

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

Se dividirmos o intervalo de a até b em n partes, de modo que \(h=\frac{b-a}{n}\), então podemos construir a seguinte tabela.

\(x\)

\(y\)

\(hf(x,y)\)

\(a\)

\(c\)

\(hf(a,c)\)

\(a+h\)

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

\(a+2h\)

\(b=a+nh\)

???

O objetivo é completar os espaços em branco na tabela, em uma linha por vez, até atingirmos ???, que é a aproximação para \(y(b)\) usando o método de Euler.

A ideia para sistemas de EDOs é semelhante.

Exemplo: Aproxime numericamente \(z(t)\) em \(t=1\) usando 4 passos do método de Euler, onde \(z''+tz'+z=0\), \(z(0)=1\), \(z'(0)=0\).

Devemos reduzir a EDO de segunda ordem a um sistema de duas EDOs de primeira ordem (usando \(x=z\), \(y=z'\)) e aplicar o 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

Portanto, \(z(1)\approx 0.65\).

Podemos também representar em um gráfico os pontos \((x,y)\) para obter uma figura da solução aproximada. A função eulers_method_2x2_plot fará isso; para usá-la, precisamos definir funções f e g que recebam um argumento com três 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 esta altura, P armazena dois gráficos: P[0], o gráfico de x versus t, e P[1], o gráfico de y versus t. Podemos visualizar os dois gráficos da seguinte forma:

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

(Para mais sobre gráficos, veja Gráficos.)

Funções Especiais

Diversos polinômios ortogonais e funções especiais estão implementadas, usando tanto o PARI [GP] como o Maxima [Max]. Isso está documentado nas seções apropriadas (“Orthogonal polynomials” and “Special functions”, respectivamente) do manual de referência do Sage (em inglês).

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.56515910399248...
sage: bessel_I(2,1.1).n()  # last few digits are random
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.56515910399248...
>>> bessel_I(Integer(2),RealNumber('1.1')).n()  # last few digits are random
0.16708949925104...

No momento, essas funções estão disponíveis na interface do Sage apenas para uso numérico. Para uso simbólico, use a interface do Maxima diretamente, como no seguinte exemplo:

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'