Дифференциально-аналитическое уравнение

Примечание

Для анимации в расчет введена переменная frame, имеющая на рис. 5.14 нулевое значение по умолчанию. Если значение переменной frame будет равно 99, то значение переменной t станет равным значению переменной tend — интегрирование будет завершено, а анимация полностью построена.

На рис. 5.20 показано развитие задачи о маятнике— ввод в модель еще одной силы — силы сопротивления среды (воздуха, в частности) с выкладыванием решения в Интернете.

В решении задачи о маятнике, показанном на рис. 5.20, эта сила пропорциональна плотности среды ра, перемноженной на площадь сечения груза f (принимается, что это шар) и скорость груза. В реальной жизни эта сила пропорциональна скорости не в первой степени, а во второй, вернее, эта степень также зависит от скорости по "хитрому" закону — ламинарное обтекание при малых скоростях и турбулентное при больших, влияние "хвоста" вихрей воздуха от предыдущего "взмаха" маятника и т. д. Но простая замена на рис. 5.20 х' (t) и у' (t) на х' (t)2 и у' (t)2 приведет к потере знака вектора скорости и к ошибке в решении задачи.

Это явление мы проиллюстрируем в другой популярной физико-математической задаче — задаче об остановке автомобиля, показанной на рис. 5.21.

На рис. 5.21 между ключевым словом Given и функцией odesoдve зажато только одно дифференциальное уравнение второго порядка, описывающее баланс сил, действующих на останавливающийся автомобиль: сила инерции (ускорение, помноженное на массу), сила трения качения (коэффициент трения, перемноженный на вес автомобиля) и сила сопротивления воздуха, зависящая от скорости в квадрате (см. выше). Здесь, в отличие от задачи о маятнике, скорость (первая производная) не меняет своего знака, а только изменяется от начального заданного значения до нулевого (остановка автомобиля). Проблемы возникают с силой трения. Если просто записать в дифференциальном уравнении остановки автомобиля, что сила трения всегда постоянна и равна весу автомобиля, перемноженному на коэффициент трения качения, то автомобиль после остановки... покатится назад. Дело в том, что сила трения равна вышеотмеченному произведению только при ненулевой скорости. Поэтому в дифференциальное уравнение на рис. 5.21 была введена встроенная в Mathcad функция sign, возвращающая 0 при нулевом своем аргументе либо 1 или -1 в противном случае. Теперь это наш автомобиль перестает катиться после остановки.

Рис. 5.21. Моделирование остановки автомобиля

Функцию sign автор пытался вставить и в дифференциальные уравнения, описывающие затухающие колебания маятника (см. рис. 5.20). Это была попытка заменить первую степень скорости в силе сопротивления воздуха на вторую. Но это срывало численное решение задачи — "угловатость" функции типа if, sign и др. не "дружит" не только с символьной, но и с численной математикой — природа не терпит острых углов. А сила трения — это довольно-таки "рваная" функция, особенно в точках старта и останова, способная загнать в тупик любого решателя — человека или машину.

Особо следует коснуться проблемы единиц измерения при решении дифференциальных уравнений, начатой в разд. 1.4.

С одной стороны, дифференциальные уравнения, как правило, описывают реальные физические явления и, следовательно, все переменные должны хранить не просто величины, а физические величины — величины с размерностями длины, времени, силы и т. д.

Но с другой стороны, встроенные в Mathcad функции решения дифференциальных уравнений типа rkfixed, odesolve и др., описанные в данной главе, не работают с размерными величинами. Это как раз те функции, которые препятствуют полновесному использованию единиц измерения в Mathcad-расчетах — проблема, поднятая разд. 1.4.

Разрешение этого противоречия показано на рис. 5.21. При решении задачи об остановке автомобиля все переменные были размерными, но в самом дифференциальном уравнении и в начальных условиях эти переменные лишались своих размерностей делением на соответствующие единицы измерения. Потом, после решения нужная единица измерения (метр) приписывалась к ответу.

Примечание

В задачах о маятнике (см. рис. 5.18 и 5.20) все переменные были безразмерными и хранили величины в основных единицах СИ (м, кг, Н и т. д.), что чревато ошибками в формулах.

Вдумчивый читатель наверно отметил тот факт, что автор не применил функцию odesolve к решению краевой проблемы в задаче об эпидемии, а стал использовать довольно сложные методы, показанные на рис. 5.9— 5.12. Дело в том, что функция odesolve решает краевую задачу на дифференциальном уравнении второго порядка, но не решает краевую задачу на двух уравнениях первого порядка (наша задача об эпидемии, например), хотя, с точки зрения математики, это — одинаковые задачи. Так на рис. 5.22 через решение краевой задачи на дифференциальном уравнении второго порядка и особые "хитрости" форматирования трехмерного графика17 нарисован висячий мост, цепи которого провисают и дают кривую, вторая производная которой равна самой функции (свойство так называвмой цепной функции). Но это "изящное" решение нельзя перенести на краевую задачу развития эпидемии.

Рис. 5.22. Краевая задача и функция odesolve

Кстати, о жесткости/нежесткости систем дифференциальных уравнений. Автор долго искал несложный пример, где бы применение "нежестких" или "жестких" алгоритмов решения ОДУ давало бы не количественную (одна кривая пошла чуть круче другой, полученной при использовании иного метода интегрирования — см. рис. 5.6, например), а качественную (кривая прошла бы совсем иначе) картину численного интегрирования. В конце концов, автор такой пример сконструировал сам — это расчет траектории движения трех небесных тел, подчиняющихся законам Ньютона — второму закону механики и закону всемирного тяготения (http://twt.mpei.ac.ru/MAS/ Worksheets/3_PIanets.mcd). Так вот, в старых версиях Mathcad при определенных начальных условиях (положение планет и их скорости в нулевой момент времени) и при ориентации на "жесткий" алгоритм удалось смоделировать довольно интересный случай перехвата одной планетой спутника у другой планеты (рис. 5.23).

Рис. 5.23. Моделирование перехвата спутника

В среде же Mathcad 14 с его гибридным (Adams/BDF) алгоритмом решения системы ОДУ первая планета стала не перехватывать спутник у другой планеты, а просто сбивать его с орбиты, и автору никак не удается скорректировать начальные условия, чтобы восстановить модель перехвата спутника. Но это так— лирическое отступление... Главное же в том, что изменение и дополнение в инструментах решение ОДУ в среде Mathcad не очень, так скажем, волнуют пользователей. Это подтверждается и коллекцией соответствующих файлов, хранящихся в Интернете, где решение ОДУ не дополняется указаниями о подстройке встроенных инструментов. На форуме, посвященном Mathcad 14, эти новинки также особо не обсуждаются. Дело в том, что тот, кто разбирается в этих тонкостях и решает реальные сложные задачи, либо составляет свои собственные алгоритмы решения, задействуя инструменты программирования Mathcad (циклы, альтернативы и др. — см. несложные примеры на рис. 5.5, где встроенная в Mathcad функция rkfixed открыта для изучения) или классические языки программирования, либо прибегает к специализированным пакетам решения подобных задач, не полагаясь на "закрытые" алгоритмы. Такая же ситуация, но более четко, вырисовывается при решении систем алгебраических уравнений блоком Given/Find— если система не решается или решается неверно (см., например, рис. 2.33), то смена настроек функции Find (см. рис. 2.22), как правило, не исправляет ситуацию. Тут приходится варьировать первым приближением или создавать индивидуальную программу решения, основанную на учете физической сущности каждого уравнения, входящего в систему, и на применении к различным уравнениям различных критериев точности ("жесткость" систем алгебраических уравнений).