Метод эйлера решение систем уравнений

Метод эйлера решение систем уравнений

Системой дифференциальных уравнений называется система вида

где x — независимый аргумент,

y_i — зависимая функция, ,

Функции y_i(x), при подстановке которой система уравнений обращается в тождество, называется решением системой дифференциальных уравнений.

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

Модифицированный метод Эйлера.

Метод Рунге-Кутта четвертого порядка.

Дифференциальным уравнением второго порядка называется уравнение вида

Функция y(x), при подстановке которой уравнение обращается в тождество, называется решением дифференциального уравнения.

Численно ищется частное решение уравнения (2), которое удовлетворяет заданным начальным условиям, то есть решается задача Коши.

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

.

(3)

Функция f₂(x, y₁, y) в систему (3) введена формально для того, чтобы методы, которые будут показаны ниже, могли быть использованы для решения произвольной системы дифференциальных уравнений первого порядка. Рассмотрим несколько численных методов решения системы (3). Расчетные зависимости для i+1 шага интегрирования имеют следующий вид. Для решения системы из n уравнений расчетные формулы приведены выше. Для решения системы из двух уравнений расчетные формулы удобно записать без двойных индексов в следующем виде:

Метод Рунге-Кутта четвертого порядка.

где h — шаг интегрирования. Начальные условия при численном интегрировании учитываются на нулевом шаге: i=0, x=x₀, y₁=y₁₀, y=y₀.

Контрольное задание по зачетной работе.

Колебания с одной степенью свободы

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

Задание. Численно и аналитически найти:

1. закон движения материальной точки на пружинке х(t),
2. закон изменения силы тока I(t) в колебательном контуре (RLC — цепи) для заданных в табл.1,2 режимов. Построить графики искомых функций.

Свободные незатухающие колебания

Затухающее колебательное движение

Предельное апериодическое движение

Вынужденное колебание без сопротивления

Вынужденное колебание без сопротивления, явление резонанса

Вынужденное колебание с линейным сопротивлением

Вынужденное колебание с линейным сопротивлением, явление резонанса

Базовая алгебра и вычисления¶

Sage может осуществлять вычисления такие, как поиск решений уравнений, дифференцирование, интегрирование и преобразования Лапласа. См. Sage Constructions , где содержатся примеры.

Решение уравнений¶

Точное решение уравнений¶

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

Можно решать уравнения для одной переменной через другие:

Также можно решать уравнения с несколькими переменными:

Следующий пример показывает, как Sage решает систему нелинейных уравнений. Для начала система решается символьно:

Для приближенных значений решения можно использовать:

(Функция n выведет приближенное значение. Аргументом для данной функции является количество битов точности)

Численное решение уравнений¶

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

С другой стороны функция find_root может использоваться для решения вышеуказанного примера в интервале \(0 :

Дифференцирование, интегрирование и т.д.¶

Sage умеет дифференцировать и интегрировать многие функции. Например, для того, чтобы продифференцировать \(\sin(u)\) по переменной \(u\) , требуется:

Для подсчета четвертой производной функции \(\sin(x^2)\) надо:

Для нахождения частных производных, как, например, для функции \(x^2+17y^2\) по \(x\) и \(y\) соответственно:

Теперь найдём интегралы: и определенные, и неопределенные. Например, \(\int x\sin(x^2)\, dx\) и \(\int_0^1 \frac\, dx\)

Для нахождения разложения на простые дроби для \(\frac<1>\) нужно сделать следующее:

Решение дифференциальных уравнений¶

Sage может использоваться для решения дифференциальных уравнений. Для решения уравнения \(x’+x-1=0\) сделаем следующее:

Для этого используется интерфейс Maxima [Max] , поэтому результат может быть выведен в виде, отличном от обычного вывода Sage. В данном случае общее решение для данного дифференциального уравнения — \(x(t) = e^<-t>(e^+C)\) .

Преобразования Лапласа также могут быть вычислены. Преобразование Лапласа для \(t^2e^t -\sin(t)\) вычисляется следующим образом:

Приведем более сложный пример. Отклонение от положения равновесия для пары пружин, прикрепленных к стене слева,

может быть представлено в виде дифференциальных уравнений второго порядка

где \(m_\) — это масса объекта i, \(x_\) — это отклонение от положения равновесия массы i, а \(k_\) — это константа для пружины i.

Пример: Используйте Sage для вышеуказанного примера с \(m_<1>=2\) , \(m_<2>=1\) , \(k_<1>=4\) , \(k_<2>=2\) , \(x_<1>(0)=3\) , \(x_<1>‘(0)=0\) , \(x_<2>(0)=3\) , \(x_<2>‘(0)=0\) .

Решение: Надо найти преобразование Лапласа первого уравнения (с условием \(x=x_<1>\) , \(y=x_<2>\) ):

Данный результат тяжело читаем, однако должен быть понят как

Найдем преобразование Лапласа для второго уравнения:

Вставим начальные условия для \(x(0)\) , \(x'(0)\) , \(y(0)\) и \(y'(0)\) , и решим уравения:

Теперь произведём обратное преобразование Лапласа для нахождения ответа:

График для ответа может быть построен параметрически, используя

Графики могут быть построены и для отдельных компонентов:

Для более исчерпывающей информации по графикам см. Построение графиков . Также см. секцию 5.5 из [NagleEtAl2004] для углубленной информации по дифференциальным уравнениям.

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

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

требуется найти приблизительное значение решения при \(x=b\) и \(b>a\) .

Из определения производной следует, что

где \(h>0\) дано и является небольшим. Это и дифференциальное уравнение дают \(f(x,y(x))\approx \frac\) . Теперь надо решить для \(y(x+h)\) :

Если назвать \(h\cdot f(x,y(x))\) «поправочным элементом», \(y(x)\) «прежним значением \(y\) » а \(y(x+h)\) «новым значением \(y\) », тогда данное приближение может быть выражено в виде

Если разбить интервал между \(a\) и \(b\) на \(n\) частей, чтобы \(h=\frac\) , тогда можно записать информацию для данного метода в таблицу.

Целью является заполнить все пустоты в таблице по одному ряду за раз до момента достижения записи . которая и является приближенным значением метода Эйлера для \(y(b)\) .

Решение систем дифференциальных уравнений похоже на решение обычных дифференциальных уравнений.

Пример: Найдите численное приблизительное значение для \(z(t)\) при \(t=1\) , используя 4 шага метода Эйлера, где \(z»+tz’+z=0\) , \(z(0)=1\) , \(z'(0)=0\) .

Требуется привести дифференциальное уравнение 2го порядка к системе двух дифференцальных уравнений первого порядка (используя \(x=z\) , \(y=z’\) ) и применить метод Эйлера:

Итак, \(z(1)\approx 0.75\) .

Можно построить график для точек \((x,y)\) , чтобы получить приблизительный вид кривой. Функция eulers_method_2x2_plot выполнит данную задачу; для этого надо определить функции f и g, аргумент которых имеет три координаты: ( \(t\) , \(x\) , \(y\) ).

В этот момент P содержит в себе два графика: P[0] — график \(x\) по \(t\) и P[1] — график \(y\) по \(t\) . Оба эти графика могут быть выведены следующим образом:

Специальные функции¶

Несколько ортогональных полиномов и специальных функций осуществлены с помощью PARI [GAP] и Maxima [Max] .

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

Задачи с начальными условиями для систем обыкновенных дифференциальных уравнений

Рассмотрим задачу Коши для системы обыкновенных дифференциальных уравнений $$ \begin \tag <1>\frac &= f_i (t, u_1, u_2, \ldots, u_n), \quad t > 0\\ \tag <2>u_i(0) &= u_i^0, \quad i = 1, 2, \ldots, m. \end $$

Используя векторные обозначения, задачу (1), (2) можно записать как задачу Коши $$ \begin \tag <3>\frac> &= \pmb(t, \pmb), \quad t > 0, \\ \tag <4>\pmb(0) &= \pmb_0 \end $$ В задаче Коши необходимо по известному решению в точке \( t = 0 \) необходимо найти из уравнения (3) решение при других \( t \).

Численные методы решения задачи Коши

Существует большое количество методов численного решения задачи (3), (4). Вначале рассмотрим простейший явный метод Эйлера и его программную реализацию. Затем будут представлены методы Рунге—Кутта и многошаговые методы.

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

Идея численных методов решения задачи (3), (4) состоит из четырех частей:

1. Вводится расчетная сетка по переменной \( t \) (время) из \( N_t + 1 \) точки \( t_0 \), \( t_1 \), \( \ldots \), \( t_ \). Нужно найти значения неизвестной функции \( \pmb \) в узлах сетки \( t_n \). Обозначим через \( \pmb^n \) приближенное значение \( \pmb(t_n) \).

2. Предполагаем, что дифференциальное уравнение выполнено в узлах сетки.

3. Аппроксимируем производные конечными разностями.

4. Формулируем алгоритм, который вычисляет новые значения \( \pmb^ \) на основе предыдущих вычисленных значений \( \pmb^k \), \( k 0 \) при \( \tau\to 0 \).

Явный метод Эйлера

Проиллюстрируем указанные шаги. Для начала введем расчетную сетку. Очень часто сетка является равномерной, т.е. имеет одинаковое расстояние между узлами \( t_n \) и \( t_ \): $$ \omega_\tau = \< t_n = n \tau, n = 0, 1, \ldots, N_t \>. $$

Затем, предполагаем, что уравнение выполнено в узлах сетки, т.е.: $$ \pmb^\prime (t_n) = \pmb(t_n, u(t_n)), \quad t_n \in \omega_\tau. $$

Заменяем производные конечными разностями. С этой целью, нам нужно знать конкретные формулы, как производные могут быть аппроксимированы конечными разностями. Простейший подход заключается в использовании определения производной: $$ \pmb^\prime(t) = \lim_ <\tau \to 0>\frac<\pmb(t+\tau) — \pmb(t)><\tau>. $$

В произвольном узле сетки \( t_n \) это определение можно переписать в виде: $$ \begin \pmb^\prime(t_n) = \lim_ <\tau \to 0>\frac<\pmb(t_n+\tau) — \pmb(t_n)><\tau>. \end $$ Вместо того, чтобы устремлять шаг сетки к нулю, мы можем использовать малый шаг \( \tau \), который даст численное приближение \( u^\prime(t_n) \): $$ \begin \pmb^\prime(t_n) \approx \frac<\pmb^ — \pmb^><\tau>. \end $$ Такая аппроксимация известна как разностная производная вперед и имеет первый порядок по \( \tau \), т.е. \( O(\tau) \). Теперь можно использовать аппроксимацию производной. Таким образом получим явный метод Эйлера: $$ \begin \tag <5>\frac<\pmb^ — \pmb^n> <\tau>= \pmb(t_n, \pmb^). \end $$

Четвертый шаг заключается в получении численного алгоритма. Из (5) следует, что мы должны знать значение \( y^n \) для того, чтобы решить уравнение (5) относительно \( y^ \) и получить формулу для нахождения приближенного значения искомой функции на следующем временном слое \( t_ \): $$ \begin \tag <6>\pmb^ = \pmb^n + \tau \pmb(t_n, \pmb^) \end $$

При условии, что у нас известно начальное значение \( \pmb^0 = \pmb_0 \), мы можем использовать (6) для нахождения решений на последующих временных слоях.

Программная реализация явного метода Эйлера

Выражение (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 \tag <7>\frac<\pmb^ — \pmb^n> <\tau>= \pmb(t_, \pmb^). \end $$

Таким образом для нахождения приближенного значения искомой функции на новом временном слое \( t_ \) нужно решить нелинейное уравнение относительно \( \pmb^ \): $$ \begin \tag <8>\pmb^ — \tau \pmb(t_, \pmb^) — y^n = 0. \end $$

Для решения уравнения (8) можно использовать, например, метод Ньютона.

Программная реализация неявного метода Эйлера

Функция backward_euler решения системы уравнений реализована в файле euler.py:

Отметим, что для нахождения значения u[n+1] используется функция fsolve модуля optimize библиотеки scipy . В качестве начального приближения для решения нелинейного уравнения используется значение искомой функции с предыдущего слоя u[n] .

Методы Рунге—Кутта

Одношаговый метод Рунге—Кутта в общем виде записывается следующим образом: $$ \begin \tag <9>\frac<\pmb^ — \pmb^n> <\tau>= \sum_^s b_i \pmb_i, \end $$ где $$ \begin \tag <10>\pmb_i = \pmb\left( t_n + c_i\tau, \pmb^n + \tau \sum_^s a_\pmb_j \right), \quad i = 1, 2, \ldots, s. \end $$ Формула (9) основана на \( s \) вычислениях функции \( \pmb \) и называется \( s \)-стадийной. Если \( a_ = 0 \) при \( j \geq i \) имеем явный метод Рунге—Кутта. Если \( a_ = 0 \) при \( j > i \) и \( a_ \ne 0 \), то \( \pmb_i \) определяется неявно из уравнения $$ \begin \tag <11>\pmb_i = \pmb\left( t_n + c_i\tau, \pmb^n + \tau \sum_^ a_\pmb_j + \tau a_ \pmb_i \right), \quad i = 1, 2, \ldots, s. \end $$ О таком методе Рунге—Кутта говорят как о диагонально-неявном.

Одним из наиболее распространенных является явный метод Рунге-Кутта четвертого порядка: $$ \begin \tag <12>\pmb_1 & = \pmb(t_n, \pmb^n), &\quad \pmb_2 &= \pmb\left( t_n + \frac<\tau><2>, \pmb^n + \tau \frac<\pmb_1> <2>\right),\\ \pmb_3 &= \pmb\left( t_n + \frac<\tau><2>, \pmb^n + \tau \frac<\pmb_2> <2>\right), &\quad \pmb_4 &= \pmb\left( t_n + \tau, \pmb^n + \tau \pmb_3 \right),\\ \frac<\pmb^ -\pmb^n> <\tau>&= \frac<1> <6>(\pmb_1 + 2\pmb_2 + 2\pmb_3 + \pmb_4) & & \end $$

Многошаговые методы

В методах Рунге—Кутта в вычислениях участвуют значения приближенного решения только в двух соседних узлах \( \pmb^n \) и \( \pmb^ \) — один шаг по переменной \( t \). Линейный \( m \)-шаговый разностный метод записывается в виде $$ \begin \tag <13>\frac<1> <\tau>\sum_^m a_i \pmb^ = \sum_^ b_i \pmb(t_, \pmb^), \quad n = m-1, m, \ldots \end $$ Вариант численного метода определяется заданием коэффициентов \( a_i \), \( b_i \), \( i = 0, 1, \ldots, m \), причем \( a_0 \ne 0 \). Для начала расчетов по рекуррентной формуле (13) необходимо задать \( m \) начальных значений \( \pmb^0 \), \( \pmb^1 \), \( \dots \), \( \pmb^ \) (например, можно использовать для их вычисления метод Эйлера).

Различные варианты многошаговых методов (методы Адамса) решения задачи с начальными условиями для систем обыкновенных дифференциальных уравнений могут быть получены на основе использования квадратурных формул для правой части равенства $$ \begin \tag <14>\pmb(t_) — \pmb(t_n) = \int_^> \pmb(t, \pmb) dt \end $$

Для получения неявного многошагового метода используем для подынтегральной функции интерполяционную формулу по значениям функции \( \pmb^ = \pmb(t_, \pmb^) \), \( \pmb^n \), \( \dots \), \( \pmb^ \), т.е. $$ \begin \tag <15>\frac<\pmb^ — \pmb^n> <\tau>= \sum_^ b_i \pmb(t_, \pmb^) \end $$

Для интерполяционного метода Адамса (15) наивысший порядок аппроксимации равен \( m+1 \).

Для построения явных многошаговых методов можно использовать процедуру экстраполяции подынтегральной функции в правой части (14). В этом случае приближение осуществляется по значениям \( \pmb^n \), \( \pmb^ \), \( \dots \), \( \pmb^ \) и поэтому $$ \begin \tag <16>\frac<\pmb^ — \pmb^n> <\tau>= \sum_^ b_i \pmb(t_, \pmb^) \end $$

Для экстраполяционного метода Адамса (16) погрешность аппроксимации имеет \( m \)-ый порядок.

На основе методов Адамса строятся и схемы предиктор–корректор. На этапе предиктор используется явный метод Адамса, на этапе корректора — аналог неявного метода Адамса. Например, при использовании методов третьего порядка аппроксимации в соответствии с (18) для предсказания решения положим $$ \frac<\pmb^ — \pmb^n> <\tau>= \frac<1> <12>(23 \pmb^ -16\pmb^ + 5\pmb^). $$ Для уточнеия решения (см. (17)) используется схема $$ \frac<\pmb^ — \pmb^n> <\tau>= \frac<1> <24>(9\pmb^ + 19\pmb^ — 5\pmb^ + \pmb^). $$ Аналогично строятся и другие классы многошаговых методов.

Жесткие системы ОДУ

При численном решении задачи Коши для систем обыкновенных дифференциальных уравнений (3), (4) могут возникнуть дополнительные трудности, порожденные жесткостью системы. Локальные особенности поведения решения в точке \( u = w \) передаются линейной системой $$ \begin \frac