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

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

Пример. Исследуем решение задачи Коши
t
x
x
sin
3



,
0
)
0
(
,
0
)
0
(



x
x
Преобразуем уравнение 2 – го порядка в систему уравнений 1 – го порядка.
Сделав замену
 
   
 
t
y
t
x
t
y
t
x
2 1
,



, приходим к задаче
 
0
)
0
(
,
0
)
0
(
,
sin
,
2 1
3 2
2 1






y
y
t
t
y
t
d
y
d
y
dt
y
d
l
Создаем функцию правых частей системы. def f(y, t): y1, y2 = y
# вводим имена искомых функций return [y2,-y1**3+sin(t)]

245
Решаем систему ОДУ. t = np.linspace(0,50,201) y0 = [0, 0]
[y1,y2]=odeint(f, y0, t, full_output=False).T
Строим график решения. fig = plt.figure(facecolor='white') plt.plot(t,y1, ' -o',linewidth=2)
# график решения plt.grid(True)
Строим фазовую траекторию.
fig = plt.figure(facecolor ='white') plt.plot(y1,y2,linewidth=2)
# следующий рисунок слева plt.grid(True)
Как видите, количества точек, выбранных нами, недостаточно. Вы можете в команде t = np.linspace(0,50,201) увеличить это количество, например, до 1000. Можно поступить иначе – выполнить интерполяцию решения и использовать бо̀ льшее количество точек только для построения графика фазовой траектории. from scipy.interpolate import interp1d
F1=interp1d(t, y1,kind='cubic')
F2=interp1d(t, y2,kind='cubic') tnew=np.linspace(0,50,1001) fig = plt.figure(facecolor='white') plt.plot(F1(tnew),F2(tnew),linewidth=2)
# следующий рисунок справа plt.grid(True)
Поясним, как здесь выполнялась интерполяция.
Функция interp1d(xvec,yvec[,kind='linear',...]) принимает одномерные массивы xvec, yvec данных и возвращает ссылку на интерполяционную функцию. Тип интерполяции задается опцией kind, которая может принимать значения 'linear', 'nearest', 'zero', 'slinear', 'quadratic,

246
'cubic', а также быть целым положительным числом (порядок интерполяционного сплайна). Вот простой пример. from scipy.interpolate import interp1d x = np.linspace(0, 12, num=13, endpoint=True) y = np.sin(x)**2
# массив значений f = interp1d(x, y,kind='cubic')
# кубическая сплайн интерполяция xx=np.linspace(0, 12, num=49, endpoint=True)
# новый набор точек fig = plt.figure(facecolor='white') plt.plot(x, y, ' -or', xx, f(xx),' --b',linewidth=3)
# ломаная и сплайн
Пример. Исследуем решение задачи Коши для системы уравнений
)
(
sin
)
(
01 0
),
(
t
x
t
y
y
t
y
x






,
1 2
)
0
(
,
0
)
0
(


y
x
Весь код решения соберем в одну функцию fun() с подфункцией f(y,t), вычисляющей правую часть системы уравнений. def fun(): def f(y, t): y1, y2 = y return [y2,-0.01*y2-sin(y1)] t = np.linspace(0,100,501) y0 = [0, 2.1]
[y1,y2]=odeint(f, y0, t, full_output=False).T fig = plt.figure(facecolor='white') plt.plot(t,y1, linewidth =2)
# график решения x(t) слева plt.grid(True) fig = plt.figure(facecolor='white') plt.plot(y1,y2,linewidth=2)
# фазовая траектория справа plt.grid(True)
Для решения задачи выполните следующие команды import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt fun()
# вызов основной функции
На левом рисунке показан график решения x(t), а на правом – фазовая траектория.

247
■
Пример. Исследуем двухвидовую модель «хищник – жертва», впервые построенную Вольтерра для объяснения колебаний рыбных уловов. Имеются два биологических вида, численностью в момент времени t соответственно x(t) и y(t). Особи первого вида являются пищей для особей второго вида
(хищников). Численности популяций в начальный момент времени известны.
Требуется определить численность видов в произвольный момент времени.
Математической моделью задачи является система дифференциальных уравнений Лотки – Вольтерра

















y
x
d
c
dt
dy
x
y
b
a
dt
dx
где a, b, c, d – положительные константы. Проведем расчет численности популяций при
1
,
1
,
3



d
c
a
для трех значений параметра
2
,
3
,
4

b
Начальные значения положим
1
)
0
(
,
2
)
0
(


y
x
import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt def f(y, t,params): y1, y2 = y a,b,c,d=params return [y1*(a -b*y2),y2*(-c+d*y1)] t = np.linspace(0,7,71) y0 = [2, 1] fig = plt.figure(facecolor='white') plt.hold(True) for b in range(4,1, -1): params=[3,b,1,1] st='a=%d b=%d c=%d d=%d' % tuple(params)
[y1,y2]=odeint(f, y0, t,args=(params,), full_output=False).T plt.plot(y1,y2,linewidth=2, label=st) plt.legend(fontsize=12) plt.grid(True) plt.xlim(0,3)

248 plt.ylim(0,3)
Из графиков фазовых траекторий видно, что численность популяций меняется периодически.
Теперь будем менять начальные условия, оставляя параметры a=3,b=3,c=1,d=1 неизменными. def f(y, t,params): y1, y2 = y a,b,c,d=params return [y1*(a -b*y2),y2*(-c+d*y1)] t = np.linspace(0,7,71) ic = np.linspace(1.0, 3.0, 11)
# начальные значения для 1-й функции fig = plt.figure(facecolor='white') for r in ic: y0 = [r, 1.0] y = odeint(f, y0, t,args=([3,3,1,1],)) plt.plot(y[:,0],y[:,1],linewidth=2) plt.grid(True)
# следующий рисунок слева
Добавьте к предыдущему коду следующие инструкции. x1 = np.linspace(0, 3, 31) y1 = np.linspace(0, 2, 21)
X1,Y1 = np.meshgrid(x1, y1)
# создание сетки
DX1,DY1 = f([X1, Y1],t,[3,3,1,1])
LN=np.sqrt(DX1**2+DY1**2)
LN[LN == 0] = 1.
# исключение деления на 0
DX1 /= LN
# нормировка
DY1 /= LN plt.quiver(X1,Y1, DX1,DY1, LN, pivot='mid',cmap=plt.cm.jet)
Результатом работы примера будет график, показанный на предыдущем рисунке справа.

249
Пример. Решим систему дифференциальных уравнений


z
y
x
t
d
z
d
x
z
x
y
t
d
y
d
y
t
d
x
d











1 0
,
1 1
0
,
3
с начальными условиями
0
)
0
(
,
1
)
0
(
,
1
)
0
(



z
y
x
и построим ее фазовый портрет. Отличие от предыдущих задач состоит в том, что фазовая кривая расположена в трехмерном пространстве. import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D
def f(y, t): y1, y2, y3 = y return [y2,0.1*y2-y1*(y3-1)-y1**3,y1*y2-0.1*y3] y0=[1,1,0] t = np.linspace(0,25,501) fig = plt.figure(facecolor='white') ax=Axes3D(fig)
[y1,y2,y3]=odeint(f, y0, t, full_output=False).T ax.plot(y1,y2,y3,linewidth=3)
Фазовая траектория показана на следующем рисунке.
Пример. Исследуем поведения математического маятника. Пусть масса груза равна единице, а стержень, на котором подвешена масса, невесом. Тогда дифференциальное уравнение движения груза имеет вид
0
sin
2









k
где
 
t

угол отклонения маятника от положения равновесия (нижнее положение), параметр k характеризует величину трения,
l
g /
2


(g ускорение свободного падения, l – длина маятника). Для определения конкретного движения к уравнению надо добавить начальные условия
 
0 0



,
 
0 0





Преобразуем уравнение к системе ОДУ 1 – го порядка. Если обозначить





v
u
,
, то получим следующую задачу
 









)
sin
2
u
v
k
v
v
u

,
 
0 0
0
,
)
0
(





v
u
Пусть k=0.5,
10 2


, и зададим следующие начальные значения
5
,
0 0
0 0




v


import numpy as np

250 from scipy.integrate import odeint import matplotlib.pyplot as plt
Создаем функцию pend=lambda y,t: [y[1], -0.5*y[1]-10*np.sin(y[0])]
Решаем систему и строим график (следующий рисунок слева) v0=5
# начальная скорость
0 0
v



t = np.linspace(0,20,401)
[y1,y2]=odeint(pend, [0, v0], t, full_output=False).T fig = plt.figure(facecolor='white') plt.plot(t,y1,linewidth=2) plt.grid(True
)
Строим фазовую траекторию (рисунок справа) fig = plt.figure(facecolor='white') plt.plot(y1,y2,linewidth=2) plt.grid(True)
Как видно из левого графика максимальный угол отклонения маятника не превышают
2
/

и колебания маятника затухают.
Увеличим начальную скорость до 10. Для этого в предыдущем коде заменим инструкцию v0=5 на v0=10. Решаем задачу и строим график решения
(следующий рисунок слева) и фазовую траекторию (следующий рисунок справа).
Максимальное значение угла составляет примерно 9 радиан. Маятник сделал один полный оборот вокруг точки закрепления (угол отклонения увеличился на

2
), а затем колебания затухают в окрестности значения

2 радиан (для маятника угол поворота

2 представляет то же, что и 0 радиан, т.е. положение равновесия).
Построим несколько графиков угла отклонения (следующий рисунок слева) и фазовых траекторий (следующий рисунок справа), задавая различную начальную скорость v0=5,7,9,11.

251 fig = plt.figure(facecolor='white') ax1 = fig.add_subplot(1,2,1) ax1.grid(True) ax1.hold(True) ax2 = fig.add_subplot(1,2,2) ax2.grid(True) ax2.hold(True) t = np.linspace(0,20,401) for v0 in range(5,12,2):
[y1,y2]=odeint(pend, [0,v0], t, full_output=False).T ax1.plot(t,y1,linewidth=2) ax2.plot(y1,y2,linewidth=2)
Как видим, начальная скорость v0=5 и 7 недостаточна, чтобы маятник прошел верхнюю точку и сделал хотя бы один полный оборот. При начальной скорости v0=9 маятник совершает один полный оборот, а затем его колебания затухаю. При v0=11 маятник выполнил два полных оборота и только после этого его колебания стали затухать вокруг положения равновесия.
Пример. Решим задачу движения планеты вокруг Солнца под действием тяготения. Она записывается в виде системы ОДУ 2-го порядка




2
/
3 2
2 2
/
3 2
2
,
y
x
y
k
y
y
x
x
k
x










, где
   
t
y
t
x
,
- координаты движущейся планеты. Преобразуем исходную систему к системе ОДУ 1-го порядка. Для этого введем обозначения
y
z
y
z
x
z
x
z






4 3
2 1
,
,
,
. Тогда


























2
/
3 2
3 2
1 3
4 4
3 2
/
3 2
3 2
1 1
2 2
1
z
z
z
k
z
z
z
z
z
z
k
z
z
z
Для модельной задачи выберем k=1 и зададим начальные условия
 
 
 
 
1 0
,
0 0
,
0 0
,
1 0
4 3
2 1




z
z
z
z

252
Код решения ОДУ и построения графика траектории планеты приведен ниже. import numpy as np from scipy.integrate import odeint, ode import matplotlib.pyplot as plt from matplotlib.patches import Circle def f(y, t): y1, y2, y3, y4 = y return [y2,
-y1/(y1**2+y3**2)**(3/2), y4,
-y3/(y1**2+y3**2)**(3/2)] t = np.linspace(0,20,1001) y0 = [1, 0, 0, 0.4]
[y1,y2, y3, y4]=odeint(f, y0, t, full_output=False).T fig, ax = plt.subplots() fig.set_facecolor('white') ax.plot(y1,y3,linewidth=1) circle = Circle((0, 0), 0.03, facecolor='orange')
# круг ax.add_patch(circle) plt.axis('equal') plt.grid(True)
Теперь добавим код создания анимации движения планеты. import matplotlib.animation as animation fig2 = plt.figure(facecolor='white') ax = plt.axes(xlim=( -0.2, 1.2), ylim=(-0.4, 0.4) ) line, = ax.plot([ ], [ ], lw=3) circle = Circle((0, 0), 0.05, facecolor='orange') ax.add_patch(circle) pc = plt.Circle((y0[0], y0[2]), 0.02, fc='b') ax.add_patch(pc) ax.grid(True) def redraw(i): x = y1[0:i+1] y = y3[0:i+1] line.set_data(x, y) pc.center=(x[ -1],y[-1])

253 anim =animation.FuncAnimation(fig2,redraw,frames=126,interval=50) plt.show()
При запуске программы открывается графическое окно со статическим рисунком, показанным выше. Закройте его мышкой. После этого откроется второе окно с анимацией. Несколько кадров приведено на следующем рисунке.
Анимация создается с помощью функции FuncAnimation в графическом окне fig2 путем многократно вызыва функции redraw. Значение опции frames=126 подбиралось экспериментально.
Пример. Решим систему уравнений Лоренца























3 2
1 3
2 3
1 2
1 2
1
y
b
y
y
y
y
y
r
y
y
y
y
s
y
, где s, r, b некоторые параметры. Положим
3
,
25
,
10



b
r
s
и выберем следующие начальные условия:
 
,
1 0

x
 
 
10 0
,
1 0



z
y
import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D
Создаем функцию правой части системы уравнений. s,r,b=10,25,3 def f(y, t): y1, y2, y3 = y return [s*(y2 -y1),
-y2+(r-y3)*y1,
-b*y3+y1*y2]
Решаем систему ОДУ и строим ее фазовую траекторию, которая показана на следующем рисунке слева. t = np.linspace(0,20,2001) y0 = [1, -1, 10]
[y1,y2,y3]=odeint(f, y0, t, full_output=False).T fig = plt.figure(facecolor='white')
# следующий рисунок слева ax=Axes3D(fig) ax.plot(y1,y2,y3,linewidth=2)

254
Система описывает хаотический аттрактор, поэтому любые малые изменения в начальных условиях будут приводить к существенному изменению решения.
Вы это сможете заметить по меняющемуся количеству петель фазовой траектории. Например, скорректируйте начальное значение первой функции на
0.0001, т.е. вместо инструкции y0=[1, -1, 10] введите следующую строку y0=[1.0001,-1,10]
. Выполните код и получите график фазовой траектории, показанный на предыдущем рисунке справа.
Функция ode(). Вторая функция модуля scipy.integrate, которая предназначена для решения дифференциальных уравнений и систем, называется ode(). Она создает объект ОДУ (тип scipy.integrate._ode.ode).
Имея ссылку на такой объект, для решения дифференциальных уравнений следует использовать его методы. Аналогично функции odeint(), функция ode(func) предполагает приведение задачи к системе дифференциальных уравнений вида (1) и использовании ее функции правых частей. Отличие только в том, что функция правых частей func(t,y) первым аргументом принимает независимую переменную, а вторым – список значений искомых функций. Например, следующая последовательность инструкций создает объект ODE, представляющий задачу Коши. from scipy.integrate import ode f = lambda t, y: 2*t # dy/dt=2*t
ODE=ode(f)
# ссылка на объект ОДУ
Используя методы объекта ODE, можно выбрать численный метод решения.
ODE.set_integrator('vode')
Здесь указано название метода „vode‟. Он реализует неявный метод Адамса
(для нежестких задач), а для жестких задач – метод, основанный на формулах обратного дифференцирования (backward differentiation formulas, BDF). Список доступных методов можно узнать на странице справки функции ode().
С помощью метода set_initial_value() задают начальные значения.
Формат его вызова следующий:
ODE.set_initial_value(y,t=0.0).
Например,
OD.set_initial_value(1, 0) # y(0)=1

255
Здесь сказано, что начальное значение t0 независимой переменной равно нулю
(t0=0) , а начальное значение искомой функции y0=1. Если решается система, то в качестве начальных значений передается список.
Вычисления выполняются методом integrate(tnew). Например, rez=ODE.integrate(0.1)
В данном случае возвращаемое значение rez представляет одноэлементный массив, содержащий значение искомой функции в момент t=tnew. При этом используются заданные ранее начальные условия. Результат вычисления в момент t=tnew также находится в атрибуте ODE.y.
Соберем описанные инструкции в единый пример, в котором решается задача
 
1 0
,
2


y
t
t
d
y
d
(задача имеет решение
1 2


t
y
), и выполним код. from scipy.integrate import ode f = lambda t, y: 2*t
# dy/dt=2*t (функция правой части)
ODE=ode(f)
# ссылка на объект ОДУ
ODE.set_integrator('vode')
ODE.set_initial_value(1, 0)
# y(0)=1
ODE.integrate(0.1)
# вычисление решения в момент 0.1
print(ODE.y)
[ 1.0100005]
Вы можете задавать новые значения независимой переменной и получать значения искомой функции.
ODE.integrate(0.2) print(OD.y)
[ 1.0400005] print(ODE.integrate(0.3))
[ 1.0900005]
Атрибут ODE.t хранит момент времени, для которого вычислено значение решения.
ODE.t
0.5
# точка последнего вычисления
Можно построить таблицу решений. Для примера рассмотрим следующую систему 3–х дифференциальных уравнений
,
2
,
,
y
x
t
d
z
d
x
z
t
d
y
d
z
y
t
d
x
d







с начальными условиями
 
 
 
2 0
,
0 0
,
1 0



z
y
x
f=lambda t,y: [y[1] -y[2],y[2]-y[0],y[0]-2*y[1]]
OD=ode(f) y0,t0=[1,0,2], 0 r=OD.set_integrator('vode') r=OD.set_initial_value(y0, t0) t1=1.0 dt=0.1

256 while OD.successful() and OD.t <=t1: print('%.2f \t%.3f \t%.3f \t%.3f' % (OD.t, OD.y[0],OD.y[1],OD.y[2]))
OD.integrate(OD.t+dt)
0.00 1.000 0.000 2.000 0.10 0.801 0.114 2.079 0.20 0.609 0.254 2.113 0.30 0.432 0.413 2.099 0.40 0.274 0.585 2.034 0.50 0.144 0.762 1.920 0.60 0.044 0.938 1.759 0.70 -0.020 1.103 1.556 0.80 -0.046 1.250 1.316 0.90 -0.033 1.373 1.049 1.00 0.019 1.465 0.764
Здесь в условии цикла while мы использовали метод ODE.successful(), который возвращает True, если шаг интегрирования ДУ завершился успешно.
Многократно решать задачу Коши для того, чтобы получить таблицу значений решения, неэффективно. У объекта ODE имеется метод
ODE.set_solout(fout)
, который загружает в объект специальную функцию fout()
(обработчик шага). Она вызывается после каждого шага интегрирования. В теле этой функции можно запоминать значение текущего момента времени и соответствующие ему значения решения. Если функция вернет -1, то процесс интегрирования в методе integrate() завершится на текущем шаге. Чтобы процесс интегрирования не завершился, функция fout() ничего не должна возвращать (None) или возвращать 0.
Построим график решения последней задачи. ts = [ ] ys = [ ] def fout(t, y):
# функция «обработчик шага» ничего не возвращает ts.append(t)
# запоминание t ys.append(list(y.copy()))
# запоминание y f=lambda t,y: [y[1]-y[2],y[2]-y[0],y[0]-2*y[1]]
# функция правой части
OD=ode(f)
# создание объекта ОДУ
y0, t0=[3.5,1.5, 2], -4 r=OD.set_integrator('dopri5')
# метод Рунге – Кутты 4(5) порядка
# с контролем шага r.set_solout(fout)
# загрузка обработчика шага r=OD.set_initial_value(y0, t0)
# задание начальных значений ret = r.integrate(6.0)
# решение ОДУ до момента t=6
Y=np.array(ys) fig = plt.figure(facecolor ='white')
# график фазовой траектории ax=Axes3D(fig) ax.plot(Y[:,0],Y[:,1],Y[:,2],linewidth=3)

257
В нашем примере функция обработчика шага fout() добавляет в глобальные списки ts и ys значение независимого аргумента t и список трех чисел – значений решения при этом t. Затем двойной список ys мы преобразуем в двумерный массив Y, компоненты которого представляют вектора значений решения задачи.
Замечание
. Не каждый метод позволяет использование функции обработчика шага. Сейчас эту возможность поддерживают методы Рунге – Кутты 'dopri5' и 'dop853'.
Каждый из методов имеет свои собственные опции настройки.
Например, измените в предыдущем коде в инструкции r=OD.set_integrator('dopri5') имя метода на 'dop853' и выполните пример. Вы получите грубый график фазовой траектории, показанный на следующем рисунке. Чтобы построить гладкую кривую следует использовать опцию max_step, которая задает максимальное значение шага, используемое солвером. Замените инструкцию определения метода следующей командой r=OD.set_integrator('dop853', max_step=0.1)
, и снова выполните пример. В результате вы получите гладкую фазовую кривую, такую же, какая показана на предыдущем рисунке.
■
У функции ode(f[,jac=None]) имеется опция явного задания якобиана
j
i
j
i
y
f
J



правой части системы (1). Обычно функция вызывается в формате ode(f)
, а якобиан вычисляется численно. Вы можете уменьшить погрешность вычислений, если явно создадите эту функцию.

258
У методов функции ode(), кроме опции max_step , имеются другие опции, например, first_step, min_step, nsteps (максимальное количество внутренних шагов) и много других опций. Познакомьтесь с ними самостоятельно по справочной системе.