Введение в научный Python-1. Введение в научный Python
 Скачать 6.2 Mb.
|
|
Пример. Решим задачу Коши, описывающую движение тела, брошенного с начальной скоростью v 0 под углом к горизонту в предположении, что сопротивление воздуха пропорционально квадрату скорости. В векторной форме уравнение движения имеет вид g v v r m m , где ) (t r радиус – вектор движущегося тела, ) (t r v – вектор скорости тела, – коэффициент сопротивления, g m вектор силы веса тела массы m, g – ускорение свободного падения. Особенность этой задачи состоит в том, что движение заканчивается в заранее неизвестный момент времени, когда тело падает на землю. Если обозначить m k / , то в координатной форме мы имеем систему уравнений 2 2 y x x k x g y x y k y 2 2 к которой следует добавить начальные условия: 0 0 x , h y 0 (h начальная высота), cos 0 0 v x , sin 0 0 v y Положим y y y y x y x y 4 3 2 1 , , , . Тогда соответствующая система ОДУ 1 – го порядка примет вид g y y y k y y y y y y k y y y 2 4 2 2 4 4 4 3 2 4 2 2 2 2 2 1 Для модельной задачи положим h=0, k=0.01, g=9.81, v0=10, α=π/4. Создаем функции правой части системы и функцию «обработчик шага» ts = [ ] ys = [ ] def fout(t, y): # обработчик шага ts.append(t) ys.append(list(y.copy())) 259 def f(t, y): # функция правой части системы ОДУ k=0.01 g=9.81 y1, y2, y3, y4 = y return [y2, -k*y2*sqrt(y2**2+y4**2), y4, -k*y4*sqrt(y2**2+y4**2) -g] Решаем ОДУ и строим его график. tmax=1.41 # время движения, подбирается экспериментально alph=np.pi/4 # угол бросания тела v0=10.0 # начальная скорость ODE=ode(f) y0,t0=[0, v0*np.cos(alph), 0, v0*np.sin(alph)], 0 # начальные условия r=ODE.set_integrator('dopri5', max_step=0.1) # метод Рунге – Кутта r.set_solout(fout) # загрузка обработчика шага r=ODE.set_initial_value(y0, t0) # задание начальных значений ret = r.integrate(tmax) # решаем ОДУ Y=np.array(ys) fig, ax = plt.subplots() fig.set_facecolor('white') ax.plot(Y[:,0],Y[:,2],linewidth=3) # график решения ax.grid(True) Длительность полета tmax нами подбиралась «на глазок» из условия, что в конце полета вертикальная координата равна нулю. Также приходится «на глазок» определять максимальную высоту и дальность полета. Конечно, эти значения следует вычислять. Это можно реализовать в теле функции «обработчика шага» fout(). В ней нужно проверять, что вертикальная координата больше нуля, а иначе останавливать вычисления и запоминать значение момента времени. Также в обработчике можно запоминать дальность и высоту полета. В нашем примере функция fout()должна реагировать на обращение в ноль Y координаты тела (компоненты решения y3) при ее убывании (т.е. при 0 4 y y ), а также на обращения в ноль производной, соответствующей наивысшей точке траектории, т.е. на обращение компоненты решения y4 в ноль. Из физических соображений ясно, что производная обращается в ноль только в одной точке – в самой верхней точке траектории тела. Поэтому, неважно убывает или возрастает производная в этой точке. 260 Перепишем начальный код примера и функцию обработчика шага. ts = [ ] ys = [ ] FlightTime, Distance, Height =0,0,0 y4old=0 def fout(t, y): global FlightTime, Distance, Height,y4old ts.append(t) ys.append(list(y.copy())) y1, y2, y3, y4 = y if y4*y4old<=0: # достигнута точка максимума Height=y3 if y4<0 and y3<=0.0: # тело достигло поверхности FlightTime=t Distance=y1 return -1 y4old=y4 Здесь мы создали 3 глобальные переменные (FlightTime, Distance, Height ), предназначенные для хранения полетного времени, дальности и максимальной высоты полета. Также создали вспомогательную глобальную переменную y4old, которая содержит значение производной t y на предыдущем шаге вычислений. Поскольку производная t y вычисляется приближенно (вектор y4), то проверить равенство 0 y невозможно. Но можно проконтролировать изменение знака производной, проверив знак выражения y4*y4old, где y4 значение производной на текущем шаге, а y4old – значение производной на предыдущем шаге. Условие y4*y4old<=0 выполняется для наивысшей точки траектории. Поэтому, если y4*y4old<=0, то следует запомнить значение высоты y3. Второе условие y4<0 and y3<=0.0 обнаруживает первое нулевое или отрицательное значение Y координаты тела при убывании этой координаты (т.е. 0 4 y y ). Если условие выполнилось, то мы запоминаем дальность и время полета и останавливаем вычисления в методе integrate(). Напомним, что для остановки вычислений функция fout() должна вернуть -1. Функция правой части системы уравнений не меняется def f(t, y) : k=0.01 g=9.81 y1, y2, y3, y4 = y return [y2, -k*y2*sqrt(y2**2+y4**2), y4, -k*y4*sqrt(y2**2+y4**2) -g] Последующий код определяет максимально допустимый момент времени, решает систему ОДУ, печатает значение полетного времени, дальность и максимальную высоту. Затем строится траекторию движения тела. Обратите 261 внимание на инструкцию ODE.set_integrator('dopri5', max_step=0.05), которая использует опцию max_step. Если ее не задавать, то метод использует слишком большой шаг, и вы получите грубую траекторию, показанную на следующем рисунке слева. При задании опции max_step=0.05 количество точек будет больше и график, показанный на следующем рисунке справа, получается более гладким. tmax=100 # максимально допустимый момент времени alph=np.pi/4 # угол бросания тела v0=10.0 # начальная скорость ODE=ode(f) y0,t0=[0, v0*np.cos(alph), 0, v0*np.sin(alph)], 0 # начальные условия #r=ODE.set_integrator('dopri5') # след. график слева r=ODE.set_integrator('dopri5', max_step=0.05) # след. график справа r.set_solout(fout) r=ODE.set_initial_value(y0, t0) ret = r.integrate(tmax) print('Flight time = %.4f Distance = %.4f Height =%.4f '\ % (FlightTime,Distance,Height)) Y=np.array(ys) fig, ax = plt.subplots() fig.set_facecolor('white') ax.plot(Y[:,0],Y[:,2],' -or',linewidth=3) ax.grid(True) В результате вычислений мы получаем следующие значения времени движения, дальности и максимальной высоты. Flight time = 1.4157 Distance = 9.4739 Height =2.4443 Проведем в этой задаче исследование зависимости дальности полета от коэффициента сопротивления среды. Для этого мы должны переписать код функции правой части системы так, чтобы она содержала дополнительный аргумент k – коэффициент сопротивления среды. FlightTime, Distance, Height, y4old =0,0,0,0 def fout(t, y): # обработчик шага global FlightTime, Distance, Height,y4old ts.append(t) ys.append(list(y.copy())) y1, y2, y3, y4 = y if y4*y4old<=0: # достигнута точка максимума Height=y3 262 if y4<0 and y3<=0.0: # тело достигло поверхности FlightTime=t Distance=y1 return -1 y4old=y4 # функция правых частей системы ОДУ def f(t, y, k): # имеется дополнительный аргумент k g=9.81 y1, y2, y3, y4 = y return [y2, -k*y2*sqrt(y2**2+y4**2), y4, -k*y4*sqrt(y2**2+y4**2) -g] tmax=100 # максимально допустимый момент аремени alph=np.pi/4 # угол бросания тела v0=10.0 # начальная скорость K=[0.1,0.2,0.3,0.5] # анализируемые коэффициенты сопротивления y0,t0=[0, v0*np.cos(alph), 0, v0*np.sin(alph)], 0 # начальные условия ODE=ode(f) ODE.set_integrator('dopri5', max_step=0.01) ODE.set_solout(fout) fig, ax = plt.subplots() fig.set_facecolor('white') for k in K: # перебор значений коэффициента сопртивления ts, ys = [ ],[ ] ODE.set_initial_value(y0, t0) # задание начальных значений ODE.set_f_params(k) # передача дополнительного аргумента k # в функцию f(t,y,k) правых частей системы ОДУ ODE.integrate(tmax) # решение ОДУ print('Flight time = %.4f Distance = %.4f Height =%.4f ' \ % (FlightTime,Distance,Height)) Y=np.array(ys) ax.plot(Y[:,0],Y[:,2],linewidth=3,label='k=%.1f'% k) ax.grid(True) ax.set_xlim(0,8) ax.set_ylim(-0.1,2) Здесь мы использовали метод ODE.set_f_params(k), который функции правой части f(t, y, k) передает дополнительный параметр k. Обратите внимание на то, что в цикле для построения новой траектории мы каждый раз обновляем начальные условия. Графики траекторий показаны на следующем рисунке. 263 Параметры траекторий печатаются в окне исполнительной системы. Flight time = 1.2316 Distance = 5.9829 Height =1.8542 Flight time = 1.1016 Distance = 4.3830 Height =1.5088 Flight time = 1.0197 Distance = 3.5265 Height =1.2912 Flight time = 0.9068 Distance = 2.5842 Height =1.0240 Пример. Упругий мяч имеет начальное положение и скорость. Сопротивление воздуха пренебрежимо мало, но энергия движения расходуется при отскоке мяча от земли. Пусть при отскоке от земли вертикальная скорость мяча составляет 90% от вертикальной скорости в момент падения. Смоделируем такое движение. Уравнение движения (без учета сопротивления воздуха) имеет вид g r m m , где g – вектор ускорения свободного падения, имеющий направление вертикально вниз. В покоординатной форме уравнение принимает вид g y x , 0 Задачу следует дополнить начальными условиями y x v y v x h y x x 0 0 0 0 , 0 , 0 , 0 . Система распадается на два независимых уравнения, первое из которых имеет решение t v x t x x 0 0 . Это указывает на то, что в горизонтальном направлении движение происходит с постоянной скоростью x v 0 . В нашем примере положим 0 0 x . Второе уравнение также интегрируется, но мы хотим показать, как его можно решить численно, используя функцию ode(). В момент отскока t k вертикальная скорость y меняет знак и абсолютное значение, т.е. k отскока до k отскока после t y t y 9 0 ) ( Функция правой части системы ОДУ 1 – го порядка имеет вид: def f(t, y): g=9.81 y3, y4 = y return [y4,-g] Функция обработчика шага имеет вид: def fout(t, y): ts.append(t) # запоминаем моменты времени ys.append(list(y.copy())) # и значения y и vy y3, y4 = y if y4<0 and y3<=0.0: # при достижении поверхности return -1 # остановить интегрирование Создаем функцию горизонтальной координаты тела. horX=lambda t: vox*t # t v x t x x 0 , где 0 0 x 264 Вводим исходные данные. ts,ys = [ ],[ ] tmax=15 # максимально допустимый момент аремени vox=1 # начальная скорость вдоль оси x voy=10 # начальная скорость вдоль оси y Создаем ODE объект. ODE=ode(f) ODE.set_integrator('dopri5', max_step=0.01, nsteps=500) ODE.set_solout(fout) Y0, t0 =[0,voy], 0 # начальные значения (начальная высота=0) Осталось задать ODE объекту начальные условия и решить систему ОДУ. При достижении мячом земли функция fout() останавливает интегрирование, предварительно запомнив все моменты времени и координаты решения в глобальных списках ts и ys. Чтобы строить решение в последующие моменты времени, мы должны возобновить вычисления с новыми начальными условиями, соответствующими моменту отскока. Последние элементы в списках ts и ys (время и скорость) используются при задании новых начальных значений для следующего гладкого участка траектории мяча. При этом вертикальная скорость мяча в момент отскока становится положительной и составляет 90% от абсолютного значения вертикальной скорости в момент падения. На всех участках используются общие списки ts и ys. После завершения цикла все последовательные значения моментов времени находятся в списке ts, а соответствующие им значения вертикальной координаты и скорости – в списке ys. for i in range(6): # моделируем 6 отскоков ODE.set_initial_value(Y0, t0) ODE.integrate(tmax) t0=ts[-1] # новый начальный момент времени Y0=[0,0.9*abs(ys[ -1][1])] # новые начальные значения Создаем массивы горизонтальных X и вертикальных Y[:,0] координат мяча, и используем их для построения траектории мяча. t=np.array(ts) Y=np.array(ys) X=horX(t) fig, ax = plt.subplots() fig.set_facecolor('white') ax.plot(X,Y[:,0],' -',linewidth=3) ax.grid(True) 265 График траектории движения построен для вектора начальной скорости 10 , 1 , 0 0 0 y x v v v |