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

Sage может осуществлять вычисления такие, как поиск решений уравнений, 
дифференцирование, интегрирование и преобразования Лапласа. См. 
`Sage Constructions <http://www.sagemath.org/doc/constructions/>`_ , 
где содержатся примеры.

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

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

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

::

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

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

::

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

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

::

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

Следующий пример показывает, как 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(2)*sqrt(5) - 2/3],
     [p == 1, q == 8, x == 4/3*sqrt(10) - 2/3, y == -1/6*sqrt(2)*sqrt(5) - 2/3]]

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

.. link

::

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

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

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

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

    sage: theta = var('theta')
    sage: solve(cos(theta)==sin(theta), theta)
    [sin(theta) == cos(theta)]

С другой стороны функция ``find_root`` может использоваться для решения 
вышеуказанного примера в интервале :math:`0 < \phi < \pi/2`::

    sage: phi = var('phi')
    sage: find_root(cos(phi)==sin(phi),0,pi/2)
    0.785398163397448...

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

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

::

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

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

Для нахождения частных производных, как, например, для функции :math:`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

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

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

.. _section-systems:

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

Sage может использоваться для решения дифференциальных уравнений. Для 
решения уравнения :math:`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)

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

Преобразования Лапласа также могут быть вычислены. Преобразование Лапласа для 
:math:`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)
    2/(s - 1)^3 - 1/(s^2 + 1)

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

::

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

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

.. math::

    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,

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

**Пример:** Используйте Sage для вышеуказанного примера с 
:math:`m_{1}=2`, :math:`m_{2}=1`, :math:`k_{1}=4`,
:math:`k_{2}=2`, :math:`x_{1}(0)=3`, :math:`x_{1}'(0)=0`,
:math:`x_{2}(0)=3`, :math:`x_{2}'(0)=0`.

Решение: Надо найти преобразование Лапласа первого уравнения (с условием 
:math:`x=x_{1}`, :math:`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)

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

.. math:: -2x'(0) + 2s^2\cdot X(s) - 2sx(0) - 2Y(s) + 6X(s) = 0

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

::

    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

Результат:

.. math:: -Y'(0) + s^2Y(s) + 2Y(s) - 2X(s) - sy(0) = 0.

Вставим начальные условия для :math:`x(0)`, :math:`x'(0)`,
:math:`y(0)` и :math:`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)]]

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

::

    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)

Итак, ответ:

.. math:: 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)

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

::

    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)

Для более исчерпывающей информации по графикам см. :ref:`section-plot`. 
Также см. секцию 5.5 из [NagleEtAl2004]_ для углубленной информации по 
дифференциальным уравнениям.

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

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

.. math::

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

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

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

.. math::  y'(x) \approx \frac{y(x+h)-y(x)}{h},

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

.. math::   y(x+h) \approx y(x) + h\cdot f(x,y(x)).

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

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

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

============== =======================   =====================
:math:`x`      :math:`y`                 :math:`h\cdot f(x,y)`
============== =======================   =====================
:math:`a`      :math:`c`                 :math:`h\cdot f(a,c)`
:math:`a+h`    :math:`c+h\cdot f(a,c)`         ...
:math:`a+2h`   ...
...
:math:`b=a+nh` ???                             ...
============== =======================   =====================

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

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

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

Требуется привести дифференциальное уравнение 2го порядка к системе 
двух дифференцальных уравнений первого порядка (используя :math:`x=z`, 
:math:`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

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

Можно построить график для точек :math:`(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)

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

.. link

::

    sage: 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,"pari",250)
    0.56515910399248502720769602760986330732889962162109200948029448947925564096
    sage: bessel_I(1,1)
    0.565159103992485
    sage: bessel_I(2,1.1,"maxima")  # последние несколько цифр могут быть неточными
    0.167089499251049

На данный момент 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'