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

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

6.1.2
Вычисление длин, площадей и объемов.
Длина дуги плоской кривой, заданной явным уравнением
 
x
y
y

в декартовых координатах, вычисляется по формуле
 


x
d
x
y
L
x
x



1 0
2
'
1
,
(1) где x
0
и x
1
определяют начальную и конечную точки дуги.
Пример. Длина цепной линии
x
ch
y

на отрезке
1 0


x
from scipy.integrate import quad import numpy as np import sympy as smp import matplotlib.pyplot as plt x=smp.symbols('x') fn=smp.cosh(x)
# символьное выражение уравнения f=smp.sqrt(1+fn.diff(x)**2)
# подынтегральное символьное выражение g=smp.lambdify(x, f, "numpy")
# преобразование в выражение numpy
L=quad(g,0,1)
# вычисление интеграла print(L)
# значение интеграла и точность
(1.1752011936438014, 1.3047354237503154e-14)
Сравним с точным значением
 
1
sinh print(np.sinh(1))

235 1.17520119364
# грфик цепной линии fig = plt.figure(facecolor='white') #,xlim=( -2,2), ylim=(0,4)) ax = plt.axes(xlim=( -2, 2), ylim=(0, 4)) fun=smp.lambdify(x, fn, "numpy") x=np.linspace( -2,2,81) ax.plot(x,fun(x),lw=3) ax.grid(True) plt.show()
Обратите внимание на инструкцию g=smp.lambdify(x,f,"numpy"), которая из символьного выражения f создает функцию g. Она потребовалась потому, что символьные выражения нельзя использовать в функциях, выполняющих численное интегрирование и, в частности, в функции quad.
■
Длина дуги плоской кривой, заданной в параметрическом виде
 
 
t
y
y
t
x
x


,
, вычисляется по формуле
 


 


t
d
t
y
t
x
L
t
t



1 0
2 2
'
'
(2) где t
0
и t
1
– значения параметра t в начальной и конечной точках дуги.
Пример. Длина астроиды
t
y
t
x
3 3
sin
,
cos


,

2 0


t
from scipy.integrate import quad import numpy as np import sympy as smp t=smp.symbols('t') xt=smp.cos(t)**3 yt=smp.sin(t)**3 dt=smp.sqrt(xt.diff(t)**2+yt.diff(t)**2)
# подынтегральное выражение g=smp.lambdify(t, dt, "numpy")
L=quad(g,0,2*np.pi) print(L)
(6.0, 6.661338147750939e-14)
В этом коде также существенной является инструкция g=smp.lambdify(t, dt, "numpy")

236
Она из символьного выражения dt относительно символьной переменной t создает функцию g, которую можно использовать в численном интегрировании.
■
Геометрический смысл определенного интеграла
 

b
a
x
d
x
f
состоит в том, что он представляет площадь под кривой
 
x
f
y



b
x
a


Пример. Найдем площадь фигуры, ограниченной сверху параболой
2 9
x
y


, а снизу осью абсцисс.
Для этого решим уравнение
0 12 2



x
x
и найдем точки, в которых парабола пересекается с осью абсцисс. import sympy as smp x=smp.symbols('x') y=-x**2+x+12 root=smp.solve(y, x) a=root[0]; b=root[1]; print(a,b)
-3 4
Теперь вычисляем площадь. s=smp.integrate(y, (x,a,b)) print('s=',s)
s= 343/6
В этом примере интеграл удалось вычислить символьно, однако это бывает нечасто.
■
Площадь поверхности, заданной явным уравнением
 
y
x
f
z
,

. Площадь куска поверхности, который проектируется на плоскость
xOy
в виде замкнутой области D, вычисляется по формуле




D
y
x
y
d
x
d
f
f
S
2 2
1
(3)
Пример. Площадь поверхности сферы. Уравнение
2 2
2
y
x
R
z



представляет полусферу радиуса R. Учитывая, что (3) дает площадь только верхней полусферы, мы приходим к следующей формуле для площади сферы:

















R
R
y
x
R
y
x
R
y
x
D
y
x
dy
z
z
dx
y
d
x
d
z
z
S
2 2
2 2
2 2
2 2
2 2
1 2
1 2
Имеем. from scipy.integrate import dblquad import numpy as np import sympy as smp x, y=smp.symbols('x y')
R=1

237 z=smp.sqrt(R**2-x**2-y**2) ds=smp.sqrt(1+z.diff(x)**2+z.diff(y)**2)
DS=smp.lambdify((y,x), ds, "numpy") g=lambda x: -np.sqrt(R**2-x**2) h=lambda x: np.sqrt(R**2 -x**2) s=2*dblquad(DS, -R,R,g,h)[0] print(s)
12.56637061435586
print(4*np.pi*R**2)
# точное значение площади
2 4
R

12.566370614359172
■
Площадь поверхности, получаемой вращением кривой
  

b
x
a
x
f
y



вокруг оси Ox, определяется по формуле
 
 


x
d
x
f
x
f
x
d
y
y
S
b
a
b
a
x






2 2
'
1 2
'
1 2


(4)
Пример. Вычислим площадь поверхности, образованной вращением вокруг оси Ox части кривой
x
y


4 2
, отсеченной прямой x=2.
Уравнение
x
y


4 2
определяет параболу с вершиной в точке (-4, 0) и осью симметрии Ox. Поэтому для вычисления площади поверхности вращения достаточно рассмотреть одну ветвь параболы
x
y


4
на отрезке [-4,2]
(точное значение площади равно
3
/
62

). from scipy.integrate import quad import numpy as np import sympy as smp x=smp.symbols('x') y=smp.sqrt(4+x) ds=y*smp.sqrt(1+y.diff(x)**2)
DS=smp.lambdify(x, ds, "numpy") s=2*np.pi*quad(DS,-4,2)[0] print(s)
64.92624817418037
■
Плоская кривая задана параметрическими уравнениями
 
 
t
y
y
t
x
x


,


T
t
t


0
. Площадь поверхности ее вращения вокруг оси
x вычисляется по формуле
 
 


 





T
t
t
d
t
y
t
x
t
y
S
0 2
2
'
'
2

(5)

238
Пример. Дана циклоида


t
t
a
x
sin


,


t
a
y
cos
1


. Найти площадь поверхности, образованной вращением одной дуги кривой вокруг оси Ox
(точное значение площади равно
3 64 2

a
). from scipy.integrate import quad import numpy as np import sympy as smp t=smp.symbols('x') a=1 x=a*(t-smp.sin(t)) y=a*(1-smp.cos(t)) ds=y*smp.sqrt(x.diff(t)**2+y.diff(t)**2)
DS=smp.lambdify(t, ds, "numpy") s=2*np.pi*quad(DS,0,2*np.pi)[0] print(s)
67.02064327658223
■
Объем тела U в декартовых координатах Oxyz вычисляется с помощью тройного интеграла


U
z
d
y
d
x
d
V
Пример. Объем конуса высотой H и радиусом основания R. Конус ограничен поверхностью


2 2
y
x
R
R
H
z



и плоскостью z=0. В декартовых координатах его объем выражается следующим трехкратным интегралом


















R
R
x
R
x
R
y
x
R
R
H
U
z
d
y
d
x
d
z
d
y
d
x
d
V
2 2
2 2
2 2
0
Вычислим его. from scipy.integrate import nquad import numpy as np f=lambda z,y,x: 1 zbnd=lambda y,x: [0, H/R*(R -np.sqrt(x**2+y**2))] ybnd=lambda x: [ -np.sqrt(R**2 -x**2),np.sqrt(R**2-x**2)] xbnd=lambda : [-R,R]
H=1; R=1
I=nquad(f,[zbnd,ybnd, xbnd]) print(I)
(1.047197551137596, 1.4896493976954173e-08)

239
Сравним полученный результат с известным точным значением
3 2
H
R
V


print(np.pi*R**2*H/3)
1.0471975511965976
6.2
Обыкновенные дифференциальные уравнения и системы.
Задача Коши для одного дифференциального уравнения n – го порядка состоит в нахождении функции, удовлетворяющей равенству
 




1
,...,
,
,



n
n
y
y
y
t
f
y
и начальным условиям
 
 
 
 
0 0
1 0
2 0
0 1
0
,...,
,
n
n
y
t
y
y
t
y
y
t
y





Перед решением эта задача (уравнение или система уравнений) должна быть переписана в виде системы ОДУ первого порядка






t
y
y
y
f
t
d
y
d
t
y
y
y
f
t
d
y
d
t
y
y
y
f
t
d
y
d
n
n
n
n
n
,
...,
,
,
,
...,
,
,
,
...,
,
,
2 1
2 1
2 2
2 1
1 1



(1) с начальными условиями
 
 
 
0 0
0 2
0 2
0 1
0 1
,
,
,
n
n
y
t
y
y
t
y
y
t
y



Модуль scipy.integrate имеет две функции ode() и odeint(), предназначенные для решения систем обыкновенных дифференциальных уравнений (ОДУ) первого порядка с начальными условиями в одной точке
(задача Коши). Функция ode() более универсальная, а функция odeint()
(ODE integrator) имеет боле простой интерфейс и хорошо решает большинство задач.
Функция odeint(). Функция odeint()имеет три обязательных аргумента и много опций. Она имеет следующий формат odeint(func, y0, t[,args=(), ...])
Аргумент func – это имя Python функции двух переменных, первой из которых является список y=[y1,y2,...,yn], а второй – имя независимой переменной.
Функция func должна возвращать список из n значений функций


t
y
y
f
n
i
,
,
,
1
при заданном значении независимого аргумента t. Фактически функция func(y,t) реализует вычисление правых частей системы (1).
Второй аргумент y0 функции odeint() является массивом (или списком) начальных значений
0 0
2 0
1
,
,
,
n
y
y
y
при t=t0. Третий аргумент является массивом моментов времени, в которые вы хотите получить решение задачи.
При этом первый элемент этого массива рассматривается как t0.

240
Функция odeint() возвращает массив размера len(t) x len(y0).
Функция odeint() имеет много опций, управляющих ее работой. Опции rtol (относительная погрешность) и atol (абсолютная погрешность) определяют погрешность вычислений
i
e
для каждого значения
i
y
по формуле
tol
i
tol
i
a
y
r
e



. Они могут быть векторами или скалярами. По умолчанию rtol=atol=1.49012e-8.
Договоримся, что в начале всех примеров этого параграфа выполняется следующее импортирование модулей и функций. import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt
Пример. Рассмотрим задачу Коши
01 0
)
5 1
(
,
5 3






y
y
t
y
Создадим функцию, реализующую правую часть уравнения. def dydt(y, t): return -5*t**3*y
Вызываем солвер и строим график. t = np.linspace(-1.5,1.5,31)
# вектор моментов времени y0 = 0.01
# начальное значение y = odeint(dydt, y0, t)
# решение уравнения y = np.array(y).flatten()
# преобразование массива plt.plot(t, y,'-sr',linewidth=3
)
# построение графика
Обратите внимание, что начальный момент времени определяется первым значением в массиве t, а не задается отдельно.
Пример. Решим дифференциальное уравнение
 
x
y
x
y



с начальным условием
 
1 0

y
. Для этой задачи мы знаем «точное» решение
 
x
e
x
x
y




2 1
, которое сможем сравнить с приближенным.
Создаем функцию, реализующую правую часть уравнения. def dydx(y, x):
# функция правой части
 
y
x
x
y



return x - y
Вызываем солвер и строим график. x = np.linspace(0,4,41)
# массив значений независимой переменной y0 = 1.0
# начальное значение y = odeint(dydx, y0, x)
# решение уравнения y = np.array(y).flatten()
# преобразование массива fig = plt.figure(facecolor='white')
# построение графика plt.plot(x, y,' -sb',linewidth=3)

241 ax = fig.gca() ax.grid(True)
Сравним графики «точного»
x
e
x
y




2 1
и приближенного решений. y_exact = x - 1 + 2*np.exp(-x) fig = plt.figure(facecolor='white') ax = fig.gca() ax.plot(x, y, x, y_exact, 'o'); ax.grid(True)
Теперь построим график модуля разности точного и приближенного решений. fig = plt.figure(facecolor='white') y_diff =np.abs(y_exact - y) plt.plot(x,y_diff, linewidt h=2) plt.ylabel("Error") plt.xlabel("x") plt.title("Error in numerical integration"); plt.grid(True)
Как видим, абсолютная погрешность не превосходит
8 10 5


Пример. Решим задачу Коши
1
)
0
(
'
,
0
)
0
(
,
0 5
1







z
z
z
z
z
Вначале приведем ее к задаче Коши для системы ОДУ 1-го порядка.
Обозначим
z
y
z
y



2 1
,
. Тогда приходим к следующей задаче
1
)
0
(
,
0
)
0
(
,
5 1
,
2 1
1 2
2 2
1






y
y
y
y
t
d
y
d
y
dt
y
d
Создаем функцию, реализующую вычисления правых частей полученной системы ОДУ. Ее первый аргумент должен быть списком (или массивом) имен

242 искомых функций, а второй – именем независимой переменной, даже если он не используется при формировании функции. def f(y, t): y1, y2 = y
# вводим имена искомых функций return [y2,-0.2*y2-y1]
Решаем систему ОДУ. t = np.linspace(0,20,41) y0 = [0, 1] w = odeint(f, y0, t) w array([[ 0. , 1. ],
[ 0.45623694, 0.79029901],
[ 0.76275767, 0.41642039],
[ 0.86239266, -0.01890162],
[ 0.07506831, 0.1135842 ],
[ 0.11799752, 0.05551654]])
Первый столбец представляет значения функции
 
t
y
1
, а второй – значения функции
 
t
y
2
в точках вектора t.
Выделяем вектора y1 и y2 из результирующей матрицы.
y1=w[:,0]
# вектор значений решения y2=w[:,1]
# вектор значений производной
Строим график решения (значения искомой функции находятся в первом столбце матрицы w, т.е. в векторе y1) fig = plt.figure(facecolor='white')
# следующий рисунок слева plt.plot(t,y1, ' -o',linewidth=2) plt.ylabel("z") plt.xlabel("t") plt.grid(True)
Строим график фазовой траектории. Последовательность решения та же, но количество искомых точек возьмем больше. t = np.linspace(0,20,201) y0 = [0, 1]
[y1,y2]=odeint(f, y0, t, full_output=False).T fig = plt.figure(facecolor='white')
# предыдущий рисунок справа plt.plot(y1,y2,linewidth=2) plt.grid(True)

243
Чтобы сразу получит вектора решений y1 и y2 командой
[y1,y2]=odeint(...).T , мы использовали опцию full_output=False. Если full_output=True, то функция odeint() возвращает дополнительную информацию в форме словаря (см. справочную систему), и тогда операция транспонирования не будет работать, точнее ее нужно применять к первому элементу результата.
Пример. Решим дифференциальное уравнение
2 1
t
y



на отрезке


100
;
a
с начальными условиями
 
 
a
a
y
ln

,
 
a
a
y
/
1


при a=0.001. Его точное решение
t
y
ln

Приведем задачу к системе ОДУ первого порядка и решим ее с относительной погрешностью rtol=0.5·10
-4
, 10
-6
, 10
-8
(абсолютная погрешность atol задана по умолчанию). Потом решим задачу с относительной погрешностью rtol, заданной по умолчанию, и различными значениями абсолютной погрешности atol=10
-2
,10
-3
,10
-4
Построим графики приближенных решений и график точного решения. import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt def f(y, x): y1, y2 = y return [y2,-1/x**2] a=0.001 y0=[log(a),1/a] x = np.linspace(a,100,200) z=np.log(x)
# точное решение fig = plt.figure(facecolor='white') plt.hold(True) reltol=[0.5e-4,1e-6,1e-8] for tol in reltol:
[y1,y2]=odeint(f, y0, x, f ull_output=False, rtol=tol).T plt.plot(x,y1,linewidth =2)
# следующий рисунок слева plt.plot(x,z,'k',linewidth=1)
# график точного решения.
plt.grid(True)
fig = plt.figure(facecolor='white') plt.hold(True) abstol=[1e-2,1e-3,1e-4] for tol in abstol:
[y1,y2]=odeint(f, y0, x, full_output=False, atol=tol).T plt.plot(x,y1,linewidth =2)
# следующий рисунок справа plt.plot(x,z,'k',linewidth=1)
# график точного решения.

244 plt.grid(True)
На левом рисунке тонкой черной линией (нижняя кривая) показан график точного решения. Другими цветами показаны графики приближенного решения при различной относительной погрешности. Как видим, относительной погрешности 0.5·10
-4
и 10
-6
для данной задачи явно недостаточно. Хорошее приближение к точному решению получилось только при rtol = 10
-8
(график этого решения сливается с кривой точного решения). На правом рисунке средняя кривая представляет точное решение. С ней сливается график приближенного решения, соответствующего atol=10
-4
Заметим, что названия опций абсолютная и относительная погрешность не полностью отвечают определениям этих терминов. Опции rtol и atol определяют погрешность вычислений
i
e
для каждого значения
i
y
по формуле
tol
i
tol
i
a
y
r
e



. Они будут соответствовать своим названиям, если в первой группе решений вы используете опцию atol=0, например, измените строку решения следующей
[y1,y2]=odeint(f,y0,x,full_output=False,rtol=tol,atol=0).T.
Аналогично, во второй группе строку решения можно заменить такой
[y1,y2]=odeint(f,y0,x,full_output=False,atol=tol,rtol=0).T
По умолчанию обе опции «точности» равны rtol=atol=1.49012e-8.
В этом примере мы использовали функцию hold(True/False), которая включает (или выключает) режим добавления новых графиков на текущую графическую область. Если функция hold() вызывается без аргумента, то текущий режим переключается на противоположный. Заметим, что у некоторых функций типа plot() имеется аналогично действующая опция hold=True/False. Обычно режим добавления графиков включен по умолчанию.