No Image

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

СОДЕРЖАНИЕ
739 просмотров
10 марта 2020

Главная Цены Оплата Примеры решений Отзывы Ccылки Теория Книги Сотрудничество Форум
Теория / Информатика / Лекция 14. Численное решение системы дифференциальных уравнений

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

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

yi – зависимая функция, ,

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

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

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

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

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

F(x,y,у’,y")=0 (1)

y"=f(x,y,y’). (2)

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

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

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

. (3)

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

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

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

Используя векторные обозначения, задачу (1), (2) можно записать как задачу Коши $$ egin ag <3>frac> &= pmb(t, pmb), quad t > 0, \ ag <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 ) при ( au o 0 ).

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

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

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

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

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

Читайте также:  Почему вайфай не работает на айфоне

Четвертый шаг заключается в получении численного алгоритма. Из (5) следует, что мы должны знать значение ( y^n ) для того, чтобы решить уравнение (5) относительно ( y^ ) и получить формулу для нахождения приближенного значения искомой функции на следующем временном слое ( t_ ): $$ egin ag <6>pmb^ = pmb^n + au 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) используется следующий метод: $$ egin ag <7>frac<pmb^ – pmb^n> < au>= pmb(t_, pmb^). end $$

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

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

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

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

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

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

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

Читайте также:  Распечатать разлиновку для письма

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

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

В методах Рунге—Кутта в вычислениях участвуют значения приближенного решения только в двух соседних узлах ( pmb^n ) и ( pmb^ ) — один шаг по переменной ( t ). Линейный ( m )-шаговый разностный метод записывается в виде $$ egin ag <13>frac<1> < au>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
e 0 ). Для начала расчетов по рекуррентной формуле (13) необходимо задать ( m ) начальных значений ( pmb
^0 ), ( pmb^1 ), ( dots ), ( pmb^ ) (например, можно использовать для их вычисления метод Эйлера).

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

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

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

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

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

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

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

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

= sum_^ frac<partial f_i> <partial u_j>(t, w) v + ar(t), quad t > 0. end $$

Пусть ( lambda_i(t) ), ( i = 1, 2, ldots, m ) — собственные числа матрицы $$ egin A(t) = < a_(t) >, quad a_(t) = frac<partial f_i><partial u_j>(t, w). end $$ Система уравнений (3) является жесткой, если число $$ egin S(t) = frac <max_<1 leq i leq m>|Re lambda_i(t)|> <min_<1 leq i leq m>|Re lambda_i(t)|> end $$ велико. Это означает, что в решении присутствуют составляющие с сильно различающимися масштабами изменения по переменной ( t ).

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

Читайте также:  Как обработать видео в vsco

Метод называется ( A )-устойчивым, если при решении задачи Коши для системы (3) область его устойчивости содержит угол $$ egin |arg(-mu)| –>

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

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

, где x – независимая переменная, – i-ая производная от искомой функции. n – порядок уравнения. Общее решение ОДУ n–го порядка содержит n произвольных постоянных , т.е. общее решение имеет вид .

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

Ясно, что при n=1 можно говорить только о задачи Коши.

Примеры постановки задачи Коши:

Примеры краевых задач:

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

Численные методы решения задачи Коши для ОДУ первого порядка

Постановка задачи. Найти решение ОДУ первого порядка

на отрезке при условии

При нахождении приближенного решения будем считать, что вычисления проводятся с расчетным шагом , расчетными узлами служат точки промежутка [x, xn].

Целью является построение таблицы

т.е. ищутся приближенные значения y в узлах сетки.

Интегрируя уравнение на отрезке , получим

Вполне естественным (но не единственным) путем получения численного решения является замена в нем интеграла какой–либо квадратурной формулой численного интегрирования. Если воспользоваться простейшей формулой левых прямоугольников первого порядка

,

то получим явную формулу Эйлера:

, .

Зная , находим , затем т.д.

Геометрическая интерпретация метода Эйлера:

Пользуясь тем, что в точке x известно решение y(x) = y и значение его производной , можно записать уравнение касательной к графику искомой функции в точке :. При достаточно малом шаге h ордината этой касательной, полученная подстановкой в правую часть значения , должна мало отличаться от ординаты y(x1) решенияy(x) задачи Коши. Следовательно, точка пересечения касательной с прямой x = x1 может быть приближенно принята за новую начальную точку. Через эту точку снова проведем прямую , которая приближенно отражает поведение касательной к в точке . Подставляя сюда (т.е. пересечение с прямой x = x2), получим приближенное значение y(x) в точке x2: и т.д. В итоге для i–й точки получим формулу Эйлера.

Явный метод Эйлера имеет первый порядок точности или аппроксимации.

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

, .

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

Неявный метод Эйлера имеет первый порядок точности или аппроксимации.

Модифицированный метод Эйлера: в данном методе вычисление состоит из двух этапов:

Данная схема называется еще методом предиктор – корректор (предсказывающее – исправляющее). На первом этапе приближенное значение предсказывается с невысокой точностью (h), а на втором этапе это предсказание исправляется, так что результирующее значение имеет второй порядок точности.

Методы Рунге – Кутта: идея построения явных методов Рунге–Кутты p–го порядка заключается в получении приближений к значениям y(xi+1) по формуле вида

,

.

Здесь an, bnj, pn, – некоторые фиксированные числа (параметры).

При построения методов Рунге–Кутты параметры функции (an, bnj, pn) подбирают таким образом, чтобы получить нужный порядок аппроксимации.

Схема Рунге – Кутта четвертого порядка точности:

Пример. Решить задачу Коши:

.

Рассмотреть три метода: явный метод Эйлера, модифицированный метод Эйлера, метод Рунге – Кутта.

Точное решение:

Расчетные формулы по явному методу Эйлера для данного примера:

Расчетные формулы модифицированного метода Эйлера:

Комментировать
739 просмотров
Комментариев нет, будьте первым кто его оставит

Это интересно
Adblock
detector