Базовая алгебра и вычисления

Sage может осуществлять вычисления такие, как поиск решений уравнений, дифференцирование, интегрирование и преобразования Лапласа. См. Sage Constructions , где содержатся примеры.

Решение уравнений

Точное решение уравнений

Функция solve решает уравнения. Для ее использования сначала нужно определить некоторые переменные; аргументами для solve будут уравнение (или система уравнений) и переменные, для которых нужно найти решение:

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)

Можно решать уравнения для одной переменной через другие:

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)

Также можно решать уравнения с несколькими переменными:

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)

Следующий пример показывает, как Sage решает систему нелинейных уравнений. Для начала система решается символьно:

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)

Для приближенных значений решения можно использовать:

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]

(Функция n выведет приближенное значение. Аргументом для данной функции является количество битов точности)

Численное решение уравнений

Во многих случаях функция solve не способна найти точное решение уравнения. Вместо нее можно использовать функцию find_root для нахождения численного решения. Например, solve не возвращает ничего существенного для следующего уравнения:

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

theta = var('theta')
solve(cos(theta)==sin(theta), theta)

С другой стороны функция find_root может использоваться для решения вышеуказанного примера в интервале \(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...

phi = var('phi')
find_root(cos(phi)==sin(phi),0,pi/2)

Дифференцирование, интегрирование и т.д.

Sage умеет дифференцировать и интегрировать многие функции. Например, для того, чтобы продифференцировать \(\sin(u)\) по переменной \(u\), требуется:

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)

Для подсчета четвертой производной функции \(\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)

Для нахождения частных производных, как, например, для функции \(x^2+17y^2\) по \(x\) и \(y\) соответственно:

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)

Теперь найдём интегралы: и определенные, и неопределенные. Например, \(\int x\sin(x^2)\, dx\) и \(\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)

Для нахождения разложения на простые дроби для \(\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)

Решение дифференциальных уравнений

Sage может использоваться для решения дифференциальных уравнений. Для решения уравнения \(x'+x-1=0\) сделаем следующее:

sage: t = var('t')    # определение переменной t для символьных вычислений
sage: x = function('x')(t)   # определение функции x зависящей от t
sage: DE = diff(x, t) + x - 1
sage: desolve(DE, [x,t])
(_C + e^t)*e^(-t)

>>> from sage.all import *
>>> t = var('t')    # определение переменной t для символьных вычислений
>>> x = function('x')(t)   # определение функции x зависящей от t
>>> DE = diff(x, t) + x - Integer(1)
>>> desolve(DE, [x,t])
(_C + e^t)*e^(-t)

t = var('t')    # определение переменной t для символьных вычислений
x = function('x')(t)   # определение функции x зависящей от t
DE = diff(x, t) + x - 1
desolve(DE, [x,t])

Для этого используется интерфейс Maxima [Max], поэтому результат может быть выведен в виде, отличном от обычного вывода Sage. В данном случае общее решение для данного дифференциального уравнения - \(x(t) = e^{-t}(e^{t}+C)\).

Преобразования Лапласа также могут быть вычислены. Преобразование Лапласа для \(t^2e^t -\sin(t)\) вычисляется следующим образом:

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)

Приведем более сложный пример. Отклонение от положения равновесия для пары пружин, прикрепленных к стене слева,

|------\/\/\/\/\---|масса1|----\/\/\/\/\/----|масса2|
        пружина1                пружина2

может быть представлено в виде дифференциальных уравнений второго порядка

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

где \(m_{i}\) - это масса объекта i, \(x_{i}\) - это отклонение от положения равновесия массы i, а \(k_{i}\) - это константа для пружины i.

Пример: Используйте Sage для вышеуказанного примера с \(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\).

Решение: Надо найти преобразование Лапласа первого уравнения (с условием \(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()

Данный результат тяжело читаем, однако должен быть понят как

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

Найдем преобразование Лапласа для второго уравнения:

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)

Результат:

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

Вставим начальные условия для \(x(0)\), \(x'(0)\), \(y(0)\) и \(y'(0)\), и решим уравения:

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)

Теперь произведём обратное преобразование Лапласа для нахождения ответа:

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)

Итак, ответ:

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

График для ответа может быть построен параметрически, используя

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)

Графики могут быть построены и для отдельных компонентов:

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)

Для более исчерпывающей информации по графикам см. Построение графиков. Также см. секцию 5.5 из [NagleEtAl2004] для углубленной информации по дифференциальным уравнениям.

Метод Эйлера для решения систем дифференциальных уравнений

В следующем примере показан метод Эйлера для дифференциальных уравнений первого и второго порядков. Сначала вспомним, что делается для уравнений первого порядка. Дана задача с начальными условиями в виде

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

требуется найти приблизительное значение решения при \(x=b\) и \(b>a\).

Из определения производной следует, что

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

где \(h>0\) дано и является небольшим. Это и дифференциальное уравнение дают \(f(x,y(x))\approx \frac{y(x+h)-y(x)}{h}\). Теперь надо решить для \(y(x+h)\):

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

Если назвать \(h\cdot f(x,y(x))\) «поправочным элементом», \(y(x)\) «прежним значением \(y\)» а \(y(x+h)\) «новым значением \(y\)», тогда данное приближение может быть выражено в виде

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

Если разбить интервал между \(a\) и \(b\) на \(n\) частей, чтобы \(h=\frac{b-a}{n}\), тогда можно записать информацию для данного метода в таблицу.

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

???

…

Целью является заполнить все пустоты в таблице по одному ряду за раз до момента достижения записи ???, которая и является приближенным значением метода Эйлера для \(y(b)\).

Решение систем дифференциальных уравнений похоже на решение обычных дифференциальных уравнений.

Пример: Найдите численное приблизительное значение для \(z(t)\) при \(t=1\), используя 4 шага метода Эйлера, где \(z''+tz'+z=0\), \(z(0)=1\), \(z'(0)=0\).

Требуется привести дифференциальное уравнение 2го порядка к системе двух дифференцальных уравнений первого порядка (используя \(x=z\), \(y=z'\)) и применить метод Эйлера:

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)

Итак, \(z(1)\approx 0.75\).

Можно построить график для точек \((x,y)\), чтобы получить приблизительный вид кривой. Функция eulers_method_2x2_plot выполнит данную задачу; для этого надо определить функции f и g, аргумент которых имеет три координаты: (\(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)

В этот момент P содержит в себе два графика: P[0] - график \(x\) по \(t\) и P[1] - график \(y\) по \(t\). Оба эти графика могут быть выведены следующим образом:

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

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

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

Специальные функции

Несколько ортогональных полиномов и специальных функций осуществлены с помощью PARI [GAP] и Maxima [Max].

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

На данный момент Sage рассматривает данные функции только для численного применения. Для символьного использования нужно напрямую использовать интерфейс Maxima, как описано ниже:

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