Algebra di base e Analisi#

Sage sa svolgere diversi calcoli legati all’algebra di base ed all’analisi: per esempio, risoluzione di equazioni, calcolo differenziale ed integrale e trasformate di Laplace. Si veda la documentazione per le «Costruzioni di Sage» per ulteriori esempi.

Risoluzione di equazioni#

La funzione solve risolve le equazioni. Per usarla, bisogna anzitutto specificare alcune variabili; pertanto gli argomenti di solve sono un’equazione (od un sistema di equazioni), insieme con le variabili rispetto alle quali risolvere:

sage: x = var('x')
sage: solve(x^2 + 3*x + 2, x)
[x == -2, x == -1]

Si possono risolvere le equazioni rispetto ad una variabile in funzione delle altre:

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

Si può anche risolvere rispetto a diverse variabili:

sage: x, y = var('x, y')
sage: solve([x+y==6, x-y==4], x, y)
[[x == 5, y == 1]]

Il seguente esempio dell’uso di Sage per risolvere un sistema di equazioni non lineari è stato fornito da Jason Grout: per prima cosa, si risolve il sistema 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]]

Per una soluzione numerica, si può invece usare:

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

(La funzione n scrive un’approssimazione numerica, e l’argomento è il numero di bit di precisione.)

Differenziazione, Integrazione, etc.#

Sage è in grado di differenziae ed integrare molte funzioni. Per esempio, per differenziare $\sin (u)$ rispetto a $u$ , si procede come nelle righe seguenti:

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

Per calcolare la derivata quarta di $\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)

Per calcolare le derivate parziali di $x^{2} + 17 y^{2}$ rispetto a x e y, rispettivamente:

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

Passiamo agli integrali, sia indefiniti che definiti. Per calcolare $\int x \sin (x^{2}) d x$ e $\int_{0}^{1} \frac{x}{x^{2} + 1} d x$

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)

Per calcolare la decomposizione in frazioni parziali di $\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)

Risoluzione di Equazioni Differenziali#

Si può usare Sage per studiare le equazioni differenziali ordinarie. Per risolvere l’equazione $x^{'} + x - 1 = 0$ :

sage: t = var('t')    # definisce una variabile t
sage: x = function('x')(t)   # definisce x come funzione di quella variabile
sage: DE = diff(x,t) + x - 1
sage: desolve(DE, [x,t])
(_C + e^t)*e^(-t)

Questo metodo utilizza l’interfaccia di Sage per Maxima [Max], e così il suo output può essere leggermente diverso dagli altri output di Sage. In questo caso, risulta che la soluzione generale dell’equazione differenziale è $x (t) = e^{- t} (e^{t} + c)$ .

Si può anche calcolare la trasformata di Laplace; la trasformata di Laplace di $t^{2} e^{t} - \sin (t)$ è calcolata come segue:

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

Il successivo è un esempio più articolato. Lo scostamento dall’equilibrio (rispettivamente) per due molle accoppiate fissate ad un muro a sinistra

|------\/\/\/\/\---|massa1|----\/\/\/\/\/----|massa2|
         molla1                  molla2

è modellizzato dal sistema di equazioni differenziali del secondo ordine

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,

dove $m_{i}$ è la massa dell’oggetto i, $x_{i}$ è lo scostamento dall’equilibrio della massa i, e $k_{i}$ è la costante elastica della molla i.

Esempio: Usare Sage per risolvere il problema precedente 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$ .

Soluzione: Calcolare la trasformata di Laplace della prima equazione (con la notazione $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
2*((-%at('diff(x(t),t,1),t = 0))+s^2*'laplace(x(t),t,s)-x(0)*s) -2*'laplace(y(t),t,s)+6*'laplace(x(t),t,s)

Questo è di difficile lettura, ma dice che

- 2 x^{'} (0) + 2 s^{2} * X (s) - 2 s x (0) - 2 Y (s) + 6 X (s) = 0

(dove la trasformata di Laplace di una funzione in minuscolo come $x (t)$ è la funzione in maiuscolo $X (s)$ ). Calcolare la trasformata di Laplace della seconda equazione:

sage: de2 = maxima("diff(y(t),t, 2) + 2*y(t) - 2*x(t)")
sage: lde2 = de2.laplace("t","s"); lde2
(-%at('diff(y(t),t,1),t = 0))+s^2*'laplace(y(t),t,s) +2*'laplace(y(t),t,s)-2*'laplace(x(t),t,s) -y(0)*s

che significa

- Y^{'} (0) + s^{2} Y (s) + 2 Y (s) - 2 X (s) - s y (0) = 0.

Imporre le condizioni iniziali per $x (0)$ , $x^{'} (0)$ , $y (0)$ , e $y^{'} (0)$ , e risolvere le due equazioni risultanti:

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

Ora si calcola la trasformata inversa di Laplace per ottenere la risposta:

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)

Pertanto, la soluzione è

x_{1} (t) = \cos (2 t) + 2 \cos (t), x_{2} (t) = 4 \cos (t) - \cos (2 t) .

Essa può essere disegnata in forma parametrica 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)

Le singole componenti possono essere tracciate 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)

BIBLIOGRAFIA: Nagle, Saff, Snider, Fundamentals of Differential Equations, 6th ed, Addison-Wesley, 2004. (si veda § 5.5).

Metodo di Eulero per i sistemi di equazioni differenziali#

Nel prossimo esempio, si illustrerà il metodo di Eulero per le ODE di primo e secondo ordine. Per prima cosa ricordiamo l’idea di base per le equazioni di primo ordine. Dato un problema di Cauchy della forma

y^{'} = f (x, y) y (a) = c

si vuole trovare il valore approssimato della soluzione a $x = b$ con $b > a$ .

Ricordando dalla definizione di derivata che

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

dove $h > 0$ è dato e piccolo. Questo e la DE insieme danno give $f (x, y (x)) \approx \frac{y (x + h) - y (x)}{h}$ . Ora si risolve per $y (x + h)$ :

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

Se chiamiamo $h f (x, y (x))$ il «termine di correzione» (per mancanza di un termine migliore), $y (x)$ il «vecchio valore di y», e $y (x + h)$ il «nuovo valore di y», allora questa approssimazione può essere espressa come

y_{n e w} \approx y_{o l d} + h * f (x, y_{o l d}) .

Se si spezza l’intervallo da a a b in n intervalli, dimodoché $h = \frac{b - a}{n}$ , allora si possono registrare le informazioni per questo metodo in una tabella.

$x$	$y$	$h f (x, y)$
$a$	$c$	$h f (a, c)$
$a + h$	$c + h f (a, c)$	…
$a + 2 h$	…
…
$b = a + n h$	???	…

L’obiettivo è riempire tutti gli spazi vuoti della tavella, una riga alla volta, finché si arriva al valore ???, che è il metodo di approssimazione di Eulero per $y (b)$ .

L’idea per sistemi di ODE è simile.

Esempio: Si approssimi numericamente $z (t)$ a $t = 1$ usando 4 passi del metodo di Eulero, dove $z^{″} + t z^{'} + z = 0$ , $z (0) = 1$ , $z^{'} (0) = 0$ .

Si deve ridurre l’ODE di secondo ordine ad un sistema di due equazioni del primo ordine (usando $x = z$ , $y = z^{'}$ ) ed applicare il metodo di Eulero:

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

Pertanto, $z (1) \approx 0.75$ .

Si possono anche tracciare i punti $(x, y)$ per ottenere un grafico approssimato della curva. La funzione eulers_method_2x2_plot svolge questa funzione; per usarla, bisogna definire le funzioni f e g che prendono on argomento con tre coordinate: (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)

A questo punto, P ha in memoria due grafici: P[0], il grafico di x vs. t, e P[1], il grafico di y vs. t. Si possono tracciare entrambi come mostrato qui in seguito:

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

Funzioni speciali#

Sono implementati diversi polinomi ortogonali e funzioni speciali, usando sia PARI [GAP] che Maxima [Max]. Essi sono documentati nelle sezioni apposite («Polinomi ortogonali» e «Funzioni speciali», rispettivamente) del manuale di 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()
0.167089499251049

A questo punto, Sage ha soltanto incorporato queste funzioni per l’uso numerico. Per l’uso simbolico, si usi direttamente l’intefaccia di Maxima, come nell’esempio seguente:

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'

[GAP]

(en) The GAP Group, GAP - Groups, Algorithms, and Programming, Version 4.11; 2021, https://www.gap-system.org

[Max] (1,2)

(en) Maxima, Version 5.45; 2021, http://maxima.sf.net/