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

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

0.99702 | / \
| / \
| / \
| . .
|
| . .
|
| . .
0.49857 | --------------------.--------------.-------------------
|
| . .
| . .
|
| . .
| / \
| .. ..
| .. ..
0.00012 | ............ ...........
-3 0 3
5.5
Символьное решение дифференциальных уравнений
Ранее мы говорили, что для записи уравнений SymPy используется функция
Eq()
. Фактически запись Eq(выражение1,выражение2) создает специальный объект, который передается функциям, применяемым для решения дифференциальных уравнений (солверам).
>>> from sympy import *
>>> x = symbols('x')
>>> y=Function('y')
# объявление символьной функции
>>> pprint(dsolve(Eq(y(x).diff(x,x)+y(x),exp(x)),y(x)))
#
x
e
y
y



x e y(x) = C1*sin(x) + C2*cos(x) + --
2
Решением является выражение
x
e
x
C
x
C
2 1
cos sin
2 1


. Обратите внимание, что в решении появились две неопределенные константы C1 и C2.
Для решения обыкновенных дифференциальных уравнений (ОДУ) вначале объявляется символьная функция, которая будет представлять решение. Для этого имя этой функции определяется как символьная переменная в функции symbols() с опцией cls=Function. Например, команда f = symbols('f', cls=Function) создает одну такую функцию.
Можно использовать эквивалентную ей запись f=Function('f'). После этого выполняется инструкция решения ДУ. Объект Eq() может создаваться сразу внутри вызова функции dsolve(), или отдельно перед ее вызовом.
Решим ОДУ
x
f
f
f
sin
2





. Вначале создаем символьную функцию f
, а затем объект diffeq, представляющий уравнение.
>>> f = symbols('f', cls=Function)
>>> diffeq = Eq(f(x).diff(x, x) - 2*f(x).diff(x) + f(x), sin(x))

209
>>> pprint(diffeq)
2 d d f(x) - 2*--(f(x)) + ---(f(x)) = sin(x) dx 2 dx
Для решения используется функция dsolve(diffeq,f(x)) с первым аргументом – объектом уравнения. Вторым аргументом является искомое выражение.
>>> dsolve(diffeq, f(x))
Eq(f(x), (C1 + C2*x)*exp(x) + cos(x)/2)
Последний пример иллюстрирует обычную последовательность инструкций, используемых для решения ОДУ. Она состоит в следующем:

объявляется символьная независимая переменная x и символьная функция f;

создается объект, представляющий уравнение;

решается дифференциальное уравнение с помощью функции dsolve().
Решим несколько примеров.
>>> x=Symbol('x')
>>> f=Function('f')
>>> eq=Eq(f(x).diff(x,x)+f(x),0)
>>> pprint(dsolve(eq,f(x)))
#
0



f
f
f(x) = C1*sin(x) + C2*cos(x)
Договоримся, что далее решать ОДУ будем в IPython в режиме отображения
Latex, т.е. будем считать, что во всех последующих примерах этого параграфа включен режим «красивой» печати. init_printing(use_latex=True)
Функция dsolve() может использовать несколько различных методов решения ОДУ. Чтобы получить список методов, которые можно использовать для решения уравнения, следует использовать инструкцию classify_ode(уравнение, выражение/функция). classify_ode(eq, f(x))
('separable',
'1st_exact',
'almost_linear',
'1st_power_series',
'lie_group',
'separable_Integral',
'1st_exact_Integral',
'almost_linear_Integral')
Опция hint=‟method_name‟ функции dsolve() задает метод решения. f = Function('f') eq = sin(x)*cos(f(x)) + cos(x)*sin(f(x))*f(x).diff(x) rez=dsolve(eq, hint='1st_exact') ; rez

210 rez=dsolve(eq, hint='separable_Integral'); rez
Опции hint можно присвоить значение all . Тогда вы получите словарь решений в форме метод:решение. rez=dsolve(eq, hint='all');
Значение all_Integral также возвращает все интегралы ОДУ. Опции hint можно присвоить значение „best‟. Тогда функция вернет простейшее (по ее мнению) решение ОДУ. rez=dsolve(eq, hint='best'); rez
Опция ics позволяет задавать граничные условия. Она должна принимать форму {f(x0): y0,f(x).diff(x).subs(x, x1):y1[,...]}. В настоящий момент эта опция реализована только для метода решения в виде рядов для
ОДУ 1 – го порядка (hint='1st_power_series'). В этом случае условие должно принимать вид {f(x0): y0}. При решении в виде рядов опция n определяет порядок независимой переменной, до которого строится степенной ряд решения. Опция simplify позволяет использовать символьные алгоритмы преобразования полученного выражения.
Пример. Решим задачу Коши
 
 
2 2






u
t
u
t
u
,
 
 
1 0
,
0 0



u
u
import numpy as np from sympy import * from IPython.display import * import matplotlib.pyplot as plt init_printing(use_latex=True) var('t C1 C2') u = Function("u")(t)
# здесь u – выражение (не функция)
de = Eq(u.diff(t,t) -u.diff(t) -2*u,-2) display(de)
Функция display() отображает объекты Python в любой его оболочке. Она загружается из модуля IPython.display. Вид результата зависит от оболочки и, возможно, режима отображения. Здесь мы просто напечатали наше дифференциальное уравнение.
Теперь решаем ДУ и печатаем результат. des = dsolve(de,u)
# здесь u – выражение, а не функция display(des)

211
Найдем постоянные C1 и C2, которые удовлетворяют начальным условиям.
Вначале выполним подстановку t=0 в правую часть полученного выражения des
. Здесь метод rhs (right hand side) выбирает правую часть символьного выражения. eq1=des.rhs.subs(t,0); eq1
Дифференцируем правую часть выражения des, и делаем подстановку t=0. eq2=des.rhs.diff(t).subs(t,0) eq2
Выражения eq1 и eq2 должны равняться 0 и 1 (начальным значениям). Решаем систему уравнений относительно неизвестных C1 и C2. seq=solve([eq1,eq2 -1],C1,C2) seq
Строим результирующее выражение, которое представляет решение нашей задачи, и «красиво» его печатаем. rez=des.rhs.subs([(C1,seq[C1]),(C2,seq[C2])])
F = Lambda(t,rez) display(Latex('$u(t) = ' + str(latex(F(t))) + '$'))
Здесь функция sympy.Lambda(x,expr) создает символьную лямбда – функцию аналогично лямбда – функции Python 'lambda x: expr'. Функция нескольких переменных создается с использованием синтаксиса
Lambda((x,y,...),expr).
Функция sympy.latex() преобразует символьное выражение в строку
TeX кода. А функция Latex() (из модуля IPython.display) преобразует строку, содержащую TeX код, в объект, который может отобразить функция display().
Затем преобразуем символьное выражение в лямбда–функцию, поддерживающую работу с массивами модуля numpy. f=lambdify(t, rez, "numpy")
Теперь, используя функции модуля matplotlib.pyplot, строим график решения. x = np.linspace(0,5,100) plt.grid(True) plt.xlabel('Time t seconds',fontsize=12) plt.ylabel('$f(t)$',fontsize=16) plt.plot(x,f(x),color='#008000', linewidth=3) plt.show()

212
■
Приведем примеры решения систем дифференциальных уравнений. При этом, для обозначения производных в уравнениях можно использовать как функцию
Derivative()
, так и метод diff(). Заметим также, что в функции dsolve указывать имена искомых функций не требуется – это будет выполняться автоматически.
Пример. Решим систему ОДУ
 
 










t
y
t
y
t
y
x
t
x
sin sin
2
# включен режим init_printing(use_latex=True)
from sympy import * t = symbols('t') x, y = symbols('x, y', function=True) eq = (Eq(Derivative(x(t),t),x(t)*y(t)*sin(t)), \
Eq(Derivative(y(t),t),y(t)**2*sin(t))) rez=dsolve(eq) display(list(rez))
Пример. Решим систему ОДУ
 
 
 
 









x
y
x
z
x
z
x
y
from sympy import * x = symbols('x') y, z = symbols('y, z', function=True) eq = (Eq(Derivative(y(x),x), -z(x)), Eq(Derivative(z(x),x), -y(x)))
rez=dsolve(eq)
display(list(rez))
Пример. Решим систему ОДУ
 
 
 
















z
y
t
z
z
y
t
y
z
x
t
x
t = symbols('t') x,y,z = symbols('x,y,z', function=True) eq1=Eq(x(t).diff(t),-x(t)+z(t))

213 eq2=Eq(y(t).diff(t),-y(t)-z(t)) eq3=Eq(z(t).diff(t),y(t) -z(t)) rez=dsolve((eq1,eq2,eq3)) display(rez[0]) display(rez[1]) display(rez[2])
■
В пакете sympy имеется модуль sympy.solvers.pde, содержащий функции, предназначенные для символьного решения дифференциальных уравнений в частных производных (ДУЧП). Поскольку таких уравнений известно немного, то возможности этого модуля ограничены. Основной функцией этого пакета является функция sympy.solvers.pde.pdsolve(eq,...), где eq является любым допустимым ДУЧП.
Пример. Решим ДУЧП
y
x
y
u
x
u







1
from sympy import * from sympy.solvers.pde import pdsolve f = Function('f') u = f(x, y) ux = u.diff(x) uy = u.diff(y) eq=Eq(ux+uy,1/(x*y)) pdsolve(eq)
Общее решение дифференциального уравнения в частных производных, если оно вообще может быть найдено, должно включать в себя произвольные функции. В последнем решении такой функцией является


y
x
F

Пока неизвестная функция и ее производные входят в уравнение линейно можно надеяться получить общее решение. f = Function('f') u = f(x, y) ux = u.diff(x) uy = u.diff(y) eq=Eq(x*ux+y*uy,exp(x*y)) pdsolve(eq)
Здесь F(...) произвольная функция, а Ei(x) - интегральная показательная функция.

214
Если общее решение не может быть найдено, то система выводит сообщение
NotImplementedError: psolve: Cannot solve ...
В пакете SymPy имеется большое количество функций, которые предназначены для решения специальных типов обыкновенных дифференциальных уравнений. Однако, в настоящее время, возможности символьного решения дифференциальных уравнений весьма ограничены и на практике чаще используются функции, реализующие численные методы. С ними мы познакомимся в следующей главе.
5.6
Совместное использование символьной и численной математики
Имеется несколько причин включать символьные вычисления в числовой код.
Например, это случается тогда, когда мы хотим сравнить известное аналитическое решение с решением, полученным на сетке. Довольно часто перед выполнением численных расчетов, нам надо выполнить какие–либо символьные преобразования, например, вычислить производные аналитически заданной функции. То же касается вычисления интегралов или дифференциальных уравнений, если конечно их решения можно получить в аналитическом виде. В любом случае мы можем вставлять символьные выражения в код, и использовать их не только для вычисления значений.
В этом параграфе мы опишем способы использования возможностей пакета SymPy совместно с другими пакетами научного Python. Они сводятся к преобразованию символьных выражений в функции Python или в функции, способные работать с данными другого модуля. Встречаются ситуации, когда из функции Python нужно получить символьное выражение.
Пусть дано символьное выражение expr.
>>> x=symbols('x')
>>> expr=x**4+4
Преобразуем его в функцию Python.
>>> fn=lambda t : expr.subs({x:t})
>>> fn(1)
5
Достаточно выполнить подстановку аргумента лямбда – функции t в символьное выражение expr в теле функции.
Подстановка subs() в символьное выражение удобна, если его значение надо вычислить в одной точке. Но если вы намерены вычислять значение выражения во многих точках, то имеется более эффективное решение –
функция lambdify(). Она возвращает лямбда – функцию и при этом конвертирует SymPy имена (функций) в имена функций подходящего научного модуля, например, numpy. Формат вызова этой функции следующий: fun=lambdify(args, expr, modules=None[,...])
Здесь args является именем аргумента лямбда – функции или кортежем имен.
Второй параметр expr является символьным выражением, которое будет преобразовываться. Третий параметр modules, если он задан, может быть

215 строковым именем одного из модулей: „math‟, ‟mpmath‟, ‟numpy‟,
„numexpr‟, ‟sympy‟. Третий аргумент может быть также словарем с парами замен вида sympy_name : func_name. Если третий аргумент не задан, то выполняется замена (насколько это возможно), функциями модуля math, numpy или mpmath (в таком порядке).
>>> import numpy as np
>>> from sympy import *
>>> x=symbols('x')
>>> expr = sin(x)+cos(x)
>>> g= lambdify(x, expr, "numpy")
>>> f=np.arange(5)
>>> g(f)
array([ 1. ,1.38177329, 0.49315059, -0.84887249, -1.41044612])
Таким образом, если мы хотим, чтобы преобразованная из символьного выражения функция, работала с массивами модуля numpy, то ее следует преобразовать в лямбда – функцию, используя опцию modules=‟numpy‟ (в модуле numpy свои элементарные и специальные функции, которые умеют работать с массивами).
Рассмотрим теперь обратную задачу преобразования функции Python в символьное выражение.
Пусть выражение, возвращаемое функцией, содержит только стандартные математические операции над аргументом функции. Тогда для получения символьного выражения, вычисляемого по той же формуле, достаточно вызвать функцию с символьным аргументом. Например, пусть дана функция Python.
>>> fn=lambda t: t**4+4
>>> fn(2)
20
Для получения символьного выражения достаточно вызвать эту функцию с символьным аргументом.
>>> x=symbols('x')
>>> z=fn(x)
>>> type(z)

>>> factor(fn(x))
# разложение на множители символьного выражения
(x**2 - 2*x + 2)*(x**2 + 2*x + 2)
Как видите, результатом вызова функции fn с символьным аргументом является символьное выражение, с которым можно работать как с любым символьным выражением. В частности, результат можно дифференцировать.
>>> def f(x): return x**4
>>> x=S('x')
>>> f(x).diff(x)
4*x**3
>>> f(x).diff(x,2)
12*x**2

216
Однако этот прием не работает, когда выражение, возвращаемое функцией, содержит функции, вызываемые из модулей. Например, такое простое решение не работает, когда функция пользователя содержит функции модуля math. Дело в том, что подстановка символьных переменных в качестве аргументов функций модуля math приводит к ошибке. Но в модуле SymPy имеется функция lambdify(), которая умеет преобразовывать различные функции, в частности модуля math, в аналогичные функции модуля SymPy (или любого другого модуля, содержащего одноименные функции).
>>> from math import sin
>>> fn=lambda t: sin(t)**2
>>> fn(1)
0.7080734182735712
>>> from sympy import *
>>> x=symbols('x')
>>> g= lambdify(x, fn(x), "sympy")
>>> g(x)
sin(x)**2
>>> g(x).diff(x)
2*sin(x)*cos(x)
Таким образом, функция lambdify() позволяет преобразовывать функции
Python и символьные выражения друг в друга.
Иногда мы сталкиваемся с задачей преобразования строкового выражения в Python функцию. Для этого можно использовать функцию eval().
>>> eq = "x**2"
>>> func = lambda x: eval(eq)
>>> func(5)
25
А для преобразования строкового выражения в символьное подходит функция sympify().
>>> from sympy import *
>>> f=sympify('x**4'); f x**4
>>> x=S('x')
>>> f.diff(x), f.diff(x,x)
(4*x**3, 12*x**2)
Если символьное выражение содержит только знаки математических операций, то для его преобразования в строку можно использовать функцию str().
>>> f=x**3 -2*x**2+x
>>> str(f)
'x**3 - 2*x**2 + x'
Для выполнения кода, который содержится в строке, вы можете применить функцию eval(), но для этих целей имеется специальная функция exec().
>>> mc = 'print("hello world")'
>>> exec(mc)
hello world

217
Функцию eval() следует применять тогда, когда нужно получить значение.
>>> a=eval('2+2');a
4
Пример. Построить график символьного выражения одной символьной переменной, а также графики его первой и второй производных.
Первое решение состоит в использовании графической функции plot() пакета
SymPy.
>>> import sympy as sp
>>> x=sp.symbols('x')
>>> f0=x**4 - 15*x**2 - 10*x + 24
>>> f1=sp.diff(f0,x);f1 4*x**3 - 30*x - 10
>>> f2=sp.diff(f0,x,2); f2 6*(2*x**2 - 5)
>>> sp.plot(f0,f1,f2,(x, -5,5))
Поскольку управление графикой в графических функциях SymPy ограничено, то, возможно, вы пожелаете использовать графику других пакетов, а для этого потребуется преобразование символьного выражения.
Второе решение состоит в преобразовании символьного выражения и его производных в лямбда – функцию и построении ее графика функцией plot() модуля mpmath.
>>> import sympy as sp
>>> import mpmath as mp
>>> x=sp.symbols('x')
>>> f0=x**4 - 15*x**2 - 10*x + 24
>>> f1=sp.diff(f0,x)
>>> f2=sp.diff(f0,x,2)
>>> g0=sp.lambdify(x,f0,'mpmath')
>>> g1=sp.lambdify(x,f1,'mpmath')
>>> g2=sp.lambdify(x,f2,'mpmath')
>>> mp.plot([g0,g1,g2],[ -5,5])
Здесь мы создали символьное выражение f0= x**4-15*x**2-10*x+24 и вычислили его первую и вторую производные f1, f2. Затем символьные выражения, представляющие исходную функцию f0 и ее первую и вторую производные, преобразовали в лямбда – функции с помощью функции sympy.lambdify(), которой третьим аргументом указали желаемый тип

218 результата („mpmath‟). Для построения графиков использовали функцию plot() модуля mpmath.
Третье решение состоит в преобразовании символьного выражения и его производных в лямбда – функцию и построении ее графика функцией plot() модуля matplotlib.pyplot.
>>> import matplotlib.pyplot as plt
>>> import sympy as sp
>>> import numpy as np
>>> x=sp.symbols('x')
>>> f0=sp.sin(x)*sp.exp( -x**2)
>>> f1=sp.diff(f0,x)
>>> f2=sp.diff(f0,x,2)
>>> h0=sp.lambdify(x,f0,'numpy')
>>> h1=sp.lambdify(x,f1,'numpy')
>>> h2=sp.lambdify(x,f2,'numpy')
>>> t = np.linspace( -3, 3, 121)
>>> plt.plot(t,h0(t),'b',t,h1(t),'r',t,h2(t),'g', linewidth=3)
>>> plt.show()
Здесь, имея символьные выражения f0, f1 и f2, мы построили лямбда –
функции h0(x), h1(x), h2(x). Обратите внимание на третий аргумент функций lambdify(), которому присвоено значение 'numpy'. Это нужно для того, чтобы имена функций sin и exp, используемые в символьных выражениях f0, f1 и f2, были преобразованы в соответствующие функции модуля numpy. Затем лямбда – функции h0(x), h1(x), h2(x) мы используем для создания массивов значений, которые нужны для построения графиков функцией matplotlib.pyplot.plot().