Базовая алгебра и вычисления¶
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(10) - 2/3], [p == 1, q == 8, x == 4/3*sqrt(10) - 2/3, y == -1/6*sqrt(10) - 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)
-1/(s^2 + 1) + 2/(s - 1)^3
Приведем более сложный пример. Отклонение от положения равновесия для пары пружин, прикрепленных к стене слева,
|------\/\/\/\/\---|масса1|----\/\/\/\/\/----|масса2|
пружина1 пружина2
может быть представлено в виде дифференциальных уравнений второго порядка
где \(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)
Данный результат тяжело читаем, однако должен быть понят как
Найдем преобразование Лапласа для второго уравнения:
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
Результат:
Вставим начальные условия для \(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)
Итак, ответ:
График для ответа может быть построен параметрически, используя
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] для углубленной информации по дифференциальным уравнениям.
Метод Эйлера для решения систем дифференциальных уравнений¶
В следующем примере показан метод Эйлера для дифференциальных уравнений первого и второго порядков. Сначала вспомним, что делается для уравнений первого порядка. Дана задача с начальными условиями в виде
требуется найти приблизительное значение решения при \(x=b\) и \(b>a\).
Из определения производной следует, что
где \(h>0\) дано и является небольшим. Это и дифференциальное уравнение дают \(f(x,y(x))\approx \frac{y(x+h)-y(x)}{h}\). Теперь надо решить для \(y(x+h)\):
Если назвать \(h\cdot f(x,y(x))\) «поправочным элементом», \(y(x)\) «прежним значением \(y\)» а \(y(x+h)\) «новым значением \(y\)», тогда данное приближение может быть выражено в виде
Если разбить интервал между \(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
Итак, \(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).n(250)
0.56515910399248502720769602760986330732889962162109200948029448947925564096
sage: bessel_I(1,1).n()
0.565159103992485
sage: bessel_I(2,1.1).n()
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'