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

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

Пример. Нарисовать в одном графическом окне символьный график и график matplotlib
, не прибегая к преобразованиям выражений.
%matplotlib qt import sympy as sp import numpy as np import matplotlib.pyplot as plt x=sp.symbols('x') p1=sp.plot(x**3-x,(x,-1.5,1.5))
# sympy график fig=plt.gcf() fig.set_facecolor('white') ax=fig.gca()
# получение осей sympy графика

219 t=np.linspace( -1.5,1.5,16) z=t**3-t ax.plot(t,-z,'-sr', linewidth=3)
# matplotlib график ax.grid(True) ax.axis('equal')
Вначале мы построили график символьного выражения x**3-x. Затем получили ссылку на объект axes этого графика. Этот объект принадлежит модулю matplotlib.pyplot и имеет метод plot() этой библиотеки.
Пример. Построить график функции
x
y
cos

и графики ее сумм Тейлора различного порядка.
%matplotlib qt import sympy as sp import numpy as np import matplotlib.pyplot as plt x = symbols('x') f = sp.cos(x) fvec = np.vectorize(lambda val: f.subs(x, val).evalf()) z = np.linspace( -8,8, num=300)
fig = plt.figure(facecolor='white') plt.plot(z, fvec(z), linewidth=5,color='r',label='cos(x)') ax = fig.gca() for order in range(3, 20, 2): s = f.series(x, 0, order) svec = np.vectorize(lambda val: s.removeO().subs(x, val)) ax.plot(z, svec(z), ' -', label='n=%i' % order, linewidth=3)
ax.set_xlim(-8, 11) ax.set_ylim(-2, 2) ax.legend(fontsize=12) ax.grid(True)
Пример. График функции и ее касательной.
Уравнение касательной к кривой
 
x
f
y

в точке
0
x
имеет вид
 
 

0 0
0
'
x
x
x
f
x
f
y



. Чтобы его создать нужно использовать символьное дифференцирование, и правую часть уравнения задавать как символьное выражение. Для построения графика будем использовать функцию matplotlib.pyplot.plot
, которая не умеет работать с символьными

220 выражениями. Поэтому требуется преобразование символьных выражений, представляющих уравнения кривой и ее касательной, к виду, пригодному для функции plot. Ниже приведен код, учитывающий эту особенность. import numpy as np import matplotlib.pyplot as plt from sympy import * x = symbols('x') x0=1.25
# точка касания y=-exp(-x)*sin(pi*x)
# символьное уравнение кривой y1=y.diff(x) y0=y.subs(x,x0) y10=y1.subs(x,x0)
# значение производной в точке func=lambdify(x, y, "numpy")
# уравнение кривой tangent=lambdify(x, y0+y10*(x-x0), "numpy")
# уравнение касательной dlt=0.5
# полуширина (по оси x) отрезка касательной
Xt=np.linspace(x0 -dlt,x0+dlt,21)
# если производная ноль, то символьное уравнение
# касательной вырождается в константу (число)
if y10==0:
Yt=y0*np.ones(np.size(Xt)) else:
Yt=tangent(Xt)
Xf=np.linspace(-1,2,51)
Yf=func(Xf) plt.plot(Xf,Yf,linewidth=3,c='r') plt.plot(Xt,Yt, linewidth=3,c='b') plt.ylim([1.1*Yf.min(),1.1*Yf.max()]) plt.grid(True) plt.show()
Отметим еще раз, что существенным моментом этого кода является использование функции lambdify для преобразования двух символьных выражений в функции func и tangent модуля numpy. func=lambdify(x, y, "numpy") tangent=lambdify(x, y0+y10*(x-x0), "numpy")
Другим ньюансом нашего кода является проверка условия if y10==0 ....
Если производная y10 в точке x0 обращается вноль, то инструкция tangent=lambdify(x, y0+y10*(x-x0), "numpy") создаст не функцию tangent(x)
, а константу. Тогда инструкция Yt=tangent(Xt) создаст не

221 вектор значений Yt, а одно число, и график не будет нарисован (вы получите сообщение об ошибке).
Пример. Касательная и нормаль к эллипсу, заданному неявным уравнением
0 1
2 2
2 2



b
y
a
x
Если кривая задана неявным уравнением
 
0
,

y
x
F
и точка


0 0
, y
x
принадлежит кривой, то уравнение касательной к кривой в этой точке имеет вид






0
,
'
,
'
0 0
0 0
0 0




y
y
y
x
F
x
x
y
x
F
y
x
,
(1) а уравнение нормали






0
,
'
,
'
0 0
0 0
0 0




x
x
y
x
F
y
y
y
x
F
y
x
(2)
В следующем коде в начале создается символьное выражение
 
y
x
F ,
. Абсцисса точки касания x0 задается, а ордината y0 определяется путем решения уравнения


0
,
0

y
x
F
. Учитывая, что функция solve может вернуть список из одного или двух корней, для y0 мы выбираем последнее значение, поскольку оно есть всегда. Значения


0 0
, y
x
F
x

и


0 0
, y
x
F
y

находятся символьным дифференцированием с последующей подстановкой координат точки касания


0 0
, y
x
. После этого, с помощью функции lambdify мы создаем numpy функцию
 
y
x
F ,
и функции правых частей уравнений (1) и (2), которые и используются при построении графика. import numpy as np import matplotlib.pyplot as plt from sympy import * x,y = symbols('x y') a=2; b=1
F=x**2/a**2+y**2/b**2 -1
Fx=F.diff(x)
Fy=F.diff(y) x0=1.5
# x координата точки касания y0=solve(F.subs(x,x0),y)[ -1]
# находим y координаты точек касания
# и выбираем последнюю из списка
Fx0=Fx.subs([(x,x0),(y,y0)])
# значение производной Fx в точке касания
Fy0=Fy.subs([(x,x0),(y,y0)])
# значение производной Fy в точке касания
# создание numpy функций из символьных выражений
Fun=lambdify((x,y), F, "numpy")
Tan=lambdify((x,y), Fx0*(x -x0)+Fy0*(y-y0), "numpy")
Norm=lambdify((x,y), Fx0*(y -y0)-Fy0*(x-x0), "numpy")
# построение графиков эллипса, касательной и нормали t=np.linspace( -a,a,31)
Xf,Yf=np.m eshgrid(t,t) plt.figure() plt.contour(Xf, Yf, Fun(Xf,Yf), [0] , linewidths=3,colors='b' ) tx=np.linspace(x0 -1.0,x0+1.0,21)

222 y0=float(y0)
# преобразование символьного значения в числовое ty=np.linspace(y0 -1,y0+1,21)
Xt,Yt=np.meshgrid(tx,ty) plt.contour(Xt, Yt, Tan(Xt,Yt), [0] , linewidths=2,colors='r' ) plt.contour(Xt, Yt, Norm(Xt,Yt), [0] , linewidths=2,colors='g' )
#plt.gcf().gca().axis([ -3,3,-3,3]); plt.gcf().gca().axis('equal'); plt.gcf().set_facecolor('white') plt.show()
Полученное изображение эллипса, касательной и нормали к нему, приведено на следующем рисунке.
6.
Научные вычисления с пакетом SciPy
Пакет SciPy является библиотекой математических процедур. Многие из них являются Python оболочками, обеспечивающими доступ к библиотекам, написанным на других языках программирования. Поскольку функции этих библиотек хранятся в откомпилированном виде, то эти процедуры обычно работают очень быстро. В этой главе мы переходим к описанию возможностей модулей SciPy.
6.1
Численное интегрирование.
6.1
.1 Вычисление интегралов.
Многие прикладные задачи сводятся к определению интегралов – одинарных, двойных, тройных или многократных. Однако аналитически вычислить удается немногие из них. Тогда интеграл можно найти численно. В этом параграфе мы рассмотрим некоторые из функций пакета scipy.integrate, предназначенные для численного интегрирования.
Приведем список наиболее часто используемых функций интегрирования из модуля scipy.integrate.
Функция
Описание
quad однократное численное интегрирование dblquad вычисление двойного интеграла tplquad вычисление тройного интеграла nquad вычисление n – кратного интеграла cumtrapz вычисление первообразной
Вычислим определенный интеграл



1 1
sin
x
d
x
e
x

223
Первым шагом является создание функции, вычисляющей подынтегральное выражение def f( x ): return np.exp(-x)*np.sin(x)
Функция f(x) должна содержать комбинацию стандартных математических функций, имена которых имеются в модуле numpy и scipy (или их подмодулях, например, в scipy.special). Она также может содержать вызовы функций пользователя.
Для вычисления интеграла предназначена функция quad. Первым аргументом ей нужно передать ссылку на функцию f (имя функции), а вторым и третьим — нижний (a) и верхний (b) пределы интегрирования.
I=quad(f, a, b)
Таким образом, вычислить интеграл можно следующим образом: from scipy.integrate import quad import numpy as np def f( x ): return np.exp(-x)*np.sin(x)
I=quad(f, -1, 1) print(I
)
(-0.6634936666312413, 1.2783578503385453e-14)
Обычно (как в этом примере) функция quad возвращает кортеж, содержащий значение интеграла и абсолютную погрешность. Однако при задании опции full_output=True функция quad будет возвращать словарь с дополнительной информацией.
По умолчанию функция quad вычисляет приближенное значение интеграла с абсолютной и относительной погрешностями epsabs=1.49e-08, epsrel=1.49e-08
. Здесь epsabs и epsrel имена опций, задающих эти погрешности.
Функция f(x,...) может принимать много аргументов, но интегрирование будет выполняться только по первому. Для передачи подынтегральной функции дополнительных аргументов используется опция args
, которой следует передавать кортеж с дополнительными аргументами. f = lambda x,y : x**2+y**2 y,err = quad(f, 0, 3, args=(1,))
# вычисление




3 0
2 2
x
d
y
x
при y=1 print(y)
12.0
Вычислим интеграл





1 0
2
x
d
c
x
b
x
a
с параметрами
1
,
2
,
3



c
b
a
. Это можно сделать, используя следующий код: from scipy.integrate import quad def f(x, a, b,c): return a*x**2+b*x+c
I = quad(f, 0, 1, args=(3,2,1))

224 print(I)
(3.0, 3.3306690738754696e-14)
Иногда удобно не передавать вспомогательные аргументы через опцию args, а использовать лямбда – оболочку вокруг подынтегральной функции. Вот, например, как можно вычислить предыдущий интеграл. def f(x, a, b,c): return a*x**2+b*x+c
I = quad(lambda x: f(x,3,2,1), 0, 1) print(I)
(3.0, 3.3306690738754696e-14)
Можно численно находить интеграл на бесконечном участке вещественной оси.
Это значит, что верхний и/или нижний пределы интегрирования могут быть бесконечными. Для обозначения пределов ±∞ используется идентификатор numpy.inf
. Например, вычислим интеграл



0 2
x
d
e
x
i1 = quad((lambda x: np.exp( -x**2)), 0, np.inf) print(i1)
(0.8862269254527579, 7.101318390915439e-09)
Этот интеграл вычисляется аналитически. Сравним значения. x=smp.symbols('x') s1=smp.integrate(smp.exp( -x**2),(x,0,smp.oo))
#
2 0
2





x
d
e
x
print(s1) print(s1.evalf())
sqrt(pi)/2 0.886226925452758
Вычислим еще один несобственный интеграл символьно и численно. from scipy.integrate import quad, dblquad, tplquad import numpy as np i1 = quad(lambda x: 1/(1+x**2), 0, np.inf)
#



0 2
1 1
x
d
x
print(i1)
(1.5707963267948966, 2.5777915205989877e-10)
Сравним с точным значением x=smp.symbols('x') s1=smp.integrate(1/(1+x**2),(x,0,smp.oo))
#
2 1
1 0
2





x
d
x
print(s1)
pi/2
Подынтегральная функция сама может содержать интеграл. Например, вычислим интеграл


 

1 0 0 2
2
x
d
y
d
y
x
x
. Это можно сделать, используя следующий код:

225 def f(y,u):
# аргумент y указываем первым return y**2+u**2 def g(x): return quad(f,0,x,args=(x,))[0]
# возвращаем только значение интеграла
I=quad(g,0,1) print('I=',I)
I=(0.3333333333333333, 3.700743415417188e-15)
Проверим результат, используя символьное интегрирование. x,y=smp.symbols('x y') si=smp.integrate(x**2+y**2,(y,0,x),(x,0,1)) print(si)
1/3
Однако такие интегралы проще вычислять, если использовать функцию dblquad
. Она применяется для вычисления повторных интегралов вида
 
 
 
 
b
a
x
h
x
g
x
d
y
d
y
x
f
,
. Функция dblquad имеет следующий синтаксис: scipy.integrate.dblquad(ffun, a, b, gfun, hfun[,...])
Здесь ffun является именем подынтегральной функции (двух переменных).
При этом имя внутренней переменной интегрирования y при создании функции ffun должно быть указано первым. Это значит, что инструкция определения функции должна иметь вид: def ffun(y,x):
# порядок аргументов важен
Аргументы a и b задают значения нижнего и верхнего пределов интегрирования по переменной x (пределы внешнего интеграла), gfun и hfun являются именами функций, определяющих нижний и верхний пределы интегрирования по переменной y (пределы внутреннего интеграла). Имеются также опции args, epsabs, epsrel, значение которых такое же, как и для функции quad.
Вот, например, как следует вычислять предыдущий интеграл. Вначале определим функции f, g и h. f=lambda y,x: x**2+y**2 g=lambda x: 0 h=lambda x: x
Обратите внимание на то, что даже если функции g и h нижнего и верхнего пределов интегрирования являются константами (как g(x) в нашем случае), их следует задавать функциями. Затем используем функцию dblquad.
I=dblquad(f, 0, 1, g, h)
#




1 0 0 2
2
x
d
y
d
y
x
x
print(I)
(0.3333333333333333, 3.700743415417188e-15)
Как и для однократного интеграла, функция dblquad возвращает два числа: значение интеграла и абсолютную погрешность (absolute uncertainty).

226
Вычислим еще один повторный интеграл
 




1 1
1 0
2 2
2 1
x
d
y
d
y
x
x
f=lambda y,x: np.sqrt(1-x**2-y**2) g=lambda x: 0 h=lambda x: np.sqrt(1 -x**2)
I=dblquad(f, -1, 1, g, h)
#
 




1 1
1 0
2 2
2 1
x
d
y
d
y
x
x
print(I)
(1.0471975511965979, 1.162622832669544e-14)
Этот интеграл представляет объем четверти шара (интеграл по полукругу от функции
2 2
1
y
x
z



) и нам известно его точное значение
3
/

. Сравним это значение с полученным приближенным значением интеграла. print(np.pi/3)
1.0471975511965976
Как видим, значения практически совпадают.
Функцию dblquad можно использовать для вычисления двойных интегралов, предварительно «обнулив» функцию вне области интегрирования.
Однако этот подход имеет свои недостатки.
Например, чтобы проинтегрировать функцию
 
2 2
1
,
y
x
y
x
f



по области единичного круга можно проинтегрировать функцию
 










1
,
0 1
,
1
,
2 2
2 2
2 2
y
x
y
x
y
x
y
x
F
по прямоугольнику
1 1
,
1 1






y
x
from scip y.integrate import dblquad import numpy as np def F(y,x): if x**2+y**2<1: return 1-x**2-y**2 else: return 0 g=lambda x: -1 h=lambda x: 1
I = dblquad( F, -1,1,g,h, epsabs=1e-6, epsrel=1e -6) print(I)
(1.5707950459207884, 1.2636992237532934e-06)
Сравним с точным значением
2
/

print(np.pi/2)
1.5707963267948966
Обратите внимание на то, что абсолютная и относительная погрешность нами задана меньшей, чем значения по умолчанию. Использование значений по умолчанию приводит к сообщению о медленной сходимости интеграла и получению неправильного результата (проверьте). Задание значений

227 epsabs=1e-7, epsrel=1e-7 сообщения не вызывает, но результат
(1.3121922832006239, 1.5164007161061477e-08) все еще неверный.
Ниже показан график «обрезанной» функции
 
y
x
F ,
Тем не менее, поверхность интегрирования может содержать ребра. Вот пример вычисления объема правильной четырехугольной пирамиды со стороной a и высотой h







2 3
1
a
h
V
. Ее изображение при
1

h
и
2

a
(ломаную, проходящую по ребрам пирамиды) можно построить следующим кодом. fig = plt.figure() ax=Axes3D(fig) x=[-1,0,0,-1,0,0,1,0,1,0] y=[0,0,-1,0,1,0,0,1,0,-1] z=[0,1,0,0,0,1,0,0,0,0] ax.plot(x,y,z,linewidth= 3) plt.show()
Уравнение поверхности пирамиды имеет вид
y
x
z



1
, а область интегрирования ограничена снизу и сверху кривыми
 
1


x
x
g
и
 
x
x
h


1
. Тогда ее объем можно найти следующим образом. def f(y,x): return 1 -np.abs(x)-np.abs(y) g=lambda x: np.abs(x) -1 h=lambda x: 1-np.abs(x)
I = dblquad(f, -1,1,g,h) print('I=',I[0])
I= 0.6666666666666667
Трехкратный интеграл


 
 
 
 



b
a
x
h
x
g
y
x
r
y
x
q
z
d
x
y
z
f
y
d
x
d
,
,
,
,
вычисляется с помощью функции tplquad, которая имеет следующий синтаксис: tplquad(func,a,b,gfun,hfun,qfun,rfun[,args=(), epsabs=1.49e-08, epsrel=1.49e-08])
Аргумент func представляет подынтегральную функцию, по крайней мере, трех переменных, задаваемых в порядке z, y, x. Аргументы a,b определяют

228 пределы интегрирования по x; gfun,hfun задают нижний и верхний пределы интегрирования по y и являются функциями одной переменной x; qfun,rfun
- определяют нижний и верхний пределы интегрирования по z и являются функциями двух переменных x,y (в указанном порядке). Смысл опций args, epsabs, epsrel такой же, как и для функции quad.
Вычислим трехкратный интеграл


 





2 1
2 x
x
y
x
y
x
z
d
z
y
x
y
d
x
d
from scipy.integrate import tplquad def f(z, y, x):
# Обратите внимание на порядок задания аргументов return x+y+z
I=tplquad(f, 1, 2, lambda x: x, lambda x: 2*x, lambda x, y: x - y, lambda x, y: x + y) print(I)
(40.0, 4.440892098500626e-13)
Вычислим интеграл







T
z
y
x
z
d
y
d
x
d
I
3 1
, распространенный на тетраэдр T,
(трехмерный симплекс) ограниченный плоскостями
0
,
0
,
0



z
y
x
и
1



z
y
x
. Тетраэдр Т показан на следующем рисунке.
Вначале преобразуем интеграл по области в трехкратный. Имеем.




 













1 0
1 0
1 0
3 3
1 1
x
y
x
T
z
y
x
z
d
y
d
x
d
z
y
x
z
d
y
d
x
d
I
Тогда код вычисления интеграла, стоящего в правой части может быть следующим. def f(z, y, x): return 1/(1+x+y+z)**3
I=tplquad(f, 0, 1, lambda x: 0, lambda x: 1-x, lambda x, y: 0, lambda x, y: 1-x-y) print(I)
(0.03407359027997266, 3.782928446046958e-16)
Для рассмотренного интеграла известно точное значение


5 2
ln
8 16 1


Сравним его с полученным приближенным значением.
Is=(8*np.log(2)-5)/16

229 print(Is)
0.03407359028
Вычислим объем шара единичного радиуса. Для этого вычислим объем части шара O, расположенной в первом октанте


0
,
0
,
0



z
y
x
, а затем умножим результат на 8. При этом удобно использовать сферическую систему координат.








1 0
2
/
0 2
/
0 2
2
sin
8
sin
8 8












d
d
d
d
d
d
dxdydz
V
O
O
Тогда имеем. def dv(theta,phi,ro): return ro**2*np.sin(theta)
V=8*tplquad(dv,0,1, lambda phi: 0, lambda phi: np.pi/2, lambda ro, phi: 0, lambda ro,phi: np.pi/2)[0] print(V)
4.1887902047863905
Сравним с известным значением объема единичного шара

3 4
print(4/3*np.pi)
4.1887902047863905
Для вычисления n – кратных интегралов в модуле scipy.integrate имеется функция nquad. Она имеет следующий синтаксис: nquad(func, ranges, args=None, opts=None)
Здесь func имя функции, интеграл от которой вычисляется. Аргумент ranges является последовательностью (iterable object), каждый элемент которой должен быть списком (или кортежем) из двух выражений (или двух чисел).
Вместо пары выражений можно использовать функцию, которая возвращает список из двух значений. Пара значений в ranges[0] задает пределы интегрирования по x0, ranges[1] – по x1, и так далее. Например для четырехкратного интеграла, порядок следования аргументов функций должен быть таким, как написано ниже.










 
 




3 3
3 2
3 2
3 2
1 0
3 2
1 0
3 2
1 3
2 1
,
,
,
,
0 3
2 1
0
,
,
1 2
3
,
,
,
h
g
x
h
x
g
x
x
x
h
x
x
x
g
x
x
h
x
x
g
x
d
x
x
x
x
f
x
d
x
d
x
d
Т.о., при создании подынтегральной функции f аргументы следует указывать в порядке от x0 – переменной внутреннего интегрирования, до xn – переменной внешнего интегрирования. def f(x0,x1,x2,x3):
Аргументы функций, задающих пределы интегрирования, тоже нужно указывать в приведенном порядке. Например, порядок аргументов функции g0 должен быть следующим: def g0(x1,x2,x3):

230
Границы интегрирования задаются списком пар ranges. Пары состоят из нижнего и верхнего пределов интегрирования, и перечисляются от внутреннего интеграла к внешнему. Каждая пара может быть списком из двух констант, списком из двух функций или функцией, возвращающей пару выражений.
Например, вычисление указанного четырехкратного интеграла может быть выполнено следующей инструкцией: nquad(f, [[g0,h0],[g1,h1],[g2,h2],[g3,h3]])
Функция func, кроме переменных интегрирования, может принимать дополнительные аргументы, которые передаются ей через опцию args функции nquad.
Вычислим интеграл
 
 

0 1 5
x
d
y
d
y
e
y
x
from scipy.integrate import nquad import numpy as np def f(y, x): return np.exp(-x*y) / y**5
I=nquad(f, [[1, np.inf],[0, np.inf]]) print(I)
(0.2000000000189363, 1.3682975855986121e-08)
Обратите внимание на то, что порядок аргументов функции f должен соответствовать порядку, в котором передаются пределы интегрирования: от внутренней переменной к внешней. Внутреннему интегралу по аргументу y здесь соответствуют пределы [1,np.inf], а внешнему (по переменной x) – соответствуют пределы [0,np.inf]. Если хотя бы один из пределов интегрирования непостоянен, то оба элемента пары должны быть функциями.
Вычислим интеграл






!
2
!
!
!
1 1
1 0
1 0
1 0
,
0









 






r
q
p
r
q
p
y
d
y
x
y
x
x
d
y
d
x
d
y
x
y
x
x
r
q
p
y
x
y
x
r
q
p
Имеем. from scipy.integrate import nquad import numpy as np def f(x, y, p, q, r): return x**p*y**q*(1-x-y)**r def ybounds(x): return [0, 1 -x] def xbounds(): return [0, 1]
I=nquad(lambda y,x: f(x,y,1,2,3), [ybounds, xbounds]) print(I)
(0.0002976190476190476, 1.039774842432043e-17) print(12/np.math.factorial(8))
# точное значение
0.00029761904761904765

231
Обратите внимание на порядок аргументов лямбда – функции. Он соответствует порядку передачи пределов интегрирования – от внутреннего интеграла к внешнему.
Рассмотрим еще один двойной интеграл
96 1
2
/
1 0
2 1
0

 

x
d
y
d
y
x
x
(его значение известно). Имеем. f=lambda y, x: x*y xbnd=lambda : [0, 0.5] ybnd=lambda x: [0, 1-2*x]]
I=nquad(f, [ ybnd, xbnd]) print(I)
(0.010416666666666668, 4.101620128472366e-16)
Заметим, что в этом примере инструкцию nquad можно использовать в виде
I=nquad(f,[ybnd,[0,0.5]])
, и не создавать функцию xbnd.
С помощью функции nquad вычислим тройной интеграл.


















1 0
1 0
1 0
2 2
2 1
0
,
9
,
0 2
2 2
2 2
2 2
2 2
4 9
4 9
x
y
x
z
y
x
z
y
x
z
d
z
y
x
z
y
x
y
d
x
d
z
d
y
d
x
d
z
y
x
z
y
x
Имеем. def f(z, y, x): return x*y*z/np.sqrt(9*x**2+4*y**2+z**2) def zbnd(y,x): return [0, np.sqrt(1 -x**2-y**2)] def ybnd(x): return [0,np.sqrt(1 -x**2)] def xbnd(): return [0,1]
I=nquad(f,[zbnd,ybnd, xbnd]) print(I)
(0.012222222222280263, 9.273154430624787e-09)
Известно, что объем n – мерного симплекса
0
,
,
0
,
0 2
1



n
x
x
x
,
1 2
1




n
x
x
x
равен
!
1
n
(изображения тетраэдра, являющегося трехмерным симплексом, показано ранее). Определим численно объем четырехмерного симплекса
0
,
0
,
0
,
0 4
3 2
1




x
x
x
x
,
1 4
3 2
1




x
x
x
x
. Для этого вычислим следующий интеграл:











1 0
1 0
1 0
4 1
0 3
2 1
1 3
2 1
2 1
x
x
x
x
x
x
x
d
x
d
x
d
x
d
V
Имеем. from scipy.integrate import nquad import numpy as np def f(x4, x3, x2,x1): return 1

232 def lim4(x3,x2,x1): return [0, 1 -x1-x2-x3] def lim3(x2,x1): return [0, 1 -x1-x2] def lim2(x1): return [0, 1 -x1] def lim1(): return [0, 1]
V=nquad(f,[lim4,lim3,lim2,lim1]) print(V)
(0.041666666666666664, 1.1030064411555582e-14)
Сравним со значением
!
4
/
1
print(1/np.math.factorial(4))
0.041666666666666664
До сих пор мы вычисляли интеграл на фиксированном отрезке [a,b], но часто нам нужно вычислять интеграл
 

x
a
t
d
t
f
с переменным верхним пределом, который с точностью до константы совпадает с первообразной
 

x
d
x
f
. В пакете scipy.integrate есть функция cumtrapz, которая вычисляет подобный интеграл. Она имеет следующий синтаксис: scipy.integrate.cumtrapz(y, x=None, dx=1.0, axis=-1,initial=None)
Здесь y – одномерный массив значений интегрируемой функции, x – одномерный массив абсцисс, dx – расстояние между абсциссами (используется только, когда x=None), axis – задает ось интегрирования (по умолчанию используется значение –1, т.е. последняя ось). Параметр initial задает число, используемое в качестве первого значения возвращаемого массива.
Обычно оно должно быть равно нулю. По умолчанию ему присваивается значение None, что означает, что значение x[0] не используется, и возвращаемый массив на единицу короче, чем массив значений функции y.
Если массив x не передается (массив y вы должны передавать всегда), то аргумент dx определяет расстояние (по оси x) между элементами массива y.
Обратите внимание, что результат не зависит от начальной абсциссы x[0], и она не требуется.
С помощью функции cumtrapz построим график функции

x
t
d
t
0 2
cos
. Значение этого интеграла известно
4 2
sin
2
x
x

, и мы сравним его с численным решением. Имеем. from scipy.integrate import cumtrapz import matplotlib.pyplot as plt import numpy as np x = np.linspace(0, 6, num=25)

233 y=np.cos(x)**2 yint = cumtrapz(y, x, initial=0)
#

x
t
d
t
0 2
cos exact=lambda x: x/2+np.sin(2*x)/4
#
4 2
sin
2
/
x
x

[p1,p2]=plt.plot( x, exact(x), 'b -', x, yint, 'ro' ) p1.set_linewidth(3) p2.set_markersize(8) fig=plt.gcf() fig.set_facecolor('white') plt.show()
■
С помощью функции cumtrapz построим график функции


x
t
t
d
e
0 2
. Этот интеграл не выражается через элементарные функции, и с точностью до множителя совпадает с интегралом ошибок
 
x
erf
А именно,
 
x
t
d
e
x
t
erf
2 0
2




from scipy.integrate import cumtrapz import matplotlib.pyplot as plt import numpy as np import scipy.special as sp x = np.linspace(0, 3, num=16) y=np.exp(-x**2) yint = cumtrapz(y, x, initial=0) exact=lambda x: np.sqrt(np.pi)/2*sp.erf(x)
[p1,p2]=plt.plot( x, exact(x), 'b -', x, yint, 'ro' ) p1.set_linewidth(3) p2.set_markersize(8) fig=plt.gcf() fig.set_facecolor('white') plt.show()
Построим график функции
 


t
d
t
x
F
x


0
cos
10
cos
. Для этого интеграла не существует представления ни через элементарные, ни через специальные функции. Тем не менее, график этой функции можно построить. Имеем.

234 from scipy.integrate import cumtrapz import matplotlib.pyplot as plt import numpy as np x1 = np.linspace(0, 6, num=61) y1=np.cos(10*np.cos(x1)) yint1 = cumtrapz(y1, x1, initial=0) x2 = np.linspace(0, -6, num=61) y2=np.cos(10*np.cos(x2)) yint2 = cumtrapz(y2, x2, initial=0)
# объединяем массивы
X=np.append(x2[-1:0:-1],x1)
Y=np.append(yint2[-1:0:-1],yint1) p1=plt.plot( X, Y, 'b -', linewidth=3 ) fig=plt.gcf() fig.set_facecolor('white') fig.gca().grid(True) plt.show()