Численное решение обыкновенных дифференциальных уравнений (ОДУ) в Python
Рассмотрены приемы решения обыкновенных дифференциальных уравнений (ОДУ) с помощью модуля scipy.integrate языка Python
Краткое описание модуля scipy.integrate
Модуль scipy.integrate имеет две функции ode() и odeint(), которые предназначены для решения систем обыкновенных дифференциальных уравнений (ОДУ) первого порядка с начальными условиями в одной точке (т.е. задача Коши).
Функция ode() более универсальная, а функция odeint() (ODE integrator) имеет более простой интерфейс и хорошо решает большинство задач.
Функция odeint() имеет три обязательных аргумента и много опций. Она имеет следующий формат
Решение одного ОДУ
Допустим надо решить диф. уравнение 1-го порядка
Получилось что-то такое:
Решение системы ОДУ
Пусть теперь мы хотим решить (автономную) систему диф. уравнений 1-го порядка
Выходной массив w состоит из двух столбцов — y1(t) и y2(t).
Также без труда можно построить фазовые траектории:
Задачи с начальными условиями для систем обыкновенных дифференциальных уравнений
Рассмотрим задачу Коши для системы обыкновенных дифференциальных уравнений $$ \begin
Используя векторные обозначения, задачу (1), (2) можно записать как задачу Коши $$ \begin
Численные методы решения задачи Коши
Существует большое количество методов численного решения задачи (3), (4). Вначале рассмотрим простейший явный метод Эйлера и его программную реализацию. Затем будут представлены методы Рунге—Кутта и многошаговые методы.
При построении численных алгоритмов будем считать, что решение этой дифференциальной задачи существует, оно единственно и обладает необходимыми свойствами гладкости.
Идея численных методов решения задачи (3), (4) состоит из четырех частей:
1. Вводится расчетная сетка по переменной \( t \) (время) из \( N_t + 1 \) точки \( t_0 \), \( t_1 \), \( \ldots \), \( t_
2. Предполагаем, что дифференциальное уравнение выполнено в узлах сетки.
3. Аппроксимируем производные конечными разностями.
4. Формулируем алгоритм, который вычисляет новые значения \( \pmb
Явный метод Эйлера
Проиллюстрируем указанные шаги. Для начала введем расчетную сетку. Очень часто сетка является равномерной, т.е. имеет одинаковое расстояние между узлами \( t_n \) и \( t_
Затем, предполагаем, что уравнение выполнено в узлах сетки, т.е.: $$ \pmb^\prime (t_n) = \pmb
Заменяем производные конечными разностями. С этой целью, нам нужно знать конкретные формулы, как производные могут быть аппроксимированы конечными разностями. Простейший подход заключается в использовании определения производной: $$ \pmb^\prime(t) = \lim_ <\tau \to 0>\frac<\pmb(t+\tau) — \pmb(t)><\tau>. $$
В произвольном узле сетки \( t_n \) это определение можно переписать в виде: $$ \begin
Четвертый шаг заключается в получении численного алгоритма. Из (5) следует, что мы должны знать значение \( y^n \) для того, чтобы решить уравнение (5) относительно \( y^
При условии, что у нас известно начальное значение \( \pmb
Программная реализация явного метода Эйлера
Выражение (6) может быть как скалярным так и векторным уравнением. И в скалярном и в векторном случае на языке Python его можно реализовать следующим образом
При решении системы (векторный случай), u[n] — одномерный массив numpy длины \( m+1 \) (\( m \) — размерность задачи), а функция F должна возвращать numpy -массив размерности \( m+1 \), t[n] — значение в момент времени \( t_n \).
Таким образом численное решение на отрезке \( [0, T] \) должно быть представлено двумерным массивом, инициализируемым нулями u = np.zeros((N_t+1, m+1)) . Первый индекс соответствует временному слою, а второй компоненте вектора решения на соответствующем временном слое. Использование только одного индекса, u[n] или, что то же самое, u[n, :] , соответствует всем компонентам вектора решения.
Функция euler решения системы уравнений реализована в файле euler.py:
Строка F_ = lambda . требует пояснений. Для пользователя, решающего систему ОДУ, удобно задавать функцию правой части в виде списка компонент. Можно, конечно, требовать чтобы пользователь возвращал из функции массив numpy , но очень легко осуществлять преобразование в самой функции решателе. Чтобы быть уверенным, что результат F будет нужным массивом, который можно использовать в векторных вычислениях, мы вводим новую функцию F_ , которая вызывает пользовательскую функцию F «прогоняет» результат через функцию assaray модуля numpy .
Неявный метод Эйлера
При построении неявного метода Эйлера значение функции \( F \) берется на новом временном слое, т.е. для решении задачи (5) используется следующий метод: $$ \begin
Таким образом для нахождения приближенного значения искомой функции на новом временном слое \( t_
Для решения уравнения (8) можно использовать, например, метод Ньютона.
Программная реализация неявного метода Эйлера
Функция backward_euler решения системы уравнений реализована в файле euler.py:
Отметим, что для нахождения значения u[n+1] используется функция fsolve модуля optimize библиотеки scipy . В качестве начального приближения для решения нелинейного уравнения используется значение искомой функции с предыдущего слоя u[n] .
Методы Рунге—Кутта
Одношаговый метод Рунге—Кутта в общем виде записывается следующим образом: $$ \begin
Одним из наиболее распространенных является явный метод Рунге-Кутта четвертого порядка: $$ \begin
Многошаговые методы
В методах Рунге—Кутта в вычислениях участвуют значения приближенного решения только в двух соседних узлах \( \pmb
Различные варианты многошаговых методов (методы Адамса) решения задачи с начальными условиями для систем обыкновенных дифференциальных уравнений могут быть получены на основе использования квадратурных формул для правой части равенства $$ \begin
Для получения неявного многошагового метода используем для подынтегральной функции интерполяционную формулу по значениям функции \( \pmb
Для интерполяционного метода Адамса (15) наивысший порядок аппроксимации равен \( m+1 \).
Для построения явных многошаговых методов можно использовать процедуру экстраполяции подынтегральной функции в правой части (14). В этом случае приближение осуществляется по значениям \( \pmb
Для экстраполяционного метода Адамса (16) погрешность аппроксимации имеет \( m \)-ый порядок.
На основе методов Адамса строятся и схемы предиктор–корректор. На этапе предиктор используется явный метод Адамса, на этапе корректора — аналог неявного метода Адамса. Например, при использовании методов третьего порядка аппроксимации в соответствии с (18) для предсказания решения положим $$ \frac<\pmb
Жесткие системы ОДУ
При численном решении задачи Коши для систем обыкновенных дифференциальных уравнений (3), (4) могут возникнуть дополнительные трудности, порожденные жесткостью системы. Локальные особенности поведения решения в точке \( u = w \) передаются линейной системой $$ \begin
Пусть \( \lambda_i(t) \), \( i = 1, 2, \ldots, m \) — собственные числа матрицы $$ \begin
Для численное решения жестких задач используются вычислительные алгоритмы, которые имеют повышенный запас устойчивости. Необходимо ориентироваться на использование \( A \)-устойчивых или \( A(\alpha) \)-устойчивых методов.
Метод называется \( A \)-устойчивым, если при решении задачи Коши для системы (3) область его устойчивости содержит угол $$ \begin
5++ способов в одну строку на Python решить первую задачу Проекта Эйлера
Однажды меня посетила мысль, а что если попробовать решить первую задачу Проекта Эйлера всевозможными способами, но с условием, что решение должно быть в одну строку. В итоге получилось более пяти однострочных решений с применением Filter, Map, Reduce, Generator Expression и т.д. В этой статье я покажу то, к чему я пришёл.
Это моя первая статья. Стоит отнестись к ней настороженно. Уникальные решения будут оформлены в отдельные пункты. Менее уникальные — в подпункты.
Условие задачи
Если выписать все натуральные числа меньше 10, кратные 3 или 5, то получим 3, 5, 6 и 9. Сумма этих чисел равна 23.
Найдите сумму всех чисел меньше 1000, кратных 3 или 5.
00 — Базовое решение
Прежде чем перейти непосредственно к однострочным решениям, разумно было бы упомянуть сначала стандартное, классическое решение:
Перебираем последовательность чисел от 1 до 999. Если перебираемое число делится на 3 или на 5 без остатка от деления, то прибавляем каждое такое число в заранее объявленную переменную result .
01 — Generator Expression. Выражение-генератор
Числа из последовательности от 1 до 999, делящиеся на 3 или на 5 без остатка от деления, собираются в генератор. Затем функция sum() складывает содержимое генератора.
01.a — List Comprehension. Выражение на основе списка
В отличии от предыдущего, здесь выражение дополнительно помещается в список. Стоило упомянуть этот вариант, так как он довольно часто встречается в различных статьях.
01.b — Set Comprehension. Выражение на основе множества
Тоже, что и в предыдущем, но вместо списка здесь множество.
02 — Filter
Функция filter схожа по принципу работы с выражением-генератором. Функция лямбда применяется к каждому элементу последовательности чисел от 1 до 999. Все числа последовательности, делящиеся на 3 или на 5 без остатка от деления, возвращаются, затем суммируются функцией sum() .
03 — Map
Перебираемые числа последовательности от 1 до 999, делящиеся на 3 или 5 без остатка от деления, остаются без изменений, все остальные числа заменяются на ноль. Полученная последовательность суммируется функцией sum() .
04 — Reduce
Из всей подборки, этот вариант «очень не очень». Как по степени реализации, так и по времени выполнения(но об этом попозже).
Если в reduce указан инициализатор(в нашем случае ноль), то он становится накопителем. К нему по очереди прибавляются только те числа из последовательности от 1 до 999, которые делятся на 3 или на 5 без остатка от деления. Если из функции reduce убрать инициализатор ноль, то инициализатором станет крайний левый элемент последовательности.
05 — Однострочное решение на основе множества
Самое элегантное решение, как по красоте написания, так и по времени выполнения.
Последовательность чисел от 1 до 999, кратную трём, помещаем во множество и объединяем со множеством, содержащим в себе последовательность чисел от 1 до 999, кратную пяти. Содержимое, полученного множества суммируем функцией sum() .
05.a — Ещё одно однострочное решение на основе множества
Похожий вариант на предыдущий, но, если использовать фигурные скобки, то последовательность чисел от 1 до 999, кратную трём и последовательность чисел от 1 до 999, кратную пяти, нужно распаковывать.
05.b — И ещё одно однострочное решение на основе множества
Создаём множество, с последовательностью чисел от 1 до 999, кратную трём и присоединяем к нему последовательность чисел от 1 до 999, кратную пяти. Затем функцией sum() суммируем.
05.c И последнее однострочное решение на основе множества
По аналогии с предыдущими. Распаковываем последовательности чисел в списки. Складываем списки. Оборачиваем во множество. Затем суммируем функцией sum() .
Смотрим на скорость выполнения каждого однострочного решения
Если проверить скорость выполнения каждого однострочного решения в командной строке, при помощи timeit, получим следующие значения в микросекундах:
Методика расчёта: python -m timeit «выражение»
Быстрее всего справились с задачей последние четыре варианта.
Заключение
Всего получилось 5 уникальных + 5 не уникальных решений. Благодаря этой задаче у меня появилось более устойчивое понимание работы функций Filter, Map, Reduce. И если раньше я недоумевал, почему функцию Reduce убрали из основного модуля, то теперь я не сомневаюсь в правильности этого решения.
В статье я старался отойти от популярного шаблона повествования «точность на грани бесполезности». Где предложения набиты под завязку «тяжёлыми» терминами, а из знакомого там только союзы и предлоги. Не уверен, что у меня получилось.
http://slemeshevsky.github.io/num-mmf/ode/html/._ode-FlatUI001.html
http://habr.com/ru/post/585176/