Как решать уравнения в maple

REDMOND

Предисловие

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

С Maple (на кафедре была 6-я версия, а у лоточников домой была куплена 8-я) познакомился ещё студентом, когда начинал работать над будущей кандидатской под крылом моего первого (ныше покойного) научного руководителя. Были и добрые люди, что помогли на самом первом этапе разобраться с пакетом и начать работать.

И вот так постепенно на его плечи была переложена большая часть вычислительной работы по подготовке диссертации. Диссертация была защищена, а Maple навсегда остался надёжным помошником в научном труде. Часто бывает необходимо быстро оценить какую-нибудь задачу, составить уравнения, исследовать их аналитически, быстро получить численное решение, построить графики. В этом отношении Maple просто незаменим для меня (ни в коем разе не хочу обидеть приверженцев других пакетов).

Сделать всё то, что будет предложено читателю под катом, меня подвигла задача принесенная ученицей (приходится ещё заниматься и репетиторством) со школьной олимпиады. Условие задачи таково:

Если не придираться к некоторонной туманности условия, то задача достаточно проста, а её решение, полученное путем довольно громоздких для школьника выкладок, в общем виде дает результат

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

Это послужило катализатором для того, чтобы взять да и откопать свои старые задумки, накопленные ещё со времен работы в оргкомитете Всероссийской Олимпиады студентов по теоретической механике — три года подряд занимался там подготовкой задач компьютерного конкурса. Задумки касались автоматизации построения уравнений движений для механических систем с неудерживающими связями и трением, используя известные всем уравнения Лагранжа 2 рода

поборов стереотип многих преподавателей о том, что уравнения эти неприменимы к системам с неудерживающими связями и трением.

Что касается Maple, то его библиотека для решения задач вариационного исчисления дает возможность быстро получить уравнения Эйлера-Лагранжа, решение которых минимизирует действие по Гамильтону, что применимо для консервативных систем

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

Так как расматриваемые задачи не относятся к классу консервативных, то автором была предпринята попытка самостоятельно реализовать автоматизацию построения и анализа уравнений движений. Что из этого вышло, изложено под катом

1. Метод избыточных координат

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

Учет неудерживающих связей требует от нас определения и анализа величины их реакций, поэтому необходимо так же определить их величину. Уберем указанные связи и введем дополнительно r обобщенных координат, выразив через них кинетическую энергию системы

Составим s + r уравнений движения в форме уравнений Лагранжа 2 рода

содержащие s+r неизвестных координат и r неизвестных реакций связей. Считая связи удерживающими, дополняем данную систему уравнениями связей (для простоты рассматривая геометрические связи) в виде

получаем замкнутую систему уравнений, из которой находятся значения реакций

являющиеся функциями первых s (независимых) обобщенных координат и скоростей и они могут быть расчитаны на любом шаге интегрирования уравнений движения (1). Для удерживающих связей типа «нить/поверхность» уравнения (1) и (2) надо дополнить условием освобождения от связи

а для связей с сухим трением вида

где Fj и Nj соответственно касательная и нормальная составляющая реакции; vj — проекция скорости относительного проскальзывани точки приложения реакции.

Таким образом, уравнения (1) — (4) представляют собой полную математическую модель движения рассматриваемой механической системы.

Засим с теорией можно покончить и перейти к практике

2. Maple-функции построения и анализа уравнений Лагранжа

Для решения этой задачи была написана Maple-библиотека lagrange, содержащая четыре функции

LagrangeEQs — построение уравнений движения в форме Лагранжа 2 рода

В качестве входных параметров функция принимает выражение кинетической энергии T как функцию обобщенных координат и обобщенных скоростей; массив обобщенных координат q; массив радиус-векторов точек приложения сил r и массив векторов сил F.

LinksEQs — получение уравнений дифференциальных связей из уравнений геометрических связей

Здесь надо отметить, что система уравнений геометрических связей eqs должна содержать избыточные координаты в явном виде, то есть иметь вид

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

ReduceSystem — преобразование уравнений движения с учетом уравнений связей

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

SolveAccelsReacts — решение уравнений движения относительно реакций и обобщенных ускорений

Данная функция принимает на вход систему уравнений движения eqs, преобразованную с учетом уравнений связей. Она линейна относительно вторых производных независимых координат и реакций связей. Другие входные параметры: q — вектор независимых координат; R — массив реакций, относительно которых необходимо разрешить уравнения движения.

Теперь проиллюстрируем, как применять описанное «хозяйство» в деле

3. Задача о маятнике на тонкой нерастяжимой нити

Расчетная схема будет такой. В качестве обобщенной координаты выбираем угол наклона нити к вертикали.

Поскольку нить — неудерживающая связь, нас будет интересовать её реакция, а значит введем дополнительную, избыточную координату r(t).

Приступаем. Чистим память и подключаем библиотеку линейной алгебры

Подключаем библиотеку lagrange

Определяем вектор обобщенных координат, вычисляем координаты и скорость груза, а так же кинетическую энергию системы

На выходе получаем выражение для кинетической энергии (для вставки сюда использована функция latex(), генерирующая результат в LaTeX-нотации)

Формируем массив сил и массив координат точек их приложения

Скармливаем всё функции LagrangeEQs()

получая на выходе уравнения движения

Нетрудно убедится, что функция отработала нормально — для иллюстрации специально выбрана не слишком громоздкая задача.

Далее задаем уравнение связи — пока нить натянута, справедливо условие

REDMOND

преобразуем систему с учетом этого условия и находим реакцию связи

Сила натяжения нити равна

Система (5) — (7) является полной системой уравнений движения груза, с учетом возможности провисания нити. Теперь подготовим её к численному интегрированию. Для начала разрешим её относительно ускорений, передав в SolveAccelsReacts() уравнения (5) и (6), вектор обобщенных координат и пустой массив реакций

получая на выходе

Для численного моделирования, хоть это и не спортивно, напишем отдельный код, дабы не забивать голову читателя длительной обработкой полученной системы напильником. Тем более что моделирование будет иметь свои особенности.

Готовим исходные данные и систему уравнений движения

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

Теперь проверяем «школьное» решение задачи

В итоге, получаем результат, приведенный на скриншоте. Скорость груза в момент удара соответствует приведенному в предисловии значению, и видно, что до провисания нити груз движется по окружности, а после провисания нити движется как свободная точка под действием силы тяжести, по параболе.

Замечу, что погрешности попадания в гвоздь — вынужденная мера: в полярных координатах, которые были использованы, задача имеет особенность, понятную из уравнения (8). Поэтому r(t) сравнивалось не с нулем, а с величиной eps достаточно малой, чтобы получить решение, и достаточно большой, чтобы численный решатель fsolve() не сходил с ума. Однако это нисколько не умаляет практической ценности изложенных результатов.

Вместо заключения

Возможно, читатель упрекнет меня, что я стреляю из пушки по воробьям. Однако, хочется заметить, что всё сложное начинается с простого, а большая наука — с малых задач.

Научный форум dxdy

Условия для уравнения — как правильно задать в Maple

Здравствуйте. Ранее в Maple ставил изредка простые задачи для популярных уравнений типа теплопроводности или колебания — чисто, чтобы посмотреть, как это выглядит. Опирался в этом деле на Help самой программы (раз, два, три, четыре) и до сегодняшнего дня не парился — не особо умел, да и не было нужды. Поэтому попрошу навести прицел на ошибку или указать, что не так.

Есть уравнение в частных производных (pde) на функцию $q(X,T)$c мнимой единицей в качестве коэффициента.
Есть на плоскости $(X,T)$кривая с началом в $(0,0)$, допустим, полукубическая парабола $X(T)=\pm a (-b T)^<3/2>$, где константы $a, b$известны. Полагается рассмотрение отрицательных $T$.

Условия такие:
(bc2). От левой ветви (не включительно) этой кривой до $X=0$(включительно) функция равна некоторой функции $q(X,T^<*>)=q_1 (X,T^<*>)$, от $X=0$(включительно) до правой ветви (не включительно) — другой функции $q(X,T^<*>)=q_3 (X,T^<*>)$.
(bc1). Левее левой ветви и правее правой полагаем, что $q=0$.

Для кривой $X(T)=\pm a (-b T)^<3/2>$ рассматриваем разные значения $T=T^<*>$, соответственно меняем и краевые условия, подставляя в них выбранное $T=T^<*>$. То есть, одну и ту же задачу решаем несколько раз для разных значений $T=T^<*>$.

И вот, что я скармливаю Maple, выбрав $T^<*>=-9$ и подставив его в уравнение кривой и в нужные функции. Подразумевается, что вычислять надо не всюду, а для конечного отрезка $X$от $-1000$до $1000$.

На что получаем ответ:

То есть, условия надо определить для каждой независимой переменной функции $q(X,T)$в одной или двух точках. Но у меня нет больше условий. Ни на производные, ни на $X$, а $T$и так в деле.
Может, взять функции $q_1, q_3$, подставить $T^<*>=-9$ в их производные, дополнить условия? Но это выглядит излишним.

Полагаю, ошибка в самой команде pdsolve или в принципиальной постановке задачи. Буду признателен, если поможете разобраться. К сожалению, сам пока не понял.

Последний раз редактировалось svv 28.05.2019, 19:37, всего редактировалось 2 раз(а).

Попробуйте в порядке эксперимента добавить два дополнительных условия:
$q(-1000, T)=0$
$q(+1000, T)=0$
Просто чтобы понять, верно ли, что Maple не хватает boundary conditions. Все условия, которые Вы до сих пор задавали — initial conditions, хотя Вы их обозначали bc1, bc2.

Последний раз редактировалось GAA 29.05.2019, 16:04, всего редактировалось 1 раз.
Исправлена опечатка: "састных" на "частных"

Попробуйте в порядке эксперимента добавить два дополнительных условия:
$q(-1000, T)=0$
$q(+1000, T)=0$
Просто чтобы понять, верно ли, что Maple не хватает boundary conditions. Все условия, которые Вы до сих пор задавали — initial conditions, хотя Вы их обозначали bc1, bc2.

Честно говоря, тут не понял, что не устраивает Мэйпл. И на $X$условие, и на $T$.

НУШ
$-i \varepsilon q_T=\varepsilon^2 q_<XX>-2 \left | q \right |^2 q, \ (\varepsilon \ll 1)$
с малым параметром $\frac<1><100>$:
$-i (\frac<1><100>) q_T=(\frac<1><100>)^2 q_<XX>-2 \left | q \right |^2 q$.
С комплексными коэффициентами могут быть сложности?

Решение множества дифференциальных уравнений в Maple с граничными условиями

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

Я попытался воспроизвести свою проблему в простом случае, который выглядит следующим образом:

Таким образом, они отвязаны на данный момент и довольно просты в решении. Тем не менее, я застрял с Кленом по какой-то причине, чего я действительно не понимаю.

Для следующих ситуаций (как видно из вышеизложенного), Maple дает мне решение

    Решение двух уравнений без граничных условий. Решение первого уравнения с граничным условием.

Теперь, когда я совмещаю оба, я получаю сообщение об ошибке

Ошибка (в dsolve) находила следующие уравнения, не зависящие от неизвестных входной системы:

Если я явно определяю P и Q как неизвестные (которые я взял в качестве подсказки) (см. Последнюю строку кода)

Ошибка (в dsolve) получила указание как неизвестную проблему и обнаружила связанную функцию (s) , не зависящую от . Укажите функциональность неизвестного явно

Заметки

    Если я определяю Q как функцию y только в e2, то операторы dsolve работают, но это не то, что я хочу Я заглянул в пдсольве, но пока не повезло. Т.е. я пробовал pdsolve(); Однако это не дает сообщений об ошибках, но также не выводит. Я несколько раз искал сообщения об ошибках, которые я получил, но это еще не привело меня к сожалению. Если y задает P в e1 как функцию от x и y P(x,y) и задает cond:=P(x,0)=0 , то dsolve не дает никакого вывода.

То, что я пытаюсь решить

Фактическая проблема, которую я пытаюсь решить, заключается в следующем.

Как этот набор можно решить в Maple в одной команде? Я знаю, что ответ должен быть

Любые советы/советы очень ценятся!

Проблема здесь полностью процедурная и математическая.

Рассмотрим первую ошибку:

Ошибка (в dsolve) находила следующие уравнения, не зависящие от неизвестных входной системы:

Это говорит вам, что у вас есть постоянная связь. Буквально первая производная от вашей функции Q(x,y) является некоторой константой C и далее пытается найти следующую производную бомбу, потому что вы пошли так низко, как можете. Чтобы быть полностью педантичным, оценка должна быть нулевой, но это не похоже на то, что Maple поддерживает.

Второе сообщение об ошибке говорит вам, что когда вы определяете P и Q как ваши неизвестные, у вас все еще нет вашего истинного неизвестного, y , связанного с решателем. Это также приводит к трудностям Maple.

Я думаю, что если бы мы увидели ваше определение для Q(x,y) , мы бы очень быстро обнаружили, что, по крайней мере, согласно Maple, функция фактически не зависит от какого-либо значения x . Это то, что вызывает часть ваших трудностей.

Кроме того, я не понимаю, почему вы называете diff вообще, поскольку принятие производной и второй производной от двух функций не похоже на то, что вы хотите сделать. Возможно, их следует заменить dsolve , с вашими условиями, указанными здесь?

REDMOND

Добавить комментарий

Ваш адрес email не будет опубликован. Обязательные поля помечены *