Введение в научный Python-1. Введение в научный Python

• Название: Введение в научный Python
• Дата: 09.01.2020
• Размер: 6.2 Mb.
• Имя файла: Введение в научный Python-1.pdf
• Тип: Документы; #103229
• страница: 16 из 22

evaluate float), который возвращает его десятичное представление (аналогично методу n()).
>>> sqrt(12).evalf()
# функция sqrt() модуля sympy
3.46410161513775
>>> E.evalf()
2.71828182845905
Метод evalf([n,...]) допускает использование аргумента, который задает точность результата (n = количество значащих цифр)
>>> sqrt(12).evalf(40)
3.464101615137754587054892683011744733886
>>> pi.evalf(20)
3.1415926535897932385
Функция N(expr,args) эквивалента команде simplify(expr).evalf(args), где simplify(expr) означает упрощение символьного выражения expr.
>>> N(sqrt(2),30)
1.41421356237309504880168872421
Сравните результаты следующих команд.
>>> N(sqrt(2)*pi)
4.44288293815837
>>> float(sqrt(2)*pi)
4.442882938158366
>>> complex(sqrt(2)*pi)
(4.442882938158366+0j)
Функция nsimplify(value,list,tolerance=1e-9) преобразует десятичное число value в ближайшее рациональное или в комбинацию рациональных чисел и констант, перечисленных в списке list.
>>> nsimplify(0.33333,tolerance=0.001)
1/3
>>> nsimplify(6.28,[pi],tolerance=0.01)
2*pi
>>> nsimplify(atan(1),[pi]) pi/4
Следует помнить об одной особенности вещественной арифметики. Обычно она не возвращает точный результат. Рассмотрим следующий пример.
>>> from sympy import *
>>> one=S('one')
>>> one = cos(1)**2 + sin(1)**2
>>> one.evalf()
# равно 1 1.00000000000000

183
Но
>>> (one-1).evalf()
# должно быть равно 0
-0.e-124
Если вы знаете, что результат содержит погрешность вычислений, то ее можно удалить с помощью опции chop=True метода evalf(). В этом случае очень маленькое значение вещественной или мнимой части результата заменяется нулем.
>>> (one - 1).evalf(chop=True)
0
Помните, что выражение, записанное в виде строки, отличается от символьного выражения. Функция sympify (не путать с simplify), о которой говорилось выше, используется для преобразования строки в символьное выражение.
>>> str = "x**4-8*x**2+16"
>>> expr = sympify(str)
>>> expr x**4-8*x**2+16
>>> expr.subs(x, 2)
0
>>> factor(expr)
# разложение на множители символьного выражения
(x-2)**2*(x+2)**2
Еще одним приемом, о котором следует знать, это способ извлечения частей символьных выражений.
>>> x,y,z=symbols('x y z')
>>> z=x+y+1
Атрибут args возвращает кортеж аргументов функции верхнего уровня.
>>> z.args
(1, x, y)
>>> z=x**2*y**2/4
>>> z.args
(1/4, x**2, y**2)
Метод as_ordered_terms() преобразует выражение в упорядоченный список операндов.
>>> z=x+y+1
>>> z.as_ordered_terms()
[x, y, 1]
После выделения кортежа или списка вы можете по индексу получить объект любого операнда.
После выполнения инструкции from sympy import * становится доступен специальный символ oo (две буквы „o‟) – бесконечность, с которым тоже можно выполнять некоторые операции.
>>> oo+1 oo
>>> 1000000< oo
True

184
>>> 1/oo
0
Символ oo (бесконечность) в основном используется функцией limit(), а также функцией integrate() при задании пределов интегрирования.
5.2 Алгебраические вычисления
Алгебраические преобразования. С символьными выражениями можно выполнять алгебраические преобразования. Для этого следует использовать подходящие функции пакета SymPy. Например, функция factor раскладывает алгебраическое выражение на множители.
>>> x,y = symbols('x y')
>>> factor(x**2-y**2)
(x - y)*(x + y)
Можно ли x
4
+4 разложить на множители без комплексных чисел. Давайте проверим.
>>> factor(x**4+4)
(x**2 - 2*x + 2)*(x**2 + 2*x + 2)
>>> expr=sum([ x**i for i in range(0, 6)]); expr x**5 + x**4 + x**3 + x**2 + x + 1
>>> factor(expr)
(x + 1)*(x**2 - x + 1)*(x**2 + x + 1)
Функция expand обратная к функции factor. Она раскрывает скобки в алгебраическом выражении.
>>> f=(a+b)**3
>>> expand(f) a**3 + 3*a**2*b + 3*a*b**2 + b**3
>>> expand((1 -x)*(1+x+x**2+x**3+x**4+x**5))
-x**6 + 1
Формула
 
y
x
y
x
log log log


справедлива для вещественных положительных x и y, но для отрицательных и комплекснозначных x и y она неверна. Поэтому следующее разложение не срабатывает.
>>> expand(log(x*y)) log(x*y)
Можно проигнорировать подобное ограничение и заставить функцию выполнить привычное преобразование.
>>> expand(log(x* y), force=True) log(x) + log(y)
Многие другие функции, выполняющие алгебраические преобразования, поддерживают опцию force=True. Однако лучшим решением является наложение подходящих ограничений на символьные переменные.
>>> var('a,b', positive=True)
(a, b)
>>> expand(log(a*b)) log(a) + log(b)

185
>>> expand(prod([(x -i) for i in range(1,5)])) x**4 - 10*x**3 + 35*x**2 - 50*x + 24
Здесь использована функция prod() модуля sympy, которая вычисляет произведение своих аргументов.
>>> from sympy import prod, S
>>> prod([S(5), S(3)])
15
>>> type(_)

>>> prod(range(1,5))
24
>>> prod([4, 5], 2)
40
В выражении собрать коэффициенты при одинаковых степенях переменной можно с помощью функции collect(expr,var).
>>> x,y,z=symbols('x y z')
>>> expr = x*y + x - 3 + 2*x**2 - z*x**2 + x**3
>>> collect(expr, x)
x**3 + x**2*(-z + 2) + x*(y + 1) - 3
>>> a,b,c,d=symbols('a b c d')
>>> collect(a*x+b*x+c*y+d*y,x) c*y + d*y + x*(a + b)
>>> collect(a*x+b*x+c*y+d*y,(x,y)) x*(a + b) + y*(c + d)
Функция cancel() приводит рациональное выражение к общему знаменателю и сокращает числитель и знаменатель на общий множитель, если он есть.
>>> cancel(a/b+c/d)
(a*d + b*c)/(b*d)
>>> cancel((x**2 + 2*x + 1)/(x**3 + x**2))
(x + 1)/x**2
>>> cancel((x - 7)/(x +3)-(x-12)/(x-2))
50/(x**2 + x - 6)
Наоборот, функция apart() выполняет разложение дробей.
>>> apart((4*x**3 + 21*x**2 + 10*x + 12)/(x**4 + 5*x**3 + 5*x**2 + 4*x))
(2*x - 1)/(x**2 + x + 1) - 1/(x + 4) + 3/x
Универсальная функция simplify() пытается упрощать алгебраические выражения, используя различные комбинации приведенных выше функций.
>>> from sympy import *
>>> var('x')
>>> simplify(sin(x)**2 + cos(x)**2)
1
>>> a,b = symbols('a b')
>>> z=simplify((a**4 - 2*a**2*b**2 + b**4)/(a -b)**2);z a**2 + 2*a*b + b**2
Как видите, функция simplify не до конца справилась со своей задачей.
Поэтому ей требуется помощь.

186
>>> factor(z)
(a + b)**2
Имеются функции, предназначенные для преобразования тригонометрических выражений.
>>> expr=sin(2*x)+cos(3*x)
>>> expand_trig(expr)
2*sin(x)*cos(x) + 4*cos(x)**3 - 3*cos(x)
>>> tr=expand_trig(cos(4*x) -2*cos(x)*cos(3*x)); tr
-2*(4*cos(x)**3-3*cos(x))*cos(x)+8*cos(x)**4-8*cos(x)**2+1
>>> trigsimp(tr)
2*(-cos(2*x) + 1)**2 + 3*cos(2*x) - cos(4*x) - 3
>>> simplify(tr)
-cos(2*x)
>>> trigsimp(sin(x)**4 - 2*cos(x)**2*sin(x)**2 + cos(x)**4) cos(4*x)/2 + 1/2
Символьные выражения являются объектами и имеют различные методы
(аналогичные функциям), позволяющие их преобразовывать.
>>> a,b,c,d,x,y,z = symbols('a b c d x y z')
>>> expr=6+11*x-3*x**2-2*x**3
>>> z=expr.factor();z
-(x - 2)*(x + 3)*(2*x + 1)
>>> z.expand()
-2*x**3 - 3*x**2 + 11*x + 6
>>> ((a+b+c)**2).expand() a**2 + 2*a*b + 2*a*c + b**2 + 2*b*c + c**2
>>> w=a/b+c/d
>>> u=w.together(); u
(a*d + b*c)/(b*d)
>>> u.simplify() a/b + c/d
>>> z=a*x+b*x+c*y+d*y
>>> z.collect(x)
c*y + d*y + x*(a + b)
>>> z.collect((x,y)) x*(a + b) + y*(c + d)
Комплексная арифметика. Для ввода комплексных чисел можно использовать стандартную нотацию, используя для мнимой единицы символ
„I‟ (заглавная буква I).
>>>
I*I
-1
>>> a,b =symbols('a b',real=True)
>>> z=symbols('z')
>>> z=a+b*
I
Функции re(z) и im(z) вычисляют вещественную и мнимую часть комплексного числа z.

187
>>> re(z ) a
>>> im(z ) b
Функция Abs(z) и arg(z) вычисляют модуль и аргумент комплексного числа z.
>>> Abs(z) sqrt(a**2 + b**2)
>>> arg(1+2*
I
) atan(2)
Функция conjugate(z) выполняет комплексное сопряжение.
>>> conjugate(z) a - I*b
Арифметические операции с комплексными числами выполняются обычным образом.
>>> z1=4+3*I
>>> z2=1-2*I
>>> z1+z2 5 + I
>>> simplify(z1*z2)
10 - 5*I
>>> simplify(z1/z2)
-2/5 + 11*I/5
>>> var('w')
>>> w=z**3;
>>> w=expand(w);
>>> collect(w,I)
a**3 - 3*a*b**2 + I*(3*a**2*b - b**3)
>>> simplify(re(w)) a*(a**2 - 3*b**2)
>>> simplify(im(w)) b*(3*a**2 - b**2)
Cкажем несколько слов об отображении формул. Среда выполнения IPython и
IPython console в Spyder поддерживают инструкцию init_printing(...), которая позволяет настраивать способ отображения математических выражений в своем окне. Наиболее часто используются опции: pretty_print= True / False use_unicode = True / False use_latex = True / False
Например, инструкция init_printing(use_latex=True) включает Latex режим отображения. С другой стороны, инструкция init_printing(pretty_print=False) включает стандартный текстовый режим отображения.
В нашем пособии для представления результатов вычисления мы, обычно, используем текстовый режим. Но иногда, для отображения формул,

188 будем использовать Latex режим (в среде IDLE он недоступен). Вот, например, как выглядит решение квадратного уравнения в режиме Latex. from sympy import * init_printing(use_latex=True) a,b,c,x = symbols('a b c x') solve(a*x**2+b*x+c,x)
Вернемся к комплексной арифметике. SymPy умеет работать с комплексными функциями. var('x',real=True) exp(I*x).expand(complex=True)
I*sin(x) + cos(x)
В IPython, включив Latex режим, вы получите результат в привычной математической записи. init_printing(use_latex=True) exp(I*x).expand(complex=True)
re(_)
Решение уравнений. Большинство алгебраических преобразований мы можем выполнить самостоятельно. Однако решение уравнений является боле сложной задачей. И в этом вопросе помощь SymPy может быть значительной. Здесь мы рассматриваем функции пакета, с помощью которых можно решать алгебраические уравнения символьно. В тех случаях, когда функции пакета
SymPy не могут получить решение, вам следует обращаться к функциям других пакетов, реализующих подходящие численные алгоритмы.
Для записи уравнений в SymPy не используются знаки равенства ‘=’ или двойного равенства ‘==’, а используется функция Eq(expr1,expr2). Если можно проверить, равны ли выражения expr1 и expr2, то функция Eq() возвращает True.
>>> from sympy import *
>>> Eq(3, 4)
False
У функции Eq() есть опция evaluate=False, которая принудительно может отменить проверку равенства.
>>> Eq(3, 4,evaluate=False)
Eq(3, 4)
>>> x,y=symbols('x y')
>>> Eq(x,y)
Eq(x, y)
Функции, которые предназначены для решения уравнений, используют функцию Eq() для записи этих уравнений. Кроме того, в SymPy все функции
«решатели» в том месте, где может стоять функция Eq(), а стоит символьное выражение expr, автоматически полагают, что оно равно нулю (т.е.

189 используется Eq(expr,0)). Поэтому в «решателях» вместо уравнения вида
Eq(expr1,expr2) можно использовать символьное выражение expr1-expr2.
Функцией, предназначенной для символьного решения алгебраического уравнения, является solveset(). Она имеет следующий синтаксис. solveset(equation,variable=None[,domain=S.Complexes]), где equation может быть записано в форме Eq() или быть выражением, которое полагается равным нулю. Второй аргумент задает имя искомой переменной, третий – определяет область, в которой ищутся корни
(вещественная ось или комплексная плоскость). Результат функция возвращает в виде множества.
>>> solveset(Eq(x**2, 4), x)
{-2, 2}
>>> solveset(Eq(x**2 - 7*x+12, 0), x)
{3, 4}
>>> solveset(Eq(x**2 - 5*x+12, 0), x)
{5/2 - sqrt(23)*I/2, 5/2 + sqrt(23)*I/2}
>>> type(_)

Для одного уравнения результатом будет объект типа FiniteSet.
>>> solveset(x**4 -1, x)
{-1, 1, -I, I}
>>> type(_)

В некоторых случаях результат может иметь другой тип.
>>> solveset(x**2 - x**2, x, domain=S.Reals)
(-oo, oo)
>>> type(_)

>>> solveset(cos(x) - 1, x, domain=S.Reals)
ImageSet(Lambda(_n, 2*_n*pi), Integers())
Если решения нет, то возвращается пустое множество.
>>> solveset(exp(x), x)
EmptySet()
Если решения найти не удалось, то функция возвращает ConditionSet.
В пакете есть функция real_roots(), которая возвращает только вещественные корни.
>>> x=S('x')
>>> real_roots(x**2-4)
[-2, 2]
>>> p = (x-3)**2*(x-2)*(x-1)*x*(x+1)*(x**2 + x + 1)
>>> real_roots(p)
[-1, 0, 1, 2, 3, 3]
В последнем примере два комплексных корня множителя (x**2 + x + 1) не были учтены функцией real_roots().

190
>>> real_roots(x**2+x+1)
[]
Функция roots() возвращает символьное выражение для корней полинома одной переменной. Результат возвращается в форме словаря с парами корень:кратность.
>>> a,b,c,x=symbols('a,b,c,x')
>>> roots(x**2 -c**2,x)
{-c: 1, c: 1}
>>> roots(x**3 - 5*x**2 + 3*x + 9)
{3: 2, -1: 1}
Если последнее выражение разложить на множители, то получим
>>> factor(x**3 - 5*x**2 + 3*x + 9)
(x - 3)**2*(x + 1)
Имеется функция nroots(p[,n=15,...]), которая находит приближенное значение корней полинома p с точностью n значащих цифр.
>>> nroots(x**2-2, 20)
[-1.4142135623730950488, 1.4142135623730950488]
Для построения аналитических решений алгебраических уравнений и систем предназначена функция solve(...). Первым аргументом она также принимает уравнение (или список уравнений) в форме Eq() или выражение (список выражений), которое полагается равным нулю.
>>> x,y = symbols('x y')
>>> solve(x**2+x-30,x)
[-6, 5]
Второй аргумент – имя искомой переменной необязателен, но его лучше всегда указывать. Рассмотрим пример.
>>> a,b,c,x = symbols('a b c x')
>>> eq = Eq(a*x**2+b*x+c, 0)
>>> solve(eq)
[{a: -(b*x + c)/x**2}]
Уравнение решено относительно переменной a. Поэтому лучше всегда указывать имя искомой переменной.
>>> solve(eq,x)
[(-b+sqrt(-4*a*c+b**2))/(2*a),-(b+sqrt(-4*a*c+b**2))/(2*a)]
В следующем примере решим систему линейных уравнений.







7 3
0 2
y
x
y
x
>>> solve([x-2*y,3*x+y-7],x,y)
{y:1,x:2}
Форма возвращаемого результата зависит от поставленной задачи и может быть списком, словарем, множеством и т.д. Кроме того, форма результата может быть задана опциями. Например, опция dict=False отменяет вывод результата в форме словаря (и включает вывод в форме списка).
>>> solve([Eq(x**2-y**2,0),Eq(x+y,2)],x,y,dict=False)
[(1, 1)]

191
Функция solve() может решать некоторые трансцендентные уравнения.
>>> solve(exp(x**2+x -6) - 1, x)
[-3, 2]
>>> solve(exp(x) + 1, x)
# комплексный корень
[I*pi]
Функция solve() имеет много вариантов использования, с которыми вы можете познакомиться самостоятельно по справочной системе.
Для решения систем линейных уравнений предназначена функция linsolve(...). Запишем последнюю систему уравнений в виде








0 7
3 0
2
y
x
y
x
,
Функции linsolve(...) можно передавать левые части уравнений этой системы.
>>> linsolve([x-2*y,3*x+y-7],x,y)
{(2, 1)}
Другой способ решения этой задачи с помощью функции linsolve() состоит в использовании расширенной матрицы системы







7 1
3 0
2 1
>>> linsolve(Matrix(([1,-2,0],[3,1,7])),( x,y))
{(2, 1)}
Матричная алгебра. Символьные матрицы создаются функцией Matrix() и являются объектами, у которых имеются атрибуты и методы.
>>> from sympy import *
>>> a,b,c,d = symbols('a b c d')
>>> M=Matrix([[a,b],[c,d]])
При работе с символьными матрицами следует убедиться, что они (матрицы) отображаются в вашей среде корректно. Например, в IPython, если просто напечатать имя матрицы в режиме init_printing(use_latex=True), происходит сбой и матрица не отображается. Для печати матрицы ее имя следует передать аргументом функции «милой» печати pprint(имя_матрицы).
Можно также отключить
Latex режим следующей командой init_printing(use_latex=False). Тогда, чтобы напечатать матрицу, в командной строке можно просто ввести ее имя и выполнить инструкцию (без использования функций печати).
Доступ к элементам символьных матриц выполняется по индексам с использованием квадратных скобок.
>>> M[0,0] a
Использование одномерного списка в аргументе функции Matrix() приводит к созданию вектора–столбца (матрицы с одним столбцом).
>>> v1=Matrix([1,2,3]); v1

192
Matrix([
[1],
[2],
[3]])
После создания матрицы вы можете вызывать методы, которые будут возвращать затребованные преобразования.
>>> M.det()
# вычисление определителя a*d - b*c
>>> M.T
# транспонирование матрицы
Matrix([
[a, c],
[b, d]])
>>> B=simplify(M.inv()) ; B
# вычисление обратной матрицы
Matrix([
[ d/(a*d - b*c), -b/(a*d - b*c)],
[-c/(a*d - b*c), a/(a*d - b*c)]])
>>> simplify(M*B)
# проверка
Matrix([
[1, 0],
[0, 1]])
>>> M=Matrix([[1,4],[2,3]])
>>> M.eigenvals()
# вычисление собственных чисел
{5: 1, -1: 1}
# с.число 5 имеет кратность 1; с.число -1 имеет кратность 1
>>> M.eigenvects()
[(-1, 1, [Matrix([[-2],[ 1]])]),
( 5, 1, [Matrix([[1],[1]])])]
Здесь видно, что собственному числу -1 с кратностью 1 соответствует собственный вектор (-2,1); собственному числу 5 с кратностью 1 соответствует собственный вектор (1,1).
Операции между символьными матрицами выполняются с помощью привычных операторов.
>>> x=symbols('x:2:2'); x
(x00, x01, x10, x11)
>>> y=symbols('y:2:2') ; y
(y00, y01, y10, y11)
>>> A=Matrix([[x[0],x[1]],[x[2],x[3]]]); A
Matrix([
[x00, x01],
[x10, x11]])
>>> B=Matrix([[y[0],y[1]],[y[2],y[3]]])
>>> A+B
Matrix([
[x00 + y00, x01 + y01],
[x10 + y10, x11 + y11]])
>>> A*B

193
Matrix([
[x00*y00 + x01*y10, x00*y01 + x01*y11],
[x10*y00 + x11*y10, x10*y01 + x11*y11]])
>>> A**2
Matrix([
[ x00**2 + x01*x10, x00*x01 + x01*x11],
[x00*x10 + x10*x11, x01*x10 + x11**2]])
>>> A**( -1)
# вычисление обратной матрицы
Matrix([
[1/x00+x01*x10/(x00**2*(x11-x01*x10/x00)),
-x01/(x00*(x11- x01*x10/x00))],
[ -x10/(x00*(x11 - x01*x10/x00)), 1/(x11 - x01*x10/x00)]])
>>> A=Matrix([[1,2],[3,4]])
>>> B=Matrix([[1,-1],[-1,3]])
>>> A/B
# умножение на обратную матрицу (A*B
-1
)
Matrix([
[ 5/2, 3/2],
[13/2, 7/2]])
>>> A**(-1)
# вычисление обратной матрицы
Matrix([
[ -2, 1],
[3/2, -1/2]])
Многие операции можно выполнить, используя методы объектов «матрица».
>>> v1=Matrix([1,2,3]);
>>> v2=Matrix([1,1,1])
>>> v1.dot(v2)
# скалярное произведение
6
>>> v1.cross(v2)
# векторное произведение
Matrix([
[-1],
[ 2],
[-1]])
В модуле SymPy имеются конструкторы специальных матриц: нулевых (zeros), единичных (ones), диагональных (diag) и т. д, а также большое количество других функций, выполняющих матричные вычисления.