CC BY
116
Электронный журнал Cloud of Science. 2015. T. 2. № 2
http:/ / cloudofscience.ru ISSN 2409-031X
Движения планет: расчет и визуализация в среде Mathcad или Часы Кеплера
В. Ф. Очков1, Е. П. Богомолова1, Д. А. Иванов1, К. Писачич2
1Национальный исследовательский университет «Московский энергетический институт», 111250, Москва, ул. Красноказарменная, 14
University North, Mechanical engineering, Varazdin 104. brigade 3 42000 Varazdin, Croatia
e-mail: [email protected], [email protected]
Аннотация. В статье показано, как можно использовать математический пакет Mathcad для решения задачи о движении планет (материальных точек) под действием гравитационных сил на плоскости и в трехмерном пространстве. Задача проецируется на небесную механику. Обсуждаются вопросы пределов точности численных методов решения обыкновенных дифференциальных уравнений. Водится понятие «Часов Кеплера» с неравномерным движением стрелок, иллюстрирующих второй закон Кеплера.
Ключевые слова: Mathcad, движение планет, Часы Кеплера.
Не только специализированные пакеты, такие, например, какие описаны в [1], но и универсальные математические программы позволяют решать задачи о движении материальных точек с заданной массой (планет, спутников, комет и т. д.) согласно второму закону Ньютона m ■ a = 2 F и закону всемирного тяготения F = G ■ m ■ m2 / r2.
Так, например, в среде математической программы Mathcad есть функция Odesolve, численно решающая системы обыкновенных дифференциальных уравнений (ОДУ) [2], к которым сводится задача о движении материальных точек. Если использование этой функции дополнить средствами анимации [3], то можно получить довольно интересные и поучительные решения, связанные с темой статьи.
Рассмотрим несколько примеров.
Пример 1 (тестовый). Великолепная восьмерка.
Движение трех материальных точек по орбите в виде знака бесконечности (или «великолепной восьмерки» по терминологии [1]).
В. Ф. Очков, Е. П. Богомолова, Д. А. Иванов, К. Писачич
В работе [4] дано аналитическое решение этого частного случая движения трех материальных точек1. Мы же сейчас решим эту задачу в среде Mathcad численно и сравним результаты численного и аналитического решений.
Для того чтобы в среде Mathcad 15 численно решить систему ОДУ, нужно ввести в расчет ключевое слово Given (Дано), записать саму систему уравнений, используя апостроф (он вводится аккордом Ctrl+F7) для обозначения первой производной или двойной апостроф для второй производной2, а также ввести начальные условия — положение материальных точек и их скорости в начальный момент времени (на старте t = 0). В нашем случае дифференциальные уравнения движения материальных точек записываются парами — по горизонтальному (х) и вертикальному (у) направлениям.
На рис. 1 можно видеть систему из шести ОДУ, записанную для трех материальных точек (черная, красная и синяя3): произведение массы на ускорение (левая часть уравнения) равно сумме двух гравитационных сил, действующих на рассматриваемую материальную точку от ее двух «соседок». Вернее, не самих ускорений и сил, а их проекций в горизонтальном (х) и вертикальном (у) направлениях4. Каждое уравнение можно разделить на массу соответствующей планеты, оставив слева только вторую производную. Кроме того, можно убрать переменную G (гравитационная постоянная), так как в наших выкладках она принимается (пока — см. окончание статьи) за единицу. Это упрощает и ускоряет расчет. Но мы намеренно оставляем эти переменные в уравнениях, чтобы максимально сохранить физику задачи [5]: слева произведение массы точки на ее ускорение, справа силы, действующие на точку — планету или спутник.
1 Для двух точек аналитические решения есть для всех случаев, сводящихся к эллипсу (окружности), параболе и гиперболе. Для трех и более точек аналитические решения есть для ограниченного числа случаев.
2 Тут можно применять и оператор взятия производной d /dt.
3 Среде Mathcad 15 есть возможность давать переменным разные цвета, чем мы и воспользовались в расчете вместо того, чтобы приписывать переменным индексы: 1, 2, 3 и т. д. С другой стороны, статья может быть опубликована и в черно-белом издании, где цвет переменных пропадет. Да и цветное издание может читать дальтоник, что тоже будет не очень хорошо. Когда-то давно первый автор послал статью о цвете в программах в один журнал. Статью не опубликовали — рецензент сказал, что цветные дисплеи и принтеры еще не скоро появятся и появятся ли вообще. А сам он (рецензент) — дальтоник. Через пять лет после попытки публикации этой статьи появился язык Visual Pascal, где цвет использовался для выделения программистских конструкций: встроенных, пользовательских операторов, ремарок, сообщений об ошибках и т. д.
4 Некоторые читатели, увидев куб, а не квадрат в знаменателе дробей на рис. 1 и не видя числитель, могут тут вспомнить старый анекдот: «Если б на голову Ньютона упало не яблоко, а более весомый предмет (кокос, например), то в законе всемирного тяготения расстояние возводилось бы не во вторую, а в третью степень».
Рисунок 1. Запись системы ОДУ в среде Mathcad 15
Блок решения системы ОДУ в среде Mathcad 15 заканчивается вызовом функции Odesolve5 (рис. 2), имеющей следующие аргументы:
- вектор имен функций, являющихся решением системы ОДУ (в среде Mathcad Prime за именами функции нужно в скобках указать ее аргумент — см. ниже);
- переменная, по которой ведется численное интегрирование (аргумент искомых функций; в среде Mathcad Prime эту переменную не указывают);
- конечное (правое) значение интервала интегрирования («левое» значение — это нуль);
- число точек, по которому ведется табулирование искомых функций (необязательный параметр — по умолчанию имеем 1000 точек).
Кроме того, нажатием правой кнопки мыши можно выбрать для функции Odesolve один из четырех возможных методов численного решения задачи:
- Adams/BDF — смешанный алгоритм Адамса и обратной формулы дифференцирования;
- Fixed (интегрирование методом Рунге-Кутты с фиксированным шагом — на сайте http://twtmas.mpei.ac.ru/mas/worksheets/Euler.mcd можно увидеть
5 В среде Mathcad есть и другие инструменты численного решения дифференциальных уравнений (обыкновенных и в частных производных), но в данной статье они не рассматриваются.
В. Ф. Очков, Е. П. Богомолова, Д. А. Иванов, К. Писачич
сравнение этого метода с методом Эйлера; см. также конец статьи, где эти методы описываются и сравниваются);
- Adaptive (интегрирование методом Рунге-Кутты с переменным шагом — см. http://twtmas.mpei.ac.ru/mas/worksheets/rkadapt.mcd);
- Radau — алгоритм RADAUS для жестких систем ОДУ.
Рисунок 2. Настройка функции Odesolve на нужный метод решения систем обыкновенных дифференциальных уравнений
Если трем материальным точкам задать начальные условия, взятые из [4] и показанные на рис. 3, то решением системы ОДУ (см. рис. 1), будут функции, графическое отображение которых (след-в-след) показано на рис. 4а.
m х0 VD х'о У'о"1 m х0 УО х'о УО m х0 Уо х'о y'üj
1 0.97000436 -0.24308763
1 -0.97000436 0.24308753
0.93240737 0.86473146 ^
2 2 0.93240737 0.86473146
2 2 о 0 о -0.93240737 -0.86473146.,
Рисунок 3. Начальные условия для задачи о «великолепной восьмерке» при движении трех планет
Сравнение «телеметрических данных, снятых с орбиты» и показанных на рис. 4а, с данными из [4] дает полное совпадение результатов, что может свидетельствовать об адекватности нашей модели, зафиксированной в уравнениях рис. 1. Но если быть точными, то происходит следующее. На первых примерно ста «витках» орбиты не наблюдается отклонения численного решения от аналитического. Но при увеличении значения t происходит некое «размывание» орбиты вследствие накапливаемой ошибки численного решения задачи (рис. 46). Эта особенность хо-
рошо просматривается на втором примере, где в модель (см. рис. 1) добавлены новые уравнения и новые слагаемые у существующих для дополнительных двух точек — для коричневой и зеленой.
¿3 РГау Апгтзйоп
....................\ ,,--..............^ V J ■ 1 1 1
1 = 3.267 у(1) = 0.573 = 1 269 = 0 715 Уа!егу ОсЬУш + МаШсай 15
и
а)
Ь)
Рисунок 4. Кадры анимации https://www.ptcusercommunitv.сот/у1йго8/5777: а) движение трех материальных точек (планет) по орбите в форме знака бесконечности (показаны векторы скорости, который у красной планеты дополнительно разложен по осям х и у); Ь) графическое отображение потери точности численного решения системы ОДУ
Пример 2 (тоже пока тестовый). Хоровод планет.
На вершинах пятиконечной звезды помещаются пять материальных точек (планеты) с одинаковой массой. Начальные скорости точек равны по модулю и
В. Ф. Очков, Е. П. Богомолова, Д. А. Иванов, К. Писачич
имеют направления по часовой стрелке по касательной к окружности, куда вписана звезда (рис. 5а). Движение планет с такой стартовой позиции и с такими стартовыми скоростями (в работе [1] такое движение для трех точек называется хороводом — см. https://www.ptcusercommunity.com/videos/5782) отображено на рис. 5Ь.
а)
Ь)
с)
Рисунок 5. Кадры анимации хоровода пяти планет ^йрч://№'№'№.р^шегсоттит^. сот/твччазв/З01560): а) старт; Ь) полное совпадение численного решения с аналитическим; с) сбой в численном решении
«Космический хоровод», как и любой хоровод не может длиться бесконечно долго — в какой-то момент он рассыпается, а его участники разбегаются в разные стороны. Что-то подобное можно наблюдать и в нашем хороводе планет при увеличении значения переменной численного интегрирования t (рис. 5сб).
В задачу о движении материальных точек можно внести и третье измерение, переведя ее из плоскости в объем («в космическое пространство»7). Для этого необходимо и достаточно записать еще одну тройку уравнений в систему, показанную на рис. 1, и добавить начальных данных. На рис. 6 показан «хоровод» трех планет в объеме. Третья ось направлена вверх (по ней задано одинаковое начальное равномерное движение планет вверх по оси z), а кривые являются графиками трех решений. Этот график можно также рассматривать как изменение местоположения точек на плоскости (см. рис. 4а) с течением времени (такие графики называются интегральными кривыми).
Рисунок 6. Трехмерная задача о движении трех материальных точек (https://www.ptcusercommunity.com/thread/60224)
6 Такие «произведения изобразительного искусства» часто создаются детьми: сначала все идет хорошо, а потом ребенок устает рисовать правильно и малюет на почти готовом рисунке «каляку-маляку». Можно сказать, что МаШса(! «устал рисовать нам красивый космический цветок» (см. рис. 55) и закончил все «калякой-малякой» (см. рис. 5с).
7 Мы говорим «полеты в межпланетном пространстве». А правильнее нужно говорить «полеты на межпланетной плоскости». Большинство расчетов полетов естественных и искусственных космических тел ведется в двух измерениях с выбором правильной плоскости.
В. Ф. Очков, Е. П. Богомолова, Д. А. Иванов, К. Писачич
На рис. 4a фактически показана одна из трех проекций движения трех точек (см. рис. 6), когда ось z направлена на наблюдателя. Траектории, отображенные на рис. 6, можно интерпретировать как заплетание трех прядей волос в косичку, что хорошо видно при анимации этих траекторий. Если же начальные скорости для трех планет задавать разные, то можно получить очень запутанные траектории движения, из которых уже будет трудно выделить составляющие, связанные с самой физикой задачи и с потерей точности численного ее решения. Аналитическое (абсолютно точное) решение тут, повторяем, невозможно.
Но вернемся к «плоским» случаям, имея в виду, во-первых, то, что планеты Солнечной системы, например, движутся по орбитам, лежащим примерно в одной плоскости, и, во-вторых, то, что численным методам, встроенным в Mathcad, дай бог справиться с плоскими задачами...
Два наших предыдущих примера («восьмерка» и «хоровод») — это, повторяем, тесты для сравнения численных и аналитических решений, когда ошибка численных методов еще не достигла критического уровня. Адекватность решения при не слишком больших8 значениях t позволяет нам применить данный способ решения (см. рис. 1 и 2) к еще не исследованным случаям, для которых аналитические решения неизвестны и невозможны. Для этого, повторяем, необходимо (и достаточно ли — это отдельный вопрос) добавить новые уравнения с новыми слагаемыми в правой части, изменить начальные условия и выбрать подходящий метод численного решения.
Пример 3. Перехват спутника.
Красная планета в паре с синим спутником движется в сторону черной планеты. У всех разные массы: 2; 0.01 и 20 единиц массы (позже мы перейдем к килограммам и тоннам). Что произойдет с этими небесными телами в момент сближения? Для исследования этого случая мы проведем численные эксперименты и решим данную задачу, используя разные методы решения. На рис. 7 зафиксирован кадр анимации, который показывает, что в момент сближения планет спутник меняет своего «патрона» и начинает вращаться вокруг другой (черной) планеты. Красная же планета продолжает свой небесный путь в одиночестве.
Перехват спутника при создании анимации, показанной на рис. 79, был смоделирован с использованием метода (установки — см. рис. 2) Radau для функции Odesolve. Если же метод решения сменить с Radau на Adams/BDF, то мы увидим качественно другую картину: перехват спутника не состоится, а будет наблюдаться только некоторая коррекция его орбиты вокруг красной планеты вследствие приближения черной планеты к системе «красная планета — синий спутник» (рис. 8).
8 Знать бы еще, что такое в количественном выражении «не слишком большие значения».
9 На рисунке слева показано увеличение области перехвата спутника
Рисунок 7. Кадр анимации перехват спутника (https://www.ptcusercommunity. com/videos/5779)
Рисунок 8. Результаты анимации при использовании двух разных численных методов (Radau и Adams/BDF) решения системы ОДУ (https://www.ptcusercommunity. com/videos/5779)
Из рис. 8 видно, что два разных метода решения системы ОДУ (Radau и Adams/BDF) сначала (до встречи с черной планетой) дают одинаковые результаты, а потом (после встречи с черной планетой) результаты расходятся: при методе
В. Ф. Очков, Е. П. Богомолова, Д. А. Иванов, К. Писачич
Radau, повторяем, наблюдается перехват спутника черной планетой у красной планеты, а при методе Adams/BDF только коррекция его орбиты вблизи старого «патрона» — у красной планеты. Методы же Fixed и Adaptive (интегрирование с постоянным и переменным шагом) тут, как можно увидеть из соответствующей анимации, зайдя на вышеотмеченный сайт, совсем не годятся — срывы решения видны в самом начале процесса численного решения задачи (интегрирования).
Итак, мы имеем реальный физический процесс, его физическую модель, математическую модель физической модели (уравнения и начальные условия) и несколько вычислительных (численных) моделей того же процесса: расчетные модули и формулы для методов Radau, Adams/BDF, Fixed и Adaptive (см. рис. 3). Численные модели дают некоторые приближения математической модели, которая осуществляет приближение физической модели, которая, в свою очередь, приближает реальный физический процесс. На каждом этапе моделирования мы по разным объективным причинам допускаем ошибки, а поэтому необходимо отделить несущественные ошибки от принципиальных. Ошибки указанных численных моделей являются принципиальными: четыре расчетные схемы дали четыре разных результата, которые сильно зависят от метода решения. И не количественно (см. рис. 4 и 5), а качественно (см. рис. 8). Это означает, что задача так и не решена в целом. Хотя первая фаза движения одинакова для двух методов (см. рис. 8), а поэтому, скорее всего, начальные результаты можно использовать при численном моделировании решаемой задачи. Но см. [6] и конец статьи.
Анимации, отображенные в одном кадре на рис. 5-8, еще раз акцентируют наше внимание на «блеске и нищете» численных методов решения ОДУ: они позволяют худо-бедно решать задачи, которые «не по зубам» аналитическим методам, и создавать довольно занимательные и поучительные анимации. Но они же могут нам «врать прямо в глаза». Но к этому, повторяем, мы еще вернемся в конце статьи.
Тут затрагивается важный вопрос — можно ли использовать анимацию для оценки качества тех или иных методов решения систем ОДУ. Дерганье траекторий, «каляки-маляки» свидетельствуют о том, что использованный метод решения не годится, по крайней мере, для выбранного диапазона интегрирования. С другой стороны, правильность траекторий планет и спутников, их похожесть на реальные процессы может служить неким признаком добротности использованных методов решения. Но не всегда — см. рис. 8, где обе траектории синего спутника «красивы» и похожи на реальные траектории.
Пример 4. Обмен спутниками.
На рис. 9 показан еще один нетривиальный и более сложный случай — обмен спутниками. Случаи, показанные ранее на рис. 7 и 8, имеют житейские аналогии: существует некая супружеская пара (система «планета-спутник»), к которой при-
ближается некий посторонний объект, вернее, субъект. Он может либо «отбить» спутник (см. рис. 7), либо внести в жизнь этой супружеской пары некие «пертурбации» (см. рис. 8). Когда на сайте https://www.ptcusercommunity.com/videos/5783 публиковалась анимация, показанная на рис. 9, то она была нами прокомментирована так: «В жизни тоже такое встречается: две супружеские пары знакомятся, сближаются и в конце концов понимают, что когда-то они сделали не тот выбор; происходит законный (два развода и две новые свадьбы) или неформальный обмен партнерами». На этот комментарий сразу отозвался модератор сайта и попросил его убрать, так как «на этот сайт заходят и несовершеннолетние», школьники, например, решающие свои домашние задания в среде Mathcad.
Рисунок 9. Кадр анимации обмена спутниками (https://www.ptcusercommunity.com/thread/59968)
У рис. 7-9 несколько фривольная интерпретация и потому, что нельзя серьезно относиться к этим «житейским» траекториям из-за, повторяем, сомнений в адекватности численных методов, использованных для решения этих задач. Но случаи интересные и базирующиеся на «адекватной» физике. И с этим нужно согласиться. Вернее так: можно согласиться, а можно не согласиться.
Пример 5 (гибридный: проверка старого и показ нового). Хоровод (па-де-катр) четырех планет с четырьмя спутниками.
На рис. 10а показан стартовый кадр анимации системы из восьми материальных точек: четырех планет (масса каждой 0.5) и четырех спутников (масса каждого 0.0001). Стартовые условия такие: планеты помещаются в углах квадрата со сторо-
В. Ф. Очков, Е. П. Богомолова, Д. А. Иванов, К. Писачич
ной, равной единице, и им придаются скорости, равные 0.5 и направленные вверх (синяя планета), вправо (черная), вниз (красная) и влево (коричневая). Такие уже пять, а не четыре планеты без спутников кружатся в «хороводе», показанном на рис. 5. Но в нашем случае в серединах сторон квадрата дополнительно помещаются спутники с нулевой начальной скоростью.
с) ё)
Рис. 10. Кадры анимации движения 4 планет и 4 спутников: а) старт; Ь) полет спутников вокруг планет; с) вылет спутников из «хоровода»; d) траектория одного спутника (остальные три траектории невидимы)
Анимацию движения этих восьми небесных тел описать словами почти невозможно. Нужно смотреть саму анимацию — https://www.ptcusercommunity.com/ thread/84910. Каждый из спутников будет стараться поочередно «пристроиться» к каждой из четырех планет, выделывая всевозможные «па» (рис. 106). В конце «танца» эти спутники вылетят из «хоровода» (рис. 10с). На вышеотмеченном сайте можно видеть анимацию движения одного из спутников, когда три остальные спутника не показаны (рис. 10d). Там четко видно, как спутник меняет партнеров в этом «танце». Такие же «телодвижения» в это время производят и другие спутники с единственной разницей: траектории их полетов повернуты относительно центральной точки (центра масс) на 90° (рис. 10с). Остается открытым вопрос, вернутся ли назад в «хоровод» отлетевшие спутники. Численным экспериментом на этот вопрос ответить затруднительно — тут всю картину испортит погрешность вычислений (см., например, рис. 56). Попытка ответа (вернее, один из путей решения этой отдельной задачи) будет дана в дивертисменте к статье. Симметрии же полета планет и спутников, зафиксированные на рис. 10с, дают основание делать вывод о достаточной точности численного решения задачи.
Пример 6. Спутник Земли или Часы Кеплера.
Новая версия пакета Mathcad — Mathcad Prime пока не имеет анимации, но зато может решать ОДУ с единицами измерений физических величин. До этого (Mathcad 15) эти единицы игнорировались и использовались безразмерные единицы длины, скорости, массы и т. д., что чревато ошибками: где-то в системе уравнений по оплошности тройку заменили на двойку, система, тем не менее, решается, но выдает неверный ответ. В среде Mathcad Prime такая ошибка сразу будет замечена, при этом ведется пересчет единиц измерения (что отмечено ниже).
На рис. 11 показана орбита вращения спутника Земли с реальными значениями расстояний, скоростей и, главное, гравитационной постоянной G. Подобраны такие исходные данные, чтобы период вращения спутника был равен 12 часам. Кроме того, на графике добавлены сектора эллиптической орбиты, которые иллюстрируют второй закон Кеплера: площади этих секторов в равные (часовые) промежутки времени равны (с точностью сотых процента — см. рис. 116) и при численном решении задачи, что является дополнительным подтверждением верности нашей мо-дели10. Саму же анимацию можно рассматривать, как еще одни необычные математические часы [7], у которых стрелки движутся по эллиптическому циферблату с
10 Задача о двух телах имеет аналитическое решение, в частности, эллипс как траектория полета. Наш эллипс можно привести к каноническому виду и рассчитать значения его параметров — коэффициентов уравнения кривой второго порядка, осей, фокусов и прочее. Эту задачу мы попытаемся решить в дивертисменте.
В. Ф. Очков, Е. П. Богомолова, Д. А. Иванов, К. Писачич
переменной скоростью11. Такие часы будут хорошо смотреться в музеях Кеплера, планетариях, обсерваториях...
11 Попробуем составить такую формулу изобретения. Часы Кеплера представляют собой традиционные часы со стрелками и циферблатом. Отличие состоит в том, что циферблат часов Кеплера — это не окружность, а эллипс с числами 1, 2, 3. 12 (I, II, III. XII). Стрелки таких часов (часовая, минутная и секундная) движутся по кругу неравномерно. Длины стрелок также меняются по мере их движения. Концы стрелок описывают эллиптические траектории. Дополнительно на циферблате отображается сектор, охватывающий два соседних часа, между которыми в данный момент времени находится секундная стрелка. Этот сектор меняет свое расположение и форму каждые пять секунд, но площадь его остается постоянной согласно второму закону Кеплера. Вокруг места крепления стрелок прорисована окружность, изображающая Землю, а на конце секундной стрелки помещен кружок, изображающий спутник Земли. Такие часы могут висеть в школьных кабинетах математики, в планетариях, обсерваториях, музеях, привлекая внимание людей, заставляя их задуматься о законах Небесной механики, в частности, о втором законе Кеплера.
■ Play Animation
Krepier K-'S (Joe. к XI. XII , -----.II Y»l*>r OcMov * Matbcad 13 1 f6hr
пчп- 26 X \ III S«3: \ IV -■ I r(t)-v(t) dt= 5103B3krtT 2 -'shr
IX / \ 510383- 510069 ..... - 0.062% 510383
I \ ' / v 1 Г1№Г
\ /VI s<:= - r(t) v(t) dt ^ 510069 km" 2 MOhr
у r-24553 km
yS v-1S05dkpb ^^ VII m-IOlorine
VIII m- 5.972 x 1021 tonne
► 1 g |
b)
Рисунок 11. Часы Кеплера: a) Mathcad-расчет; b) кадр анимации (https://www.ptcusercommunity.com/thread/60078)
Из рис. 11b видно, что спутник в районе «восьми часов» максимально приближен к Земле (перигей) и в этой точке будет иметь максимальную скорость согласно второму закону Кеплера. В районе же двух часов (апогей) скорость спутника будет минимальной.
Вблизи Земли можно отходить от закона всемирного тяготения, переходить к потенциальному полю тяжести вблизи нашей планеты, характеризуемого ускорением свободного падения, и учитывать ее атмосферу. Такую задачу также несложно решить с привлечением Mathcad (рис. 12).
Пример 7. Прыжок парашютиста.
Постановка задачи о парашютисте возникла под влиянием знаменитого прыжка с парашютом с высоты почти 40 км (из стратосферы) австрийца Феликса Баум-гартнера (http://lenta.ru/news/2014/02/02/video/). Подобный прыжок с еще большей высоты в октябре 2014 г. повторил Алан Юстас — исполнительный директор фирмы Google (http://www.profile.ru/mir/item/87941-ispolnitelnyj-direktor-google-pobil-rekord-vysoty-dlya-pryzhkov-iz-stratosfery).
Рассмотрим такую математическую модель: парашютист в начале полета — это шар с радиусом r и массой 110 кг. После раскрытия парашюта он превращается в шар с радиусом r2. Масса парашютиста не меняется. Прыгает парашютист с высоты h, а раскрывает парашют на высоте И2. Эти исходные данные вводятся в Mathcad-расчет в виде таблицы (см. рис. 12), у которой первая строка — это имена переменных, вторая — единицы измерения, а третья — числовые значения.
В. Ф. Очков, Е. П. Богомолова, Д. А. Иванов, К. Писачич
Рисунок 12. Моделирование прыжка парашютиста
В задачу введены и дополнительные исходные данные: плотность воздуха рвозд
как функция высоты над уровнем моря (взято из Википедии), коэффициент трения парашютиста о воздух к (взято «из головы», но уточняется в процессе решения и сравнивания полученного ответа с реальной максимальной скоростью парашютиста) и предположительное время полета парашютиста ^ (это значение уточняется в процессе решения задачи). После ввода исходных данных записываются три вспомогательные функции: изменение в зависимости от высоты полета радиуса «парашютиста» г, площади его поперечного сечения Сечение и его объема Объем. Эти функции имеют вид ступенек: до раскрытия парашюта (к > к2) они возвращают одни значения, а после раскрытия — другие.
На рис. 12 можно видеть дифференциальное уравнение, «уравнивающее» силы, действующие на парашютиста. Это уравнение не что иное, как математическая запись второго закона Ньютона, гласящего, что сумма сил, действующих на тело, равна произведению его массы на ускорение.
Какие силы действуют на парашютиста?
Первая сила — сила тяжести: произведение массы на ускорение свободного падения g. Эта сила направлена вниз, поэтому-то она в уравнении со знаком минус. Вторая сила — это сила сопротивления воздуха, которой обычно пренебрегают при рассмотрении классической задачи полета (падения) камня, но которой никак нельзя пренебречь в случае с парашютистом, иначе он разобьется в лепешку. Мы примем, что эта сила пропорциональна плотности воздуха, площади поперечного сечения падающего тела и квадрату его скорости. Коэффициент пропорциональности — это безразмерная переменная к. Третья сила — это сила Архимеда — вес вытесненного парашютистом воздуха. Если парашют еще не раскрыт, то этой силой можно пренебречь, но если наш парашютист раскрыл парашют, «раздувшись» (согласно нашей модели) до шара с диаметром 5 м и массой в 110 кг, то эту силу учитывать уже нужно. Две эти силы направлены вверх, поэтому в их формулах мы видим знак плюс. Все эти три силы, повторяем, уравновешиваются силой инерции — произведением массы на ускорение.
Задачу можно усложнить — учитывать, например, изменение значения ускорения свободного падения по высоте, не шарообразную, а более сложную форму летящего парашютиста, изменение значения коэффициента к в зависимости от режима обтекания тела (ламинарный или турбулентный — все это является предметом изучения науки аэродинамики), горизонтальную составляющую полета парашютиста, связанную со скоростью самолета, из которого он выпрыгнул, и/или со скоростью ветра. Но и без этого наше решение получилось вполне правдоподобным: парашютист в свободном полете набирает максимальную скорость до 875 км/ч, затем его скорость падает вследствие роста плотности воздуха. После
В. Ф. Очков, Е. П. Богомолова, Д. А. Иванов, К. Писачич
раскрытия парашюта скорость парашютиста резко уменьшается до 11 км/ч. С такой скоростью он и приземляется после 13 минут полета. Даже в школьном курсе физики известны задачи, рассматривающие движение тел с установившейся скоростью, которое возникает при действии на тело силы, направленной противоположно скорости тела, причем модуль силы увеличивается с увеличением скорости тела. Можно снять еще одно допущение. У нас при h > h2 парашют раскрывается моментально. В реальности это происходит за одну-две секунды, что можно учесть, сделав функцию r(h) более сложной: не ступенькой, а неким пандусом.
На сайте, поддерживающем эту статью, есть расчет и анимация приземления спутника (https://www.ptcusercommunity.com/thread/84707) за счет кратковременного запуска двигателя, но без учета атмосферы Земли. Эту задачу несложно скрестить с задачей, показанной на рис. 12, и получить модель, близкую к реальности.
Пример 8. Прыжок с вышки в воду.
Модель и ее реализация в среде Mathcad прыжка парашютиста на Землю, вернее на землю (не на планету, а на грунт), показана на рис. 12. Но парашютист или космический аппарат может не только приземляться, но и приводняться. На заре космической эры советские аппараты приземлялись, а американские приводнялись. На рис. 13 показан расчет параметров прыжка человека в воду с вышки. Здесь также присутствует функция-ступенька, но связанная уже не с объемом «прыгуна», а с плотностью среды: воздух или вода. В расчете есть еще одна функция-ступенька — встроенная в Mathcad функция с именем sign, которая возвращает знак аргумента: минус единицу при отрицательном, нуль при нулевом и плюс единицу при положительном аргументе. Эта встроенная функция позволяет нам учитывать то положение, что сила сопротивления среды всегда действует в направлении, противоположном скорости. Сам же знак значения скорости у нас пропадает из-за возведения этой величины в квадрат.
На рис. 13a можно видеть графики параметров прыжка с вышки в воду. На одном из графиков отображен так называемый аттрактор: прыгун, вынырнув из воды, будет совершать затухающие колебательные движения на ее поверхности. Такой аттрактор можно увидеть на срезах деревянных досок с годовыми линиями. Анимация этой модели (см. рис. 13b) размещена по адресу https://www.ptcusercommunity.com/docs/DQC-3013.
Но в воду можно не только нырять — из воды можно также и «выныривать». Если вернуться к главной теме статьи, то тут будет уместен еще один «ракетный» пример.
а)
Ь)
Рисунок 13. Моделирование прыжка с вышки в воду: а) Mathcad-решение; Ь) кадр анимации
В. Ф. Очков, Е. П. Богомолова, Д. А. Иванов, К. Писачич
Пример 9. Старт ракеты из подводной лодки.
Космические объекты могут не только приводнятся, но и стартовать с поверхности воды, и даже из-под воды.
На рис. 14 можно видеть кадр анимации вертикального взлета и последующего падения в воду ракеты, запущенной из-под воды. Эта математическая модель отличается от предыдущих тем, что к силам, действующим на материальную точку, прибавилась сила тяги ракетного двигателя, а также тем, что масса материальной точки уменьшается вследствие расхода топлива и окислителя.
Рисунок 14. Кадр анимации полета ракеты с подводной лодки (https://www.ptcusercommunitv. сот/йос8/РОС-3570]
На рис. 14 показаны также графики изменения высоты подъема ракеты (h — красный пунктир) и ее скорости (v — синий пунктир) в зависимости от времени полета. На сайте https://www.ptcusercommunity.com/thread/60347 можно видеть анимацию еще одной подобной задачи, когда спутник изменяет свою орбиту за счет кратковременного включения бортового двигателя.
Рассмотренные инструментальные средства Mathcad позволяют легко, быстро и с достаточной точностью, отображать ну если не количественную, а только качественную сторону задачи, моделировать и анимировать различные случаи движения материальных точек и тел под действием сил гравитации, инерции, тяги двигателя и сопротивления среды. Предлагаем читателям рассмотреть такую модель: космический корабль на орбите Земли включает ракетный двигатель с известной тягой, стартует к Луне (https://www.ptcusercommunity.com/thread/84708) и после свободного полета при выключенном двигателе переходит на орбиту этого спутни-
ка Земли (https://www.ptcusercommunity.com/thread/84709). Если же тягу двигателя направить против скорости движения спутника Земли, то он «войдет в плотные слои атмосферы», сгорит или успешно приземлится, или приводнится с помощью парашюта. Элементы этих случаев движения тел мы рассмотрели выше, и их анимация есть на описываемом сайте.
При этом нужно только помнить одну важную вещь. Перефразируя известное выражение, можно сказать, что есть «ложь, наглая ложь и... численное решение дифференциальных уравнений». В этом можно еще раз убедиться, взглянув на рис. 5-8. Ни в коей мере не утверждаем, что результаты численных методов не заслуживают доверия. К ним просто нужно относиться очень осторожно. Как, например, к статистике, которая должна стоять в вышеприведенном крылатом выражении [8]. Численные методы решения дифференциальных уравнений уже давно успешно используются при осуществлении космических полетов. При этом данные телеметрии используют в том числе и для коррекции движения, которое было численно рассчитано заранее с определенной, часто недостаточной точностью. При этом и начальная скорость аппарата, заданная маршевыми двигателями, может отличаться от задуманной. Так и мы, принимая к исполнению, чью-то просьбу или приказание, понимаем, что в этой просьбе или приказании может быть вольная или невольная ложь, которую мы нивелируем в процессе выполнения просьбы, получая дополнительную информацию и принимая самостоятельные решения.
В настоящее время многие физические кабинеты школ и университетов оборудованы компьютерами с мультимедийными проекторами. В таком кабинете можно проводить физические опыты (изучать, например, колебание маятника), показывая на большом экране эксперимент так, чтобы всем было хорошо все видно. Но на этом же экране можно показывать решение дифференциального уравнения колебания маятника, сравнивать реальное физическое явление с его математической моделью, объясняя расхождения ограничениями и упрощениями модели.
Уравнения — и алгебраические, и дифференциальные — страшны для школьников и студентов не сами по себе, а методами их решения. Сейчас на компьютере решать такие уравнения можно довольно просто. Главное — составить уравнение или систему уравнений, понимая физику задачи, и правильно интерпретируя результаты расчетов. Поэтому такие задачи могут стать не «пытками», а удовольствием и для учеников, и для учителей.
На сайте http://communities.ptc.com/groups/dynamic-models-in-mathcad авторы разместили большое количество различных динамических задач с их решениями в среде Mathcad [9] и анимацией средствами Mathcad [3]. Вот некоторые из них:
- колебание одиночного маятника;
- колебание связанных маятников;
В. Ф. Очков, Е. П. Богомолова, Д. А. Иванов, К. Писачич
- вращение планет со спутниками (тема данной статьи; этим решениям посвящен отдельный подфорум https://www.ptcusercommunity.com/ groups/three-body-problem-with-animation);
- старт ракеты из подводной лодки;
- ныряние человека в воду;
- скатывание саней с горки;
- движение автомобиля;
- движение по подземному прямому тоннелю гравитационного поезда [10];
- и др.
Дивертисмент
«Как беззаконная комета в кругу расчисленных светил... »
А.С. Пушкин «Портрет»
На рис. 11 отображено движение спутника по эллиптической орбите. Но есть небесные тела («беззаконные кометы», например), которые движутся по иным траекториям, и это отображено в эпиграфе к дивертисменту.
А траектории, повторяем, бывают такие: эллипс (частный случай — окружность), гипербола и парабола (редкий переходный случай от гиперболы к эллипсу).
Из наших расчетов, основанных на численном решении ОДУ с начальными условиями12, часто можно видеть только какую-то дугу траектории, по которой трудно определить, что это — эллипс, гипербола или парабола. Но это можно понять, если вспомнить, что рассматриваемые кривые — это плоские кривые второго порядка. В уравнениях этих кривых (см. верхнюю часть рис. 17) шесть коэффициентов. Если попытаться их вычислить, то по ним (графически или через вычисление инвариантов) можно определить тип кривой.
Тут на ум приходит еще одна цитата из А.С. Пушкина, но уже из «Каменного гостя»:
Дон Гуан: Ее совсем не видно
Под этим вдовьим черным покрывалом, Чуть узенькую пятку я заметил. Лепорелло:
Довольно с вас. У вас воображенье В минуту дорисует остальное; Оно у нас проворней живописца,
12 Можно попытаться решить и краевую задачу, но это выходит за рамки статьи.
Вам все равно, с чего бы ни начать,
С бровей ли, с ног ли.
Но у нас есть не только «воображенье», но и мощные вычислительные средства решения систем алгебраических уравнений (см. рис. 17). «Узенькая пятка» — это наши точки численного решения системы ОДУ (рис. 15 и 16), а «остальное» — это вся гипербола, вернее, коэффициенты ее уравнения. А то, что это гипербола, доказывает значение инварианта D, которое меньше нуля (у эллипса оно больше нуля, а у параболы равно нулю), и соответствующие графики на рис. 18.
На рис. 15 показан блок исходных данных и блок численного решения системы ОДУ, отображающий движение материальной точки (спутника) вблизи другой материальной точки (Земли) с массой, намного превышающей массу первой точки так, чтобы перемещением второй точки можно было пренебречь13.
в Еаг1Ь Я х0 у0
"уО
| (шпе) (Ш) (М) (М) т
6,67384-10 11 5.972-1021 6370 2 К 2 Й -0.5 -7 12
Х(05)=Х0 х'(0 5) = ^ х"(0 = -6-—я
+у( О2)
у(0 5)=у„ У'(Os) = vy0 у"Ц) = -в---ЕаПЬ-уЦ) ^
Гх1 /ГХ(Л1 \ +у( о')
Рисунок 15. Блоки исходных данных и численного решения системы ОДУ
На рис. 16 показано графическое отображение решения задачи о движении небесного тела малой массы (астероида, например) вблизи Земли за 12 часов. Из рис. 16 довольно трудно понять, чью дугу (часть) мы получили — эллипса, гиперболы или параболы. Можно увеличивать предел интегрирования14 (увеличивать значение переменной /епа) и видеть, будет ли кривая загибаться (сворачиваться) в эллипс или уходить в бесконечность, помня при этом о точности решения (см. рис. выше). Но можно поступить иначе — восстановить всю кривую второго порядка по коэффициентам ее уравнения (см. рис. 17).
13 В задаче, показанной на рис. 11, мы этого не делали, 1 Но, повт( и 5 с выше).
14 Но, повторяем, нужно помнить о том, что это часто приводит к накоплению ошибки (см. рис. ^
В. Ф. Очков, Е. П. Богомолова, Д. А. Иванов, К. Писачич
Рисунок 16. Графическое отображение численного решения системы ОДУ
а1Гх +2 а)2->:.у+э22.у* +2 а1-х+2 а2.у+а0=0
[ - г '' 1. Э2 Эц]:-
1 1
2 2 т т
1 1
1
а„.х(0 Лг) +2 а12-х(0 йг)-у(0 Ьг) + э22-у(0 Лу) +2 агх(0 Лг) + 2 а2.у(0 Лг) + ас = О 2 2 э„-*(1 йг) +2 э,2-х(1 /1г) - у (1 ^г) + а22-у(1 Л/") +2 йу) + 2 а2-у(1 Ьу) + э0 = С
2 2 а1Гх(2 Лг) +2 э,2-х(2 Ьг)-у(2 /7г) + э22-у{2 Лу) +2 а,-х(2 Лг) + 2 а2.у(2 *1г) + Э£, = С 2 2 а„.х(3 Лу) +2 з,2 -зг(3 /1г).у(3 /7г) + а22.у(3 Ау) +2 агх(3 Лу) + 2 а2.у(3 Ьг)+ав=0 2 2 а11-х(4 Лу) +2 э;2-х(4 йг)-у(4 ^г)-|-э22-у(4 Ьг) +2 а1-х(4 Лу) + 2 а2-у(4 Ьу) + э0 = С 2 2 а,,-*(5 Лг) +2 а,2-х(5 йг)-у(5 ?)г) + э22-у(5 Лу) +2 а,-х(5 Лг) + 2 а2-у(5 ^г) + эо = 0 2 2 а„.х(6 Лу) +2 з,2.1-(е /1г).у(В Лг) + а22.у(6 Лу) +2 агх(в Лу) + 2 а2.у(6 Лг) + ас = О 2 2 а)(-х(7 Лу) +2 йг)-у(7 Лг) + э22-у(7 Ьг) +2 а,->;"(7 Лу) + 2 а2-у(7 Ьг) + э0 = С
2 2 а„-*(8 Лг) +2 а,2 -х(8 йг)-У(8 И + э22-у(8 Лу) +2 а,-х(8 Лг) + 2 а2-у(8 Ьг) + эо = 0 2 2 а„.х(9 Лг) +2 з,2 -зг(9 йг).у(9 Ьг) + а22.у(9 Лу) +2 агх(9 Лг) + 2 аг.у(9 Ьг)+ав=0 2 2 а,гх(10 Лу) +2 а)2.х(10 Лу) -у(10 Лу)-ьэ22-у(10 Лг) +2 аг*(10 Лу) + 2 агу(10 Ьг) + ао = 0
а,г*(11 Лг)2+2 э12-х{11 Лг).у(11 Лг) + а22.у(11 +2 а,-х( 11 Лг) + 2 а2-у(11 Лу)+а0 = О
а„.х(12 Ьг)2 +2 а12.х(12 Лу)-у( 12 Лу) + э22-у(12 Лу)2 +2 агх(12 Лу) + 2 а,.у(12 Ьг) + ао=0
(101.180524192814-10"15) Д-
{138.877437855728-10 15}
Л|'
(78.4361030434 229-10"15)
т
(2.9463522233284.10"°) — т
-3.1974809672896-10"° — т
-67.8361105562333
a)
Ь)
М :=
х(7 Ьг) 2 х{7 Лг)-у{7 Лг) у(7 (.',') 2 х(7 Лг) 2 у(7 Лг)
х(8 Лг)2 2 х(8 /?г).у(8 Лг) у (8 Лг)2 2 х(8 /)г) 2 у (8 Лг)
х(9 Лг) 2 х(9 пг) ■ у(9 Лг) у (9 Лг) 2 х(9 Лг) 2 у (9 Лг)
х(ЮЛг)2 2 х(10 Лг)-у{10 у(10Лг)2 2 х{10 йг) 2 у(10 Лу) 2 2
х(11 Лг) 2 х(11 йг)-у(11 Лг) у(11 Лг) 2 х(11 Лг) 2 у(11 Лг)
(1.49154295494464- ю-'5) 1 2 т
(2.04724825269486 1 о'5) 1 2 т
(1.15625764036381 1 2 т
(43.4333638807231 ■ ю-9) 1 т
-47 1353781403196-Ю " 1 т
с)
Рис. 17. Численное решение системы алгебраических уравнений: а) решение однородной системы из 13 уравнений с множеством пропорциональных решений; Ь) решение системы пяти уравнений с пятью неизвестными и единственным решением; с) решение системы линейных алгебраических уравнений с пятью неизвестными и единственным решением
На рис. 17 показано три решения. Первое (см. рис. 17а) можно назвать «лобовой атакой»: даны начальные приближения (предположения) к решению, записаны 13 уравнений, фиксирующих положение спутника в начальный момент времени
В. Ф. Очков, Е. П. Богомолова, Д. А. Иванов, К. Писачич
(t = 0), после часа полета, после 2 часов полета и так далее до 12 часов полета и вызвана функция Find, выдавшая значения своих аргументов, обращающих уравнения в тождества. Если дать иные начальные приближения к решению, то будет выдан иной ответ, который тоже будет правильный. При этом значения шести неизвестных коэффициентов при первом наборе начальных приближений будут пропорциональны значениям этих же неизвестных при втором наборе начальных приближений. И это понятно, поскольку все коэффициенты в уравнении кривой, приведенной на рис. 17а, можно умножить на одно и то же число и получить новое уравнение относительно x и y, равносильное исходному. Это свойство позволяет заранее зафиксировать значение одного из коэффициентов и найти оставшиеся. Таким образом, задавая различные начальные приближения для всех шести неизвестных коэффициентов, мы получаем различные, пропорциональные друг другу наборы коэффициентов, определяющих уравнение одной и той же кривой второго порядка.
На рис. 17b и 17с в задаче оставлено только пять уравнений, так как известно, что кривую второго порядка можно восстановить по пяти точкам. При этом коэффициенту a0 присвоено значение минус единицы. В этом случае решение будет единственным, и начальные приближения для оставшихся пяти коэффициентов уже не повлияют на ответ (см. рис. 17b), а саму задачу можно будет свести к решению системы линейных алгебраических уравнений, не требующему начальных приближений (см. рис. 17с). Линейную систему в среде Mathcad тоже можно решить с помощью функции Find (см. рис. 17b), но для этой задачи в среде Mathcad есть специализированные средства, например, умножение инвертированной матрицы коэффициентов при неизвестных на вектор свободных членов (см. рис. 17с).
На рис. 18, во-первых, показано, что главный инвариант D уравнения плоской кривой второго порядка меньше нуля и, во-вторых, построены две половинки гиперболы — верхняя (top, синяя) и нижняя (ир, черная). Аналитические выражения для них получены в результате решения уравнения кривой относительно переменной y и генерации двух функций с именем ytop и yup . На верхнюю ветвь гиперболы
(на нижнюю, черную половинку) нанизаны наши красные исходные точки.
Данный расчет можно продолжить, найдя, например, значения коэффициентов а и b канонического уравнения гиперболы:
x2/a2 - y2/b2 = 1.
Первым шагом такого преобразования будет нахождение угла, на который нужно повернуть оси декартового графика (см. рис. 19).
Рисунок 18. Построение гиперболы по ее уравнению
0:=гоо1 НапО)2 - Эг2 Хап(ф)-1 ,ф, 0 с/ед, 90 с1ед= 47.341 йед
Рисунок 19. Угол поворота осей
Предлагаем читателю самому закончить такое преобразование, имея под рукой все шесть коэффициентов а уравнения плоской кривой второго порядка. Угол поворота двух ветвей гиперболы мы уже нашли (см. рис. 19 — поиск одного из двух корней квадратного уравнения в интервале от 0 до 90 угловых градусов). Осталось только определить значение сдвига ветвей гиперболы, чтобы она приняла канонический вид — симметричный относительно осей и начала координат.
Пять левых красных точек, показанные на рис. 18, можно рассматривать как данные телеметрии движения некоего астероида к Земле. Ученые, имея эти данные, могут рассчитать, как дальше будет лететь этот астероид: пролетит ли он мимо Земли (а именно это показано на рис. 16 и 18) или столкнется с Землей. Для предотвращения такой глобальной катастрофы на данный космический объект мо-
В. Ф. Очков, Е. П. Богомолова, Д. А. Иванов, К. Писачич
жет быть, к примеру, послана команда для взрыва этого астероида — см. кинофильм «Армагеддон».
Когда статья была уже готова к публикации (май 2015 г.), появилось сообщение, что космический грузовик «Прогресс», запущенный к Международной космической станции, не вышел на расчетную орбиту и должен через некоторое время упасть на Землю. Это побудило авторов попытаться смоделировать такое поведение неуправляемого искусственного спутника Земли (рис. 20, 21 и 22).
На рис. 20 в среде Mathcad Prime вводятся исходные данные: гравитационная постоянная G, параметры Земли (ее масса me и радиус R), параметры искусственного спутника (его масса ms, плотность ре и высота над поверхностью Земли в начальный момент h0), коэффициент трения спутника о воздух k и время полета iend. Далее в расчете, показанном на рис. 20, вводится функциональная зависимость плотности воздуха от высоты p(h) (см. также рис. 12), рассчитывается скорость вращения спутника по круговой орбите v0 на высоте \ над уровнем Земли (первая космическая скорость) и параметры искусственного спутника: его объем, диаметр и площадь поперечного сечения в предположении, что это шар с радиусом R и плотностью ре.
Рисунок 20. Моделирование торможения спутника в плотных слоях атмосферы: начало расчета
На рис. 21 показаны операторы решения системы четырех обыкновенных дифференциальных уравнений, моделирующих движение и снижение искусственного спутника на орбите Земли. К силе космического притяжения, определяемого
законом всемирного тяготения, прибавилась сила торможения, связанная с трением о воздух. Эта сила пропорциональна (коэффициент к) площади поперечного сечения спутника ^, плотности воздуха р и квадрату скорости (см. также рис. 12, 13а и 136). При этом в расчет вставлен не квадрат скорости, а произведение скорости на абсолютное значение скорости для того, чтобы правильно определить направление этой силы — она действует в направлении, противоположном направлению вектора скорости.
Solve
x(Os) = 0 т v, (f) = /{f) ms v;(t)=-G- ms-me-x(t) -s * |'S5-p(Vx(f)!+y(f)2-R8)
МО s) = vt> (Vx(f)2+y(;)2)
у(0 =R,+h0 vy V/(t) = -G- -5—*-vy(i)-|vyM |.Ss.p(Vx(f)!+y(f)2-flJ
МО *) =om s (WW)
* [*(f>
У vx Odesolve y(f) "«(<) .Ui.w5
к км
Рисунок 21. Моделирование торможения спутника в плотных слоях атмосферы: продолжение расчета
На рис. 22 показано графическое отображение решения задачи о снижении орбиты спутника до момента касания поверхности Земли: зависимость высоты над уровнем моря h(t) от времени и фрагмент траектории орбиты с тремя витками.
На сайте https://www.ptcusercommunity.com/groups/three-body-problem-with-animation есть модель и анимация случая, обратного, показанного на рис. 33, когда за счет включения двигателя спутник переходит на новую, более высокую орбиту.
Все файлы приведенных в статье расчетов (Mathcad 15 и Mathcad Prime 3) можно скачать с сайта https://www.ptcusercommunity.com/groups/three-body-problem-with-animation. Там же размещены соответствующие анимации. Отдельные их адреса и кадры были приведены выше.
И предпоследнее для самых любопытных и дотошных
Какие алгоритмы заложены в численные методы решения обыкновенных дифференциальных уравнений?
В. Ф. Очков, Е. П. Богомолова, Д. А. Иванов, К. Писачич
Самый простой метод решения таких задач предложил еще великий Эйлер15. Если говорить о задаче с начальными условиями16 (задача Коши), то метод Эйлера реализуется следующим образом.
Рисунок 22. Моделирование торможения спутника в плотных слоях атмосферы: окончание расчета
Если известно начальное (первое) значение искомой функции, то отрезок интегрирования разбивается на п коротких отрезков длиной А и рассчитывается значение искомой функции во второй точке по значению производной в первой точке —
15 Сейчас у нас все чаще говорят не Эйлер, а как во всем мире — Ойлер. Такая метаморфоза случилась и с другим великим математиком — с Ньютоном. Сначала (по времена Ломоносова) его у нас называли Невтоном (как написано (Newton), так и произносилось), затем Ньютоном с ударением на последнем слоге, а сейчас как в Англии — с ударением на первом слоге.
16 А еще есть так называемая краевая задача, когда условия задаются на концах отрезка интегрирования. О такой задаче мы уже упоминали, и ее тоже можно было рассмотреть в статье, но статья и так уже стала очень объемной. Ниже будет кратко описан один из способов решения этой задачи.
а это правая часть дифференциального уравнения: у2 = у + у \ А. Эта операция проводится в цикле (у = + у А)17, в котором формируются два вектора, хранящие дискретные значения аргумента и искомой функции. Если число точек п достаточно большое и искомая функция не слишком «крутая»18, то метод Эйлера дает вполне приемлемые решения. Мы сейчас описали приложение метода Эйлера к одному дифференциальному уравнению первого порядка. В наших же задачах небесной механики рассматривалась система уравнений второго порядка: в математических моделях фигурировало ускорение — вторая производная пути по времени. Метод Эйлера можно применить не только к одному уравнению, но и к системе, а одно уравнение второго порядка можно преобразовать в два уравнения первого порядка. Такое решение показано на рис. 23 применительно к еще одной классической задаче, связанной законом всемирного тяготения и с силой тяжести — к задаче о колебании маятника19.
Математическая модель колебания маятника учитывает тот факт, что его угловое ускорение пропорционально произведению синуса угла отклонения нити маятника от вертикали на ускорение свободного падения, деленное на длину нити маятника. Если синус угла заменить на сам угол (а это можно сделать при значениях угла по модулю меньших, чем приблизительно 7°), то такое уравнение (уравнение математического маятника) можно решить аналитически (рис. 24), где для этого использован сайт программы Mathematica.
Из рис. 23 видно, что при п = 100 (число разбиений интервала интегрирования) решение системы двух ОДУ методом Эйлера получилось довольно грубое, значительно отклоняющееся от истинного. Мы специально на графике отметили линией «истинное» (аналитическое) решение, а численное (приближенное) отметили точками, чтобы подчеркнуть еще раз, что методом Эйлера получена не искомая функция, а набор точек, которые еще как-то нужно преобразовать в функцию, используя, например, методы интерполяции [9].
Наша функция Еи1ег (рис. 23) в случае одного дифференциального уравнения первого порядка возвращает искомые точки в виде матрицы с двумя столбцами — дискретные значения аргумента и функции. Если же мы имеем систему уравнений,
17 В программах на рис. 23 и 25 нижний индекс переменной (элемент вектора) заменен на «верхний индекс» (столбец матрицы с одной строкой). Это было сделано для того, чтобы выводилась итоговая матрица решений (см. форум https://www.ptcusercommunity.com/thread/85205).
18Под словом «крутой» можно понимать очень многое, как в молодежном сленге. В наших задачах спутники и планеты вели себя «круто» и ставили в тупик довольно сложные методы решения дифференциальных уравнений.
19 На сайте PTC Community можно найти авторскую анимацию одиночного маятника в вязкой среде, на эластичной нити, а также анимации двух, трех и четырех связанных маятников на жесткой несгибаемой нити. Маятники при этом описывают замысловатые кривые, подобные тем, какие показаны на рис. 10. И там также можно четко видеть срыв численного решения.
В. Ф. Очков, Е. П. Богомолова, Д. А. Иванов, К. Писачич
то в возвещаемой матрице появятся дополнительные столбцы. На рис. 23 специально показаны два графика (график отклонения по времени нити маятника от вертикали и график угловой скорости маятника) с «извлечением» нужных столбцов матрицы М: М<°>, М<!> иМ<°>.
Рисунок 23. Метод Эйлера, примененный для решения системы двух обыкновенных
дифференциальных уравнений
Если в решении, показанном на рис. 23, значение п увеличить со 100 до 1000 или до еще большего значения, то точки на графике сольются в линию (иллюзия аналитического решения с функцией в ответе) и лягут строго на кривую — решение будет более-менее точным. Но это соответственно увеличит время счета20. Из-за этого и по другими причинам (из-за «крутизны», т. е. очень быстрого изменения
20 Сейчас на мощных компьютерах такое изменение счета почти не чувствуется, но в те времена, когда создавались описываемые методы, и когда вообще не было компьютеров, эта временная разница была существенной и даже определяющей при выборе того или иного метода решения.
искомых функций) были предложены более «хитрые» методы численного решения ОДУ. Самый известный из них разработан в 1900 г. немецкими математиками К. Рунге и М. Куттой. Метод Рунге-Кутты (рис. 25) ведет расчет очередной точки значения искомой функции с некой «пристрелкой» на половинном значении А. Он встроен практически во все компьютерные вычислительные программы, включая и пакет Mathcad — см. выше рис. 2. Потом (в XX веке) были разработаны еще более сложные алгоритмы численного решения ОДУ — метод Рунге-Кутты с переменным шагом А, метод Adams, метод BDF, метод Radau и т. д.
|0<Ф) ^ http://www.wolframalpha.com
a"(t)=-a(t)*g/L, а(0)=а0, а'(0)=0 G
1
Input: |a"(t) = -a{t) x ^, a(0) = aO, a'(0) = o}
ODE classification: second-order linear ordinary differential equation
Alternate form: {cz"(t) + 8 0(0 - 0, aO - a(0), a'(0) - o}
Alternate form assuming aO, g, L and t are positive: (La"(t) -I- ga(t) = 0, aO = a(0), a'(0) = 0|
Differential equation solution: a( t) = aO cos| j
Рисунок 24. Аналитическое решение дифференциального уравнения второго порядка
Выбор метода решения ОДУ для конкретной задачи — это тонкая наука и высокое искусство. В функцию Odesolve, которую мы использовали в расчетах, дополнительно заложен алгоритм интерполяции дискретных значений, чтобы выдавалась не матрица дискретных значений, а именно функция (набор функций). Из-за этого иногда возникает недопонимание — некоторые считают, что если эта встроенная функция возвращает не таблицу, а искомую функцию, то получено аналитическое решение задачи.
Если же нужно решить не задачу с начальными условиями (задачу Коши), а краевую задачу, то применяют метод стрельбы. Краевая задача применительно к задаче о колебании маятника может быть сформу лирована так: известны значения
В. Ф. Очков, Е. П. Богомолова, Д. А. Иванов, К. Писачич
угла отклонения нити маятника от вертикали в два момента времени — определить положение нити в другие моменты времени. Метод стрельбы заключается в следующем. Задается недостающее начальное условие (угловая скорость нити, если иметь в виду задачу о маятнике).
Рисунок 25. Метод Рунге-Кутты, примененный для решения системы двух обыкновенных дифференциальных уравнений
Задача при этом сводится21 к известной задаче с начальными условиями, решая которую находят нужное значение на другом конце отрезка интегрирования. Если оно не совпадает с изначально заданным («промах»), то ведут «корректировку огня» — меняют в нужную сторону значение нужного параметра в начальной точке: скорость нити маятника или угол наклона пушки. Эта операция повторяется до тех
21 Как вскипятить воду в чайнике? Ответ. Нужно налить в него воды и поставить на огонь. Новая задача: в чайнике уже есть вода. Ответ математика. Воду нужно вылить. После этого задача сводится к предыдущей. Шутки шутками, но многие утром выливают воду из чайника, считая, что она за ночь испортилась, и наливают в чайник новую из крана или из фильтра.
пор, пока задача не будет решена с заданной точностью. Можно замерить положение спутника в два момента времени и по вышеописанной методике рассчитать его траекторию (орбиту — краевая задача). Можно это сделать также, замерив не только положение спутника, но и его скорость в какой-то момент времени (задача с начальными условиями). Если же нет влияния третьих космических тел, то можно замерить положение спутника в пяти точках движения (в разных моментах времени) и решить систему уже не дифференциальных, а линейных алгебраических уравнений (см. рис. 17с), имея в виду, что спутник движется строго по эллиптической (круговой), гиперболической или параболической (редкий случай) траектории.
Если же космических тел больше трех, то, повторяем, аналитических решений тут
22
нет .
На авторском сайте http://twtmas .mpei.ac.ru/mas/Worksheets/DE shot.mcd помещен интерактивный расчет методом стрельбы задачи о развитии эпидемии, сводящейся в простейшем варианте к решению системы двух обыкновенных дифференциальных уравнений первого порядка (http://twt.mpei.ac.ru/ochkov/mc8Pro.book/ 5 text.htm).
Часто можно слышать споры о том, какие пакеты численно или аналитически лучше решают системы ОДУ. Но все эти споры разбиваются о такой факт.
Математическое моделирование более-менее сложных (реальных) физических объектов часто сводится к решению, нескольких тысяч дифференциальных уравнений. И часто оказывается, что никакие разработанные и оттестированные средства, встроенные в математические пакеты Mathcad, Maple, Mathematica, ANSYS, COMSOL и др., не годятся для их решения. Тут приходится составлять индивидуальную расчетную схему со множеством допущений, упрощений, линеаризаций и других «хитростей». А это может сделать только специалист в данной области моделирования (космические полеты, энергетика, радиотехника и т. д., и т. п.) в паре с математиком.
В пакете Mathcad (в отличие от Maple и Mathematica) нет средств аналитического (символьного) решения дифференциальных уравнений. Из-за этого мы воспользовались сторонними средствами — см., например, рис. 24.
В этом (невозможность непосредственного аналитического решения ОДУ) часто упрекают пакет Mathcad, вернее, его разработчиков. Но тут нужно понимать, что более-менее сложное дифференциальное уравнение символьно не решается, а
22 Их математики и астрономы искали столетиями, пока не убедились, что их нет. Хотя — как считать. Некоторые задачи теплопроводности, сводящиеся к решению дифференциальных уравнений, не имели аналитических решений. Но потом математики придумали новые функции (функции Бесселя) и стали считать, что аналитические решения есть, но в них присутствуют эти самые функции Бесселя. По ним даже составили специальные таблицы, подобные таблицам синуса или логарифма. В принципе, для нашей задачи о движении спутника и двух планет (см. рис. 7) можно тоже придумать и затабу-лировать специальные функции, формирующие аналитическое решение этой задачи.
В. Ф. Очков, Е. П. Богомолова, Д. А. Иванов, К. Писачич
если и решается, то решение получается очень громоздким и неудобным для дальнейшего использования. В этом можно убедиться, если в исходное уравнение, записанное в окошке на рис. 24, вернуть синус, т. е. математический маятник превратить в физический.
На сайте https://www.ptcusercommunity.com/videos/1471 можно видеть анимацию (см. один ее кадр на рис. 26) численного решении ОДУ методами Эйлера (см. рис. 23) и Рунге-Кутты (см. рис. 25). На анимации видно, что «Эйлера слегка заносит на поворотах», а Рунге в паре с Куттой «четко чеканят шаг по начерченной линии». Аналогичное поведение зафиксировано и на рис. 23 и 25.
Рисунок 26. Кадр анимации сравнения метода Эйлера и метода Рунге-Кутты
На сайте https://www.ptcusercommunity.com/thread/85205 обсуждался процесс создания программ, показанных на рис. 23 и 25. Там же можно найти программу и пример анимации численного решения системы ОДУ методом Рунге-Кутты с переменным шагом (Adaptive — см. рис. 2), когда шаг меняется в зависимости от величины производных — правых частей ОДУ.
Заканчивая статью, можно привести еще одну цитату из А.С. Пушкина: «О сколько нам открытий чудных готовит просвещенья дух...», заменив в нем слово «просвещенья» на слово IT.
И самое последнее
Детство и отрочество первого автора статьи прошло в те времена, кода запускался первый спутник Земли, отправлялись зонды на Луну, был осуществлен первый полет человека космосе, первый выход человека в открытый космос в скафандре... Других особых успехов, успехов мирового уровня тогда у страны не было, В то время почти все мальчишки мечтали стать космонавтами. Автор не исключение. Но пойти учится на летчика, из которых набирались космонавты, ему не удалось — зрение подвело. Потом этот автор попал на срочную службу в советскую армию в ракетную часть, но опять же из-за зрения не в ракетчики, а в «обслугу» — был начальником котельной в воинской части, занимавшейся телеметрией ракетных пусков (http://twt.mpei .ac.ru/ochkov/Army.htm). Но теперь автор, имея под рукой Mathcad с его анимацией, может совершать виртуальные полеты в космос.
Литература
[1] Бутиков Е. И. Пакет компьютерных программ «Движение космических тел» [Электронный документ] http://butikov.faculty.ifmo.ru/Planets/index.html.
[2] Очков В. Ф. MAS на занятиях по математике, физике, информатике. // Компьютерные учебные программы и инновации. 2006. № 2. С. 73-80 (http://twt.mpei.ac.ru/ochkov/ Mathcad_12/Planet/index. html)
[3] Очков В. Ф. Живые кинематические схемы в Mathcad // Открытое образование. 2013. № 3. С. 27-33 (http://twt.mpei.ac.ru/ochkov/Mathcad-15/kinematic.html)
[4] Chenciner A., Montgomery R. A remarkable periodic solution of the three-body problem in the case of equal masses // Annals of Mathematics. 2000. Vol. 152. P. 881-901. (http://arxiv.org/abs/math/0011268)
[5] Очков В. Ф. Задачи по физике: новый подход к решению // Открытое образование. 2012. № 6 (http://twt.mpei.ac.ru/ochkov/Mathcad-15/Physic.pdf )
[6] Коробов В. И., Очков В. Ф. Химическая кинетика: введение с Mathcad/Maple/ MCS.— М. : Горячая линия-Телеком, 2009. (http://twt.mpei.ac.ru/TTHB/New-Chem-Kin/En-Ru-book.html)
[7] Очков В. Ф., Богомолова Е. П. Обратная задача коммивояжера или Необычные математические часы // Открытое образование. 2014. № 2. С. 22-28 (http://twt.mpei.ac.ru/ ochkov/SalesMan)
[8] Очков В. Ф., Богомолова Е. П. Интерполяция, экстраполяция, аппроксимация или «Ложь, наглая ложь и статистика» // Cloud of Science. 2015. T. 2. № 1. С. 61-88. (http://twt.mpei.ac.ru/ochkov/CoS_2_1.pdf)
[9] Очков В. Ф., Богомолова Е. П. Это страшное слово дифуры... // Информатика в школе. 2015. № 1. С. 55-58 (http://twt.mpei.ac.ru/ochkov/ODE.pdf)
В. Ф. Очков, Е. П. Богомолова, Д. А. Иванов, К. Писачич
[10] Ochkov V., Pisacic K. Journey from St. Petersburg to Moscow or Model of Gravity Train in Mathcad // Technical Journal University North Croatia, 2015. Vol. 9, No. 1. P. 1-5 (http://twt.mpei.ac.ru/ochkov/Tunnel.pdf)
Авторы:
Очков Валерий Федорович — доктор технических наук, профессор, профессор какфедры тепловых электрических станций Национального исследовательского университета «МЭИ» Богомолова Елена Петровна — кандидат физико-математических наук, доцент, доцент кафедры Высшей математики Национального исследовательского университета «МЭИ» Иванов Дмитрий Александрович — кандидат технических наук, доцент, доцент кафедры Общей физики и ядерного синтеза Национального исследовательского университета «МЭИ»
Писачич Катарина — инженер-исследователь, Северо-Хорвацкий университет, Хорватия
Movement of the Planets: the Calculation and Visualization in Mathcad or Kepler's Watch
V. F. Ochkov1, E. P. Bogomolova1, D. Ivanov1, K. Pisacic2
1National Research University Moscow "Power Engineering Institute", 14, Krasnokazarmennay street, Moscow, 111250
2University North, Mechanical engineering, Varazdin 104. brigade 3 42000 Varazdin, Croatia
e-mail: [email protected], [email protected]
Abstract. This article describes use of mathematical package Mathcad to solve motion of the planets (the material points) under the influence of gravitational forces on the plane and in three-dimensional space problem. The problem is projected onto the celestial mechanics. The questions of the limits of accuracy of numerical methods for solving ordinary differential equations is discussed. The concept of Kepler's clock with an uneven movement of the arrows is used to illustrate Kepler's second law. Keywords: Mathcad, Planets motions, Kepler Clock.
Reference
[1] Butikov E. I. Dvizhenie kosmicheskih tel. Available: http://butikov.faculty.ifmo.ru/ Planets/index.html. (In Rus)
[2] Ochkov V. F. (2006) Comuternye uchebnye programmy i innovacii, 2, 73-80. (In Rus)
[3] Ochkov V. F. (2013) Otkrytoe obrazovanie, 3, 27-33. (In Rus)
[4] Chenciner A., Montgomery R. (2000) Annals of Mathematics, 152, 881-901.
[5] Ochkov V. F. (2012) Otkrytoe obrazovanie, 6. (In Rus)
[6] Korobov V. I., Ochkov V. F. (2009) Himicheskaja kinetika: vvedenie s Mathcad/Maple/MCS. Moscow, Gorjachaja linija-Telekom, 2009. (In Rus)
[7] Ochkov V. F., Bogomolova E. P. (2014) Otkrytoe obrazovanie, 2, 22-28. (In Rus)
[8] Ochkov V. F., Bogomolova E. P. (2015) Cloud of Science, 2(1), 61-88. (In Rus)
[9] Ochkov V. F., Bogomolova E. P. Informatika v shkole, 1, 55-58 (In Rus)
[10] Ochkov V., Pisacic K. (2015) Technical Journal University North Croatia, 9(1), 1-5.
