Метод Эйлера с уточнением
Дата добавления: 2013-12-23 ; просмотров: 4250 ; Нарушение авторских прав
Задача Коши
Дифференциальных уравнений
Численные методы решения обыкновенных
Раздел № 7
связывающее неизвестную функцию у(х),независимую переменную х и производные у’(x), у»(x). у ( n ) (x) неизвестной функции, называется обыкновенным дифференциальным уравнением. Порядок п старшей производной называется порядком дифференциального уравнения.
В задаче Коши для дифференциального уравнения п-го порядкаискомая функция у(х),кроме самого дифференциального уравнения, должна удовлетворять начальным условиям
Методы решения задач для дифференциальных уравнений можно разбить на три типа: точные, приближенные и численные.
Точныминазывают методы, с помощью которых решение дифференциального уравнения можно выразить через известные функции (элементарные функции или интегралы от элементарных функций). Точные методы решения известны только для некоторых классов дифференциальных уравнений (линейные дифференциальные уравнения, уравнения с разделяющимися переменными и др.).
Приближенными называют методы, в которых решение находят как предел последовательности функций, являющихся элементарными или интегралами от элементарных функций. Например, метод разложения искомой функции в ряд Тейлора является приближенным методом.
Численный метод решения дифференциального уравнения – алгоритм вычисления значений искомого решения у(х) на некотором дискретном множестве значений аргумента х. При этом вычисляемые значения искомого решения у(х) являются приближенными, но могут быть и точными.
Рассмотрим задачу Коши для обыкновенного дифференциального уравнения первого порядка
Требуется найти функцию у = у(х), которая удовлетворяет уравнению (7.1) на интервале (х0, X) и начальному условию (7.2) в точке х0.
Приведем без доказательства теорему существования и единственности решения задачи Коши.
Теорема 7.1. Пусть в области R<(x, у), |х — х0| £ а, |у — у0| £ b> функция f(x, у) непрерывна. Тогда на некотором отрезке |х — x0| £ d существует решение уравнения (7.1), удовлетворяющее условию (7.2).
Если в области R функция f(x, у) удовлетворяет условию Липшица
то указанное решение единственно.
Произведем разбиение отрезка [х0, X] на п частей:
xi = x0 + ih, . (7.3)
Найдем приближенные значения решения у(х) в точках xi , i = 1, 2. n.
7.2. Метод Эйлера (ломаных)
Рассмотрим уравнение (7.1) в точках xi, i = 0, 1, . , n — 1 и заменим производную y‘(xi) разностной формулой
y‘(xi) = (7.4)
Тогда получим рекуррентную формулу метода Эйлера для вычисления приближенных значений у(xi + 1):
На рис. 7.1 дана геометрическая иллюстрация метода Эйлера (7.5). Уравнение касательной к графику решения у(х) в точке (х0, у0) имеет вид
удовлетворяющей уравнению (7.1) и начальному условию y(x1) = y1.
Таким образом, с каждым шагом i метод Эйлера (7.5) дает точки (xi, yi), которые, вообще говоря, удаляются от интегральной кривой, соответствующей точному решению задачи Коши (7.1), (7.2). Вместо интегральной кривой метод Эйлера дает ломаную, изображенную на рис. 7.1, поэтому его часто называют методом ломаных.
Формулу (7.5) можно получить и другим способом. Рассмотрим разложение искомого решения у(х) в ряд Тейлора в окрестности точки х0:
у(х) = у(х0) + у ¢ (х0)(х — х0) + + . .
Также здесь получен следующий результат – погрешность вычисления значения у1 есть величина порядка O(h 2 ), а погрешность значения уп – величина порядка O(h). Из-за большой погрешности метод Эйлера применяется редко.
Для повышения точности метода Эйлера применяют следующий прием. Сначала находят приближенное значение решения по методу Эйлера:
,
а затем уточняют его по формуле
.
Этот метод называется методом Эйлера-Коши, и рекуррентные соотношения для его реализации могут быть записаны в виде
, (7.7)
Метод Эйлера-Коши имеет погрешность порядка O(h 2 ). Геометрическая иллюстрация метода Эйлера-Коши показана на рис. 7.2. Очередное приближение метода Эйлера-Коши соответствует точке пересечения диагоналей параллелограмма, построенного на двух звеньях ломаной метода Эйлера.
Метод Эйлера-Коши является одним из частных случаев методов Рунге-Кутта.
Задачи с начальными условиями для систем обыкновенных дифференциальных уравнений
Рассмотрим задачу Коши для системы обыкновенных дифференциальных уравнений $$ \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
Модифицированный метод Эйлера с итерационным уточнением и переменным шагом Текст научной статьи по специальности « Математика»
Аннотация научной статьи по математике, автор научной работы — Трубников С. В.
Предложен новый численный метод решения задач Коши для обыкновенных дифференциальных уравнений .
Похожие темы научных работ по математике , автор научной работы — Трубников С. В.
Текст научной работы на тему «Модифицированный метод Эйлера с итерационным уточнением и переменным шагом»
МОДИФИЦИРОВАННЫЙ МЕТОД ЭЙЛЕРА С ИТЕРАЦИОННЫМ УТОЧНЕНИЕМ
И ПЕРЕМЕННЫМ ШАГОМ
Предложен новый численный метод решения задач Коши для обыкновенных дифференциальных уравнений.
Ключевые слова: численный метод, задача Коши, обыкновенное дифференциальное уравнение, переменный шаг.
Значительную часть прикладных математических задач составляют краевые задачи для дифференциальных уравнений. Поэтому одной из центральных проблем современной прикладной математики является разработка и исследование численных методов решения подобных краевых задач. Этой проблеме посвящена обширная литература (см., например, библиографию в [1], [2] и [3]). В статье [5] описан новый подход к построению численных методов решения задач Коши для обыкновенных дифференциальных уравнений, основанный на аппроксимации приближенного решения с помощью кусочно-полиномиальной интерполяции многочленами Эрмита, а также на принципе минимизации невязки. С помощью этого подхода был построен новый численный метод для решения задачи Коши, названный исправленным методом Эйлера с итерационным уточнением, который можно отнести к итерационно-разностным методам. В нем используется сетка с постоянным шагом. В данной статье описан новый численный метод с переменным шагом сетки и описана процедура автоматического выбора узлов сетки, обеспечивающих заданную точность. Этот метод также относится к итерационно-разностным методам.
1. Составная функциональная кинематическая кривая и ее свойства
Кинематические кривые введены и описаны в статье [4]. Функции, задающие составные кинематические кривые, представляют собой результат кусочно-многочленной интерполяции многочленами Эрмита 5 порядка с двумя трехкратными узлами. Уравнение составной кинематической кривой на плоскости хОу можно записать в виде:
г = г(г) = ^ Т5 (г — г +1)-Оу,_1 + ^ Т5 (г — г +1)-Q5-ji при г е [/-1,/], г = 1,2,* , т . (1)
Здесь т — заданное натуральное число, О j г = Qxj г -1 + Qyji’ ‘] ( j = 0,1,2, г = 0,1, * , т ) —
заданные векторы, Т5 (г) — интерполяционные многочлены Эрмита 5 порядка [2], которые можно записать в виде:
Т05 () = 1 -10*3 +15*4 — бг5, Т5(г)=г -бг3 + 8*4 -3*5, Т25(г) = 1 г2 -2г3 + 2г4 -2г5,
Т35 (г ) = 2 г3 — г4 + 2 г5, Т45 (г ) = -4г3 + 7г4 — 3г5, Т55 (г ) = 10г3 — 15г4 + бг5. (2)
Если ввести обозначения компонент векторной функции
г(г) = Гх (г)- 1 + Гу (г)-] , (3)
то равенство (1) можно записать в виде:
ГХ(і) = XТ! (і -і +1)‘ °х/г-1 + XГ/(і -і + !)• °х5-іі при іє[і -1,і], і = ІД* ,т • (4)
гу (і) = XТ/ (і-*’ +^ 6_УІ і-1 + XТ/ (і-/ +^ Оу5-і і при і є[/-1,/], і = 1,2,‘ , т • (5)
Одно из главных свойств кинематической кривой [4] состоит в том, что
Гх (і) = 6×0 і , _ 6×1 і, ё [Х2 ) _ 6х2 і, і = 0,1,‘ ,т ,
Введем новую кривую, задаваемую уравнением (3) и уравнениями
Гх (і ) = Хі-1 +(хі — Хі-1 У(і — і + 1) при і є[і -1> і] , і = 1,2
гу (г) = X Т](г — г + 1)-^М + ХТ7(г — г + 1)-Qy/ 5-ji при г е[ — ^ г], г = 1,2,к , т . (9)
х0 Не можете найти то, что вам нужно? Попробуйте сервис подбора литературы.
ё2гх (г — 0) . ё2гх (г + 0) .
Из условий (11), (12), (15) — (21) следует, что
У(х ) = б ёУ(хг — 0) = ^уИ г ё2У(хг — 0) = @у12 г
ёх х — х-1 ёх2 (х,- — хг-1 )2
поскольку хг е [хг -1 , хг ]. С другой стороны
у(х ) = б ёу(хг + 0) = буг1 г ёу(хг + 0) = буг 2 г
^ х+1 — х ёх2 (х,-+1 — хг )2
буг 0,=бу,0 г. , ( буг 2 г )2 = ( бу12 г )2 , г = 1,2,* , т -1. (22)
х+1 — хг хг — хг-1 (хг+1 — хг )2 (хг — хг-1 )2
Эти условия гарантируют однозначность и непрерывность функции у(х) = Гу (гх ( 1)( х)) и её производных до второго порядка включительно на [х0 , хт ] и выполнение следующих равенств.
/ ) = б ёу(Х0) = буг 1 0 а2у<х0) = Оуг2 0
у(х ) = О йУ(хт ) = буЛ т й2У(хт ) = бу/2
у(хі ) = буг 0 і = бу/0 і , й2 у (хі ) _ буг 2 і = Не можете найти то, что вам нужно? Попробуйте сервис подбора литературы.
^ ё2гту (г) . Фтх () — ёгту (г) _ ё2гтх ()
Для определения приближенного решения остается вычислить значения величин
Qугj 0 , Qу/jm ( j = 0,1,2) Qугji, Qу/ji ( г = 1,2,* , т -1, j = 0,1,2), которые удовлетворяют
условиям (22) и определяют функцию гту (г) и приближенное решение ут (х) = гту (гтх( 1)(х)).
4. Получение приближенного решения задачи Коши
Для определения неизвестных постоянных мы, так же как в работе [5], введем и будем использовать понятие невязки. Невязкой дифференциального уравнения (24) на его прибли-
женном решении ут (х) = гту (г,
ту тх мы назовем величину
ёУт (х) — I(^ ут (х)), х е [a, Ь] •
ёКт (х ) = ё 2 ут (х ) Э/ (Х ут (х )) Э/ (Х ут (х )) Фт (х )
Учитывая (30), отсюда получим значение невязки и её производной в узлах хг:
(х0 ) = ^у1° — 1 ^ ^г00 ^ Кт (хт ) = —-У (хт, ^0т )
-/(х,Оуг0, ) = -/(х,Qуl0г), г = 1,2,- , т-1, (36)
хг+1 — хг хг — хг-1
ёКт (х0 ) = Qуг 2 0 — Э/(Хo_’Qут0o) Э/ (х0, ^ут 00 ) Qуг 1 0 ёх (х1 — х0 )2 Эх Эу х1 — х0 ’
ёКт (хт ) = Qуl 2 т Э/ (хт, ^у10 т ) Э/ (х0, Qу/ 0 т ) Qу/1 т
ёх (хт — хт-1 )2 Эх Эу хт — хт-1
ёЯт (х, ) = Qуг 2 г -Э/ІХi,Qуr0i]-Э/ІХi,Q^0i] ^г1 г =
ёх (х,+1 — х, )2 Эх Эу х,+1- х,
= Qуl 2 г -Э/(х^бу/0г)-Э/(iQd0i) Qуl 1 г
(х, — х,-1 )2 Эх эу хг — хг-1
Для определения неизвестных постоянных мы будем использовать принцип предельного обнуления невязки, то есть будем требовать выполнение условий, которые бы способствовали тому, чтобы невязка стремилась к нулю, когда т ® ¥, а расстояния между узлами х, стремятся к нулю.
Потребуем, чтобы для приближения ут (х) = гту (гтх^ ^(х)) в точках х, невязка и её
производная обращались в ноль. Исходя из этих требований и (37) получаются выражения
для Qуг 1 г, Qуl 1 г, Qуг 2 г и Qуl 2 г через Qуг 0 г и Qуl 0 г:
^ 10 =(х1 — х0)-у (x0, ^ 00 ),
^ 1 г = (хг -хг-1Уу(хг,^0г^ ^ 1 г = (хг+1 -хгУу(хг,^0г^ г = 1,2,* ,т-1,
Qуl 1 т = (хт хт-1 )-1 (хт, Qуl 0 т ) . (38)
( )2 Э/(х0,^°уг00 ) ( )Э/(х0, ^°уг00 )
^ ( )2 Э/ (х,, ^г 0,),( )Э/ (х,, ^г 0, )Q , 12 1
^ ( )2 Э/(хг,Qуl0,) . ( )Э/(хг,Qуl0г) ^ 1 2 1
^ 2 г = (хг — х,-1) —-э^— ^хг — х,-1)—-Эу——^ г, г = ‘* , т — ^
Q =( )2 Э/ (хт, Qуl 0 т )+( )Э/ (хт, Qуl 0 т ) Q (39)
Qуl 2 т = чхт — хт-1) + ^т — хт-1) Эу ^^у1 1 т . (39)
Зная Qуl 0, (, = 1,2,* ,т ) и Qуг 0, (, = 0,1,* ,т -1), по формулам (38) и (39) можно найти
^ 1, (г = 1,2,* , т X ^г 1 г (г = 0,1,* ,т -1), а затем ^ 2 г (г = 1,2,* , т X Qуг 2 г
(, = 0,1,* ,т -1). Таким образом, построение приближенного решения свелось к определению ^ 0 г (г = 1,2,* , т ) и ^г 0 г (г = °Л* ,т -1).
Неизвестное значение QуГ 0 0 мы найдем из начального условия (25) и условия (30):
^ 0 0 = ут (х0 ) = ут (а) = у . (40)
Остальные значения Пуг 0 і и Пу/ 0 і будем искать последовательно минимизируя квадраты
величин невязки в серединах отрезков [хі-1 , хі ]. Координаты середин хі _ (хі-1 + хі)/2. В результате возникает цепочка задач минимизации
Эти задачи минимизации заменяют требования Ят (х,) = 0 . Значение (Ят (х, ))2 зависит от Qуг 0,-1 и Qуl 0,. При г = 1 первая из этих величин известна: Qуг 0,-1 = Qуг 00 = у . Таким
образом, (Ят (х ))2 зависит только от одной неизвестной Qуl 01. Решив первую задачу минимизации (41) для (Ят (х))2 как функции одной переменной, мы найдем Qуl 01. Используя связи (22), мы найдем Qуг 01 = Qуl 01. Тогда (Ят (щ ))2 также будет зависеть только от одной неизвестной величины Qуl 02. Решив вторую задачу минимизации (41) для (Ят(щ))2 как функции одной переменной, мы найдем Qуl 0 2. Используя связи (22), мы найдем Qуг 02 = Qуl 02. И так далее. Продолжая этот процесс мы последовательно найдем все неизвестные величины.
На г -ом шаге описанного процесса решения задач минимизации вычисляется точка минимума (Ят (х, ))2 как функции одной переменной т = Qуl 0,. Представление для Ят (х,), как функции т = Qуl 0,, мы получим, подставив в выражение для невязки (34) представление
(28) и вычислив значения Т^ (12) и
(хі хі-1 ) ^т (хі ) 22 Оуг 0 і-1 32 буг 1 і-1 32 Пуг 2 і-1 +
http://slemeshevsky.github.io/num-mmf/ode/html/._ode-FlatUI001.html
http://cyberleninka.ru/article/n/modifitsirovannyy-metod-eylera-s-iteratsionnym-utochneniem-i-peremennym-shagom