Sage может осуществлять вычисления такие, как поиск решений уравнений, дифференцирование, интегрирование и преобразования Лапласа. См. Sage 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]]
Для приближенных значений решения можно использовать:
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 может использоваться для решения вышеуказанного примера в интервале 0 < \phi < \pi/2:
sage: phi = var('phi')
sage: find_root(cos(phi)==sin(phi),0,pi/2)
0.785398163397448...
Sage умеет дифференцировать и интегрировать многие функции. Например, для того, чтобы продифференцировать \sin(u) по переменной u, требуется:
sage: u = var('u')
sage: diff(sin(u), u)
cos(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)
Для нахождения частных производных, как, например, для функции 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
Теперь найдём интегралы: и определенные, и неопределенные. Например, \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)
Для нахождения разложения на простые дроби для \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)
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)
Для этого используется интерфейс 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)
2/(s - 1)^3 - 1/(s^2 + 1)
Приведем более сложный пример. Отклонение от положения равновесия для пары пружин, прикрепленных к стене слева,
|------\/\/\/\/\---|масса1|----\/\/\/\/\/----|масса2|
пружина1 пружина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,
где 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
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)
Данный результат тяжело читаем, однако должен быть понят как
-2x'(0) + 2s^2*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
Результат:
-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)]]
Теперь произведём обратное преобразование Лапласа для нахождения ответа:
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)
Итак, ответ:
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)
Для более исчерпывающей информации по графикам см. Построение графиков. Также см. секцию 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*f(x,y(x)).
Если назвать h f(x,y(x)) “поправочным элементом”, y(x) “прежним значением y” а y(x+h) “новым значением y”, тогда данное приближение может быть выражено в виде
y_{new} \approx y_{old} + h*f(x,y_{old}).
Если разбить интервал между a и b на n частей, чтобы h=\frac{b-a}{n}, тогда можно записать информацию для данного метода в таблицу.
x | y | hf(x,y) |
---|---|---|
a | c | hf(a,c) |
a+h | c+hf(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
Итак, 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)
В этот момент P содержит в себе два графика: P[0] - график x по t и P[1] - график y по t. Оба эти графика могут быть выведены следующим образом:
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'