Дифференциальные уравнения с неявно заданным x

Калькулятор Обыкновенных Дифференциальных Уравнений (ОДУ) и Систем (СОДУ)

Порядок производной указывается штрихами — y»’ или числом после одного штриха — y’5

Ввод распознает различные синонимы функций, как asin , arsin , arcsin

Знак умножения и скобки расставляются дополнительно — запись 2sinx сходна 2*sin(x)

Список математических функций и констант :

• ln(x) — натуральный логарифм

• sh(x) — гиперболический синус

• ch(x) — гиперболический косинус

• th(x) — гиперболический тангенс

• cth(x) — гиперболический котангенс

• sch(x) — гиперболический секанс

• csch(x) — гиперболический косеканс

• arsh(x) — обратный гиперболический синус

• arch(x) — обратный гиперболический косинус

• arth(x) — обратный гиперболический тангенс

• arcth(x) — обратный гиперболический котангенс

• arsch(x) — обратный гиперболический секанс

• arcsch(x) — обратный гиперболический косеканс

Неявные методы решения дифференциальных уравнений

1. Неявные методы

Методы численного решения ОДУ, с которыми мы познакомились в первом разделе этого курса (метод Эйлера, метод средней точки и т. п.), называются «явными» методами. Однако иногда система ОДУ может стать «жесткой», а решать такие системы явными методами неудобно. В этом случае желательно изменить формулировку задачи так, что не пришлось иметь дело с жесткой системой. Но это не всегда возможно, поэтому вы должны уметь решать жесткие ОДУ. Для этого, как правило, используют методы решения, которые называются «неявными».

2. Пример жесткого ОДУ

Во-первых, в чем смысл и причина появления жестких уравнений? Давайте рассмотрим пример, который часто возникает в динамике. Предположим, что у нас есть частица, с координатами \((x(t),y(t))\) , и предположим, что мы хотим, чтобы ее \(y\) -координата всегда оставалась равной нулю. Один из способов добиться этого — добавить слагаемое \(-ky(t)\) , к производной \(\dot(t)\) , где \(k\) — большая положительная постоянная. Если \(k\) достаточно велико, то частица никогда не уйдет далеко от \(y(t)=0\) , так как слагаемое \(-ky(t)\) всегда приведет \(y(t)\) обратно к нулю. Предположим далее, что мы хотим, чтобы пользователь мог перемещать частицу как угодно вдоль оси \(x\) . Дифференциальное уравнение, описывающее движение данной системы, будет иметь вид

(Кроме того, мы предполагаем, что частица не запускается из \(y_0=0\) ). В результате частица будет сильно притягиваться к прямой \(y = 0\) , и менее сильно — к \(x = 0\) . Если решать ОДУ на достаточно продолжительном интервале времени, то частица рано или поздно попадет в точку \((0, 0)\) и останется в ней.

Теперь предположим, что для решения уравнения мы используем метод Эйлера. Если сделать шаг размера \(h\) , то получим

Если мы посмотрим на \(у\) -компоненту этого уравнения, то увидим, что при \(|1-hk|>1\) , вычисленное нами \(y_\) будет по модулю больше, чем \(|у_0|\) . Другими словами, при \(|1-hk|>1\) метод Эйлера будет неустойчив: каждый шаг приводит к увеличению \(y_\) по сравнению с предыдущим значением и приближенное решение будут все дальше отклоняться от нуля. Таким образом, для обеспечения устойчивости метода нужно, чтобы \(1-hk>-1\) или \(hk . Самый большой шаг, который мы можем сделать не нарушив устойчивости, должен быть меньше \(2/k\) .

Теперь, если \(k\) – большое число, мы вынуждены делать очень маленькие шаги. Это означает, что частица будет двигаться к \((0,0)\) мучительно медленно. Даже если взять \(y_0\) очень близким к нулю, то придется делать настолько маленькие шаги, что изменение \(x\) -координаты будет практически незаметно. Вот так выглядит жесткое ОДУ. В данном случае жесткость возникает из-за слишком большого \(k\) , призванного удержать частицу возле прямой \(у = 0\) . Позже, когда мы будем рассматривать частицы, соединенные пружинами, мы увидим то же самое явление: от жесткости пружины и происходит термин «жесткое» ОДУ. Даже если мы используем более совершенный численный метод, такой как метод Рунге-Кутта 4-го порядка, это лишь слегка улучшит ситуацию с выбором величины шага, но мы все равно будем иметь серьезные вычислительные проблемы.

Теперь, как мы уже говорили выше, нужно попытаться переформулировать свою задачу так, чтобы избежать появления жесткого ОДУ. Если же это не получится, то нужно использовать неявный метод решения ОДУ. Метод, который мы покажем ниже, является самым простым из неявных методов. Он основан на том, что шаг обычного метода Эйлера выполняется «наоборот».

3. Решение жесткого ОДУ

Пусть дано дифференциальное уравнение

Формула явного метода Эйлера

продвигает систему вперед на шаг \(h\) во времени. Для жестких систем, однако, удобнее заменить эту формулу на следующую

То есть, нам нужно вычислить \(f\) в точке, в которую мы стремимся попасть ( \(\textbf_\) ), а не в исходной ( \(\textbf_0\) ). (Если представить, что время может двигаться в обратном направление, то смысл этого уравнения очевиден. Оно говорит: «если вы находитесь в \(\textbf_\) , и сделаете шаг \(-hf (\textbf_)\) , то попадете в \(\textbf_0\) ». Так что если ваше ДУ представляет собой систему, движущуюся вспять во времени, то этот шаг имеет смысл. Это просто поиск точки \(\textbf_\) такой, что если запустить время вспять, вы в конечном итоге окажетесь в \(\textbf_0\) .) Таким образом, мы ищем точку \(\textbf_\) такую, что \(f\) , вычисленная в ней и умноженная на \(h\) , приводит к исходной точке \(\textbf_0\) . К сожалению, найти \(\textbf_\) из уравнения \eqref в общем случае невозможно, если только \(f\) не является линейной функцией.

Чтобы справиться с этим, заменим \(f (\textbf_)\) линейной аппроксимацией, основанной на разложении \(f\) в ряд Тейлора. Введем обозначение \(\Delta\textbf=\textbf_-\textbf_0\) . Подставив его в уравнение \eqref, получим

Теперь заменим \(f(\textbf_0 + \Delta\textbf)\) следующим приближением

(Заметим, что поскольку \(f(\textbf_0)\) является вектором, то производная \(f^\prime(\textbf_0)\) является матрицей.) Используя это приближение, мы можем записать \(\Delta\textbf\) как

Разделим обе части последнего соотношения на \(h\) и перепишем результат в виде

где \(I\) — единичная матрица.

Разрешая это соотношение относительно \(\Delta\textbf\) , получим

Вычисление \(\textbf_=\textbf_0+\Delta\textbf\) для неявного метода очевидно требует больших вычислительных затрат, чем при использовании явного метода, так как мы должны решать систему линейных уравнений на каждом шаге. Хотя это может показаться серьезным недостатком (в вычислительном плане), не отчаивайтесь (пока). Для многих классов задач, матрица \(f^\prime\) будет разреженной — например, для «решетки» частиц, соединенных пружинами, \(f^\prime\) будет иметь структуру, которая соответствует связям между частицами. (Обсуждение разреженности и методов решений, применимых в этом случае, см. Baraff и Witkin [1]. Основной материал в Press et al. [2] также будет полезен.) В результате, как правило, можно решить уравнение \eqref в линейном времени (т. е. за время, пропорциональное размерности \(\textbf\) ). Выигрыш в таких случаях будет весьма существенным: мы, как правило, можем делать большие шаги по времени без потери устойчивости (т. е. без расхождения, как это происходит для явного метода, если длина шага слишком велика). Время, необходимое для решения каждой линейной системы, таким образом, более чем компенсируется тем, шагов можно зачастую сделать на порядки больше, чем при использовании явных методов. (Конечно, код, необходимый для реализации всего этого, гораздо сложнее, чем в явном случае; как мы уже говорили: переделывайте свои задачи в нежесткие, а если не получается, то платите положенную цену.)

Применим теперь неявный метод для решения уравнения \eqref. В нашем случае \(f(\textbf(t))\) равно

Дифференцирование по \(\textbf\) дает

Тогда матрица \((1/h)\textbf — f^<\prime>(\textbf_0)\) будет равна

Обращая эту матрицу и умножая на \(f(\textbf_0)\) , получим

Какова же предельная длина шага в этом случае? Ответ: нет предела! Если позволить \(h\) расти до бесконечности, мы получаем следующее

Это означает, что мы достигнем \(\textbf_=\textbf_0 + (-\textbf_0)=0\) за один шаг! В общем случае для жесткого ОДУ мы не сможем сделать шаг произвольного размера, но мы сможем сделать его гораздо большим, чем если бы использовали явный метод. Дополнительные расходы на решение системы линейных уравнений компенсируются экономией, возникающий благодаря возможности сделать меньше шагов.

4. Решение уравнений второго порядка

Большинство задач динамики записывается в виде ДУ 2-го порядка

Это уравнение легко преобразуется систему ДУ 1-го порядка, добавлением новых переменных. Если мы определим \(\textbf=\dot<\textbf>\) , то сможем переписать уравнение \eqref в виде

что представляет собой систему ДУ 1-го порядка. Однако, применяя обратный (неявный) метод Эйлера к уравнению \eqref, получим линейную систему размера \(2n \times 2n\) где \(n\) — размерность \(\textbf\) . Простое преобразование позволяет уменьшить размер линейной системы до \(n \times n\) . Важно отметить, что обе системы — исходная и преобразованная — будут иметь одинаковую степень разреженности. Таким образом, решение системы меньшего размера будет выполняться быстрее.

Система \(n \times n\) , которая должна быть решена, получается следующим образом. Введем для краткости следующие обозначения \(\textbf_0=\textbf(t_0)\) и \(\textbf_0=\textbf(t_0)\) . Определяем \(\Delta\textbf\) и \(\Delta\textbf\) как \(\Delta\textbf = \textbf(t_0+h)-\textbf(t_0)\) и \(\Delta\textbf = \textbf(t_0+h)-\textbf(t_0)\) . Очередное приближение по неявному методу Эйлера, примененному к уравнению \eqref, дает в результате

Применяя к \(f\) разложение в ряд Тейлора, которое в данном контексте является функцией и \(\textbf\) и \(\textbf\) , получим приближение 1-го порядка

В этом уравнении, производная \(\partial f / \partial\textbf\) оценивается в точке \((\textbf_0, \textbf_0)\) и аналогично вычисляется \(\partial f / \partial\textbf\) . Подставляя это приближение в уравнение \eqref, получим линейную систему

Подставив во вторую строку последнего равенства соотношение \(\Delta\textbf = h(\textbf_0 + \Delta\textbf)\) (т. е. первую строку того же равенства), придем к

Вводя единичную матрицу \(I\) и перегруппировав члены, получим соотношение

из которого находим \(\Delta\textbf\) . А зная \(\Delta\textbf\) легко вычислить \(\Delta\textbf = h(\textbf_0 + \Delta\textbf)\) .

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

Ссылки

1. D. Baraff and A. Witkin. Large steps in cloth simulation. Computer Graphics (Proc. SIGGRAPH), 1998.
2. W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling. Numerical Recipes. Cambridge University Press, 1986.

Читайте также

Комментарии

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

Дифференциальные уравнения высших порядков, не содержащие переменную в явном виде

Метод решения

Рассмотрим уравнение, не содержащее независимую переменную в явном виде:
(1) .
Порядок этого уравнения понижается на единицу с помощью подстановки:

Далее считаем, что функция u зависит от переменной y , тогда:
;
;
и т. д.

В результате такой подстановки, порядок уравнения понижается на единицу.

Пример

Уравнение не содержит независимую переменную в явном виде. Делаем подстановку:
.
Считаем, что функция u зависит от переменной y . Тогда
.

Подставляем в исходное уравнение:
.
Делим на u . При имеем:
.
Это уравнение с разделяющимися переменными. Делим на и умножаем на dy . При имеем:
.
Интегрируем:
(2) .

Подставляем в (2):
.
Потенцируем:
.
Заменим постоянную интегрирования . Знак модуля сводится к умножению на ±1 . Включим ±1 в постоянную . То есть мы теперь полагаем, что может быть не только положительным, но и отрицательным числом. Тогда:
.

Выполняем преобразования:
;
.
При имеем:
;
.

Разделяем переменные:
.
Интегрируем:
(3) .

Вычисляем интеграл:

.
Подставляем в (3):
;
.
Возводим в квадрат и выполняем преобразования:
;
;
(4) .

При выводе формулы (4) мы предполагали, что
и .
Теперь рассмотрим случаи
.
Нетрудно видеть, что решение, охватывающее эти три равенства, есть
(5) ,
где C – произвольная постоянная. Тогда . Подставляя это в исходное уравнение нетрудно убедиться, что оно выполняется. Это особое решение. Добавим его в ответ.

Автор: Олег Одинцов . Опубликовано: 18-07-2013 Изменено: 27-06-2018