Algèbre de base et calcul infinitésimal

Sage peut accomplir divers calculs d’algèbre et d’analyse de base : par exemple, trouver les solutions d’équations, dériver, intégrer, calculer des transformées de Laplace. Voir la documentation Sage Constructions pour plus d’exemples.

Résolution d’équations

La fonction solve résout des équations. Pour l’utiliser, il convient de spécifier d’abord les variables. Les arguments de solve sont alors une équation (ou un système d’équations) suivie des variables à résoudre.

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]

x = var('x')
solve(x^2 + 3*x + 2, x)

On peut résoudre une équation en une variable en fonction des autres :

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

x, b, c = var('x b c')
solve([x^2 + b*x + c == 0],x)

On peut également résoudre un système à plusieurs 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]]

x, y = var('x, y')
solve([x+y==6, x-y==4], x, y)

L’exemple suivant, qui utilise Sage pour la résolution d’un système d’équations non-linéaires, a été proposé par Jason Grout. D’abord, on résout le système de façon symbolique :

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

var('x y p q')
eq1 = p+q==9
eq2 = q*y+p*x==-6
eq3 = q*y^2+p*x^2==24
solve([eq1,eq2,eq3,p==1],p,q,x,y)

Pour une résolution numérique, on peut utiliser à la place :

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

solns = solve([eq1,eq2,eq3,p==1],p,q,x,y, solution_dict=True)
[[s[p].n(30), s[q].n(30), s[x].n(30), s[y].n(30)] for s in solns]

(La fonction n affiche une approximation numérique ; son argument indique le nombre de bits de précision.)

Dérivation, intégration, etc.

Sage est capable de dériver et d’intégrer de nombreuses fonctions. Par exemple, pour dériver \(\sin(u)\) par rapport à \(u\), on procède comme suit :

sage: u = var('u')
sage: diff(sin(u), u)
cos(u)

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

u = var('u')
diff(sin(u), u)

Pour calculer la dérivée quatrième 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)

diff(sin(x^2), x, 4)

Pour calculer la dérivée partielle de \(x^2+17y^2\) par rapport à \(x\) et \(y\) respectivement :

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

x, y = var('x,y')
f = x^2 + 17*y^2
f.diff(x)
f.diff(y)

Passons aux primitives et intégrales. Pour calculer \(\int x\sin(x^2)\, dx\) et \(\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)

integral(x*sin(x^2), x)
integral(x/(x^2+1), x, 0, 1)

Pour calculer la décomposition en éléments 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)

f = 1/((1+x)*(x-1))
f.partial_fraction(x)

Résolution des équations différentielles

On peut utiliser Sage pour étudier les équations différentielles ordinaires. Pour résoudre l’équation \(x'+x-1=0\) :

sage: t = var('t')    # on définit une variable t
sage: function('x')(t)   # on déclare x fonction de cette variable
x(t)
sage: DE = lambda y: diff(y,t) + y - 1
sage: desolve(DE(x(t)), [x(t),t])
(_C + e^t)*e^(-t)

>>> from sage.all import *
>>> t = var('t')    # on définit une variable t
>>> function('x')(t)   # on déclare x fonction de cette variable
x(t)
>>> DE = lambda y: diff(y,t) + y - Integer(1)
>>> desolve(DE(x(t)), [x(t),t])
(_C + e^t)*e^(-t)

t = var('t')    # on définit une variable t
function('x')(t)   # on déclare x fonction de cette variable
DE = lambda y: diff(y,t) + y - 1
desolve(DE(x(t)), [x(t),t])

Ceci utilise l’interface de Sage vers Maxima [Max], aussi il se peut que la sortie diffère un peu des sorties habituelles de Sage. Dans notre cas, le résultat indique que la solution générale à l’équation différentielle est \(x(t) = e^{-t}(e^{t}+C)\).

Il est aussi possible de calculer des transformées de Laplace. La transformée de Laplace de \(t^2e^t -\sin(t)\) s’obtient comme suit :

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

s = var("s")
t = var("t")
f = t^2*exp(t) - sin(t)
f.laplace(t,s)

Voici un exemple plus élaboré. L’élongation à partir du point d’équilibre de ressorts couplés attachés à gauche à un mur

|------\/\/\/\/\---|masse1|----\/\/\/\/\/----|masse2|
        ressort1                ressort2

est modélisée par le système d’équations différentielles d’ordre 2

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

où \(m_{i}\) est la masse de l’objet i, \(x_{i}\) est l’élongation à partir du point d’équilibre de la masse i, et \(k_{i}\) est la constante de raideur du ressort i.

Exemple : Utiliser Sage pour résoudre le problème ci-dessus avec \(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\).

Solution : Considérons la transformée de Laplace de la première équation (avec les notations \(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)

de1 = maxima("2*diff(x(t),t, 2) + 6*x(t) - 2*y(t)")
lde1 = de1.laplace("t","s"); lde1.sage()

La réponse n’est pas très lisible, mais elle signifie que

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

(où la transformée de Laplace d’une fonction notée par une lettre minuscule telle que \(x(t)\) est désignée par la majuscule correspondante \(X(s)\)). Considérons la transformée de Laplace de la seconde équation :

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)

t,s = SR.var('t,s')
x = function('x')
y = function('y')
f = 2*x(t).diff(t,2) + 6*x(t) - 2*y(t)
f.laplace(t,s)

Ceci signifie

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

Injectons les conditions initiales pour \(x(0)\), \(x'(0)\), \(y(0)\) et \(y'(0)\) et résolvons les deux équations qui en résultent :

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

var('s X Y')
eqns = [(2*s^2+6)*X-2*Y == 6*s, -2*X +(s^2+2)*Y == 3*s]
solve(eqns, X,Y)

À présent, prenons la transformée de Laplace inverse pour obtenir la réponse :

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)

var('s t')
inverse_laplace((3*s^3 + 9*s)/(s^4 + 5*s^2 + 4),s,t)
inverse_laplace((3*s^3 + 15*s)/(s^4 + 5*s^2 + 4),s,t)

Par conséquent, la solution est

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

On peut en tracer le graphe paramétrique en utilisant

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)

t = var('t')
P = parametric_plot((cos(2*t) + 2*cos(t), 4*cos(t) - cos(2*t) ),
    (t, 0, 2*pi), rgbcolor=hue(0.9))
show(P)

Les coordonnées individuelles peuvent être tracées en utilisant

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)

t = var('t')
p1 = plot(cos(2*t) + 2*cos(t), (t, 0, 2*pi), rgbcolor=hue(0.3))
p2 = plot(4*cos(t) - cos(2*t), (t, 0, 2*pi), rgbcolor=hue(0.6))
show(p1 + p2)

Les fonctions de tracé de graphes sont décrites dans la section Graphiques de ce tutoriel. On pourra aussi consulter [NagleEtAl2004], §5.5 pour plus d’informations sur les équations différentielles.

Méthode d’Euler pour les systèmes d’équations différentielles

Dans l’exemple suivant, nous illustrons la méthode d’Euler pour des équations différentielles ordinaires d’ordre un et deux. Rappelons d’abord le principe de la méthode pour les équations du premier ordre. Etant donné un problème donné avec une valeur initiale sous la forme

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

nous cherchons une valeur approchée de la solution au point \(x=b\) avec \(b>a\).

Rappelons que par définition de la dérivée

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

où \(h>0\) est fixé et petit. Ceci, combiné à l’équation différentielle, donne \(f(x,y(x))\approx \frac{y(x+h)-y(x)}{h}\). Aussi \(y(x+h)\) s’écrit:

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

Si nous notons \(h\cdot f(x,y(x))\) le « terme de correction » (faute d’un terme plus approprié), et si nous appelons \(y(x)\) « l’ancienne valeur de \(y\) » et \(y(x+h)\) la « nouvelle valeur de \(y\) », cette approximation se réécrit

\[y_{nouveau} \approx y_{ancien} + h\cdot f(x,y_{ancien}).\]

Divisions l’intervalle entre \(a\) et \(b\) en \(n\) pas, si bien que \(h=\frac{b-a}{n}\). Nous pouvons alors remplir un tableau avec les informations utilisées dans la méthode.

\(x\)

\(y\)

\(h\cdot f(x,y)\)

\(a\)

\(c\)

\(h\cdot f(a,c)\)

\(a+h\)

\(c+h\cdot f(a,c)\)

…

\(a+2h\)

…

…

\(b=a+nh\)

???

…

Le but est est de remplir tous les trous du tableau, ligne après ligne, jusqu’à atteindre le coefficient « ??? », qui est l’approximation de \(y(b)\) au sens de la méthode d’Euler.

L’idée est la même pour les systèmes d’équations différentielles.

Exemple: Rechercher une approximation numérique de \(z(t)\) en \(t=1\) en utilisant 4 étapes de la méthode d’Euler, où \(z''+tz'+z=0\), \(z(0)=1\), \(z'(0)=0\).

Il nous faut réduire l’équation différentielle d’ordre 2 à un système de deux équations différentielles d’ordre 1 (en posant \(x=z\), \(y=z'\)) et appliquer la méthode d’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

t,x,y = PolynomialRing(RealField(10),3,"txy").gens()
f = y; g = -x - y * t
eulers_method_2x2(f,g, 0, 1, 0, 1/4, 1)

On en déduit \(z(1)\approx 0.65\).

On peut également tracer le graphe des points \((x,y)\) pour obtenir une image approchée de la courbe. La fonction eulers_method_2x2_plot réalise cela ; pour l’utiliser, il faut définir les fonctions \(f\) et \(g\) qui prennent un argument à trois coordonnées : (\(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'))

f = lambda z: z[2]        # f(t,x,y) = y
g = lambda z: -sin(z[1])  # g(t,x,y) = -sin(x)
P = eulers_method_2x2_plot(f,g, 0.0, 0.75, 0.0, 0.1, 1.0)

Arrivé à ce point, P conserve en mémoire deux graphiques : P[0], le graphe de \(x\) en fonction de \(t\), et P[1], le graphique de \(y\) par rapport à \(t\). On peut tracer les deux graphiques simultanément par :

sage: show(P[0] + P[1])

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

show(P[0] + P[1])

(Pour plus d’information sur le tracé de graphiques, voir Graphiques.)

Fonctions spéciales

Plusieurs familles de polynômes orthogonaux et fonctions spéciales sont implémentées via PARI [GAP] et Maxima [Max]. Ces fonctions sont documentées dans les sections correspondantes (Orthogonal polynomials et Special functions, respectively) du manuel de référence de Sage (Sage reference manual).

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

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

x = polygen(QQ, 'x')
chebyshev_U(2,x)
bessel_I(1,1).n(250)
bessel_I(1,1).n()
bessel_I(2,1.1).n()

Pour l’instant, ces fonctions n’ont été adaptées à Sage que pour une utilisation numérique. Pour faire du calcul formel, il faut utiliser l’interface Maxima directement, comme le présente l’exemple suivant :

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'

maxima.eval("f:bessel_y(v, w)")
maxima.eval("diff(f,w)")