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

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

с начальным условием

x = x 0 , y (x 0 )= y 0 (3.2)

Требуется найти решение уравнения на отрезке [а , b ].

Разобьем отрезок [a , b ] на n равных частей и получим последовательность х 0 , х 1 , х 2 ,…, х n , где x i = x 0 + ih (i =0,1,…, n ), а h =(b - a )/ n − шаг интегрирования.

В методе Эйлера приближенные значения у(х i +1 )  y i +1 вычисляются последовательно по формулам:

y i+1 = у i +hf(x i , y i ) (i=0,1,2…) (3.3)

При этом искомая интегральная кривая у=у(х) , проходящая через точку М 0 (х 0 , у 0 ), заменяется ломаной М 0 М 1 М 2 … с вершинами М i (x i , y i ) (i =0,1,2,…); каждое звено М i M i +1 этой ломаной, называемой ломаной Эйлера, имеет направление, совпадающее с направлением той интегральной кривой уравнения (1), которая проходит через точку М i (см. рисунок 2):

Рисунок 2. Вид ломаной Эйлера

Модифицированный метод Эйлера более точен.Сначала вычисляются вспомогательные значения искомой функции у к+1/2 в точках х к+1/2 , затем находится значение правой части уравнения (3.1) в средней точке y k+1/2 =f( xk+1/2 , y k+1/2 ) и определяют у к+ :

Тогда:
(3.4)

Формулы (3.4) − рекуррентные формулы метода Эйлера.

Для оценки погрешности в точке х к проводят вычисления у к с шагом h , затем с шагом 2 h и берут 1/3 разницы этих значений:

,

где у(х) - точное решение дифференциального уравнения.

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

3.2. Метод Рунге-Кутта

Методы Рунге-Кутта обладают следующими свойствами:

Эти методы являются одноступенчатыми: чтобы найти у к+1  нужна информация о предыдущей точке (x к  y к ) 

Методы согласуются с рядом Тейлора вплоть до членов порядка h p  где степень р различна для различных методов и называется порядковым номером или порядком метода 

Они не требуют вычисления производных от f(x  y)  а требуют вычисления самой функции

Алгоритм Рунге-Кутта третьего порядка:

(3.5)

Алгоритм Рунге-Кутта четвертого порядка:

(3.6)

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

3.3. Метод Адамса

Метод Адамса относится к многошаговым схемам решения ДУ, характеризующихся тем, что решение в текущем узле зависит от данных не в одном предыдущем или последующем узле сетки, как это имеет место в одношаговых методах, а зависит от данных в нескольких соседних узлах .

Идея методов Адамса заключается в том, чтобы для повышения точности использовать вычисленные уже на предыдущих шагах значения

Y k -1 , Y k -2 , Y k -3 …

Если используются значения в k предыдущих узлах, то говорят о k-шаговом методе интегрирования уравнения. Одним из способов построения многошаговых методов заключается в следующем. По значениям функции, вычисленным в k предшествующих узлах, строится интерполяционный полином степени (k-1) – L k -1 (x ) , который используется при интегрировании дифференциального уравнения по выражению:

Интеграл при этом выражается через квадратурную формулу:

где λ l – квадратурные коэффициенты.

Полученное таким образом семейство формул называется явной k -шаговой схемой Адамса . Как видно, при k =1 в качестве частного случая получается формула Эйлера.

Например, для формулы 4 порядка имеем:

(3.7)

y ( p ) k +1 – “прогноз” ,вычисленный с использованием значений в предыдущих точках, f ( p ) k +1 –приближенное значение функции,вычисленное в точке получения прогноза, y ( c ) k +1 – «коррекция» прогнозного значения, y k +1 – искомое значение по Адамсу.

Достоинство такого метода решения ДУ заключается в том, что в каждой точке рассчитывается только одно значение функции F(x,y). К недостаткам можно отнести невозможность старта многошагового метода из единственной начальной точки, так как для вычислений по k -шаговой формуле необходимо значение значения функции в k узлах. Поэтому приходится (k-1) решение в первых узлах x 1 , x 2 , …, x k-1 получать с помощью какого-либо одношагового метода, например метода Рунге-Кутта 4–го порядка.

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

4. Краткое описание программы на C++ и представление результатов ее выполнения

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

ОДУ можно записать в виде:
F x , y , y , y ,..., y
где x – независимая переменная,				y i - i -ая производная от

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

c 1 , c 2 ,..., c n , т.е. общее решение имеет вид y x , c 1 , c 2 ,..., c n . Для выделения единственного решения необходимо задать n дополнительных условий. В зависимости от способа задания

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

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

dy x 2 y 3					y 1 1;
d 2 y dy						y 1 1,
dx 2 dx xy ,							y 1 0 .

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

	d 2 y			y sin x ,					y 0 1,	y 1 0
	dx 2
	d 3 y				d 2 y				y 1 0,		y 3 2 .
		x x dx 2				dx ,
										y 1 0,
Решить такие									аналитически удается лишь для

некоторых специальных типов уравнений, поэтому применение приближенных методов решения является необходимостью.

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

Требуется найти решение y (x ) ОДУ первого порядка

f x, y

на отрезке x 0 , x n при условии

y x0 y0 .

Приближенное решение будем искать в узлах расчетной

xi x0 ih,

i 0,1,..., n с

xn x0

Необходимо найти

приближенные

значения в

узлах сетки

y i =y (x i ). Результаты расчетов занесем в таблицу

Интегрируя

уравнение на

отрезке x i , x i

1 , получим

x i 1

y i 1

yi f x, y dx .

Для того, чтобы найти все значения y i , нужно каким-то образом

вычислить интеграл, стоящий в правой части (5.4). Применяя различные квадратурные формулы, будем получать методы решения задачи (5.2), (5.3) разного порядка точности.

Метод Эйлера

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

Явный метод Эйлера имеет первый порядок аппроксимации. Реализация метода. Поскольку x 0 , y 0 , f x 0 , y 0

известны, применяя (5.5) последовательно, определим все y i : y 1 y 0 hf x 0 , y 0 , y 2 y 1 hf x 1 , y 1 , ….

Геометрическая	интерпретация
(рис. 5.1.):
Пользуясь тем, что в точке x 0 известно решение y x 0 y 0
и значение его производной y x 0 dy				f x0 , y0 ,
			x x0

записать уравнение касательной к графику искомой функции
			f x0 , y0		y y0	f x0 , y0 x x0 .
достаточно			шаге h	ордината		y1 y0 hf x0 , y0

касательной, полученная подстановкой в правую часть значения x 1 x 0 h , должна мало отличаться от ординаты y x 1 решения

y x задачи Коши. Следовательно, точка x 1 , y 1 пересечения касательной с прямой x x 1 может быть приближенно принята

за новую начальную точку. Через эту точку снова проведем
прямую y y 1 f x 1 , y 1 x x 1 ,	которая приближенно отражает
поведение касательной к y x	использовать x i 1 прямоугольников: f x, y dx hf xi 1 , yi 1 , то неявный метод Эйлера y i 1 y i hf x i 1 , y i 1 , i 0,1,..., n 1. Этот метод называют неявным, поскольку для вычисления неизвестного значения yi 1 y xi 1 по известному значению требуется решать уравнение, в общем случае нелинейное. Неявный метод Эйлера также имеет первый порядок аппроксимации. Модифицированный метод Эйлера В данном методе вычисление y i 1 состоит из двух этапов: ~ y i 1 y i hf x i , y i , y i 1 f xi , yi f xi 1 , yi 1 Данная схема называется также методом предикторкорректор. Это английское название, означающее «предсказатьисправить». Действительно, на первом этапе приближенное значение предсказывается с первым порядком точности, а на

Известно, что обыкновенное дифференциальное уравнение первого порядка имеет вид: .Решением этого уравнения является дифференцируемая функция, которая при подстановке в уравнение обращает его в тождество. График решения дифференциального уравнения (рис 1.) называетсяинтегральной кривой.

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

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

Задача Коши заключается в отыскании функции , удовлетворяющей исходному уравнению и начальному условию. Обычно определяют решение задачи Коши на отрезке, расположенном справа от начального значения, т. е. для.

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

Численные методы не позволяют найти решение в общем виде, зато они применимы к широкому классу дифференциальных уравнений.

Сходимость численных методов решения задачи Коши. Пусть – решение задачи Коши. Назовем погрешностью численного метода функцию , заданную в узлах сетки. В качестве абсолютной погрешности примем величину.

Численный метод решения задачи Коши называется сходящимся , если для него при. Говорят, что метод имеет-ый порядок точности, если для погрешности справедлива оценка, – константа, .

Метод Эйлера

Простейшим методом решения задачи Коши является метод Эйлера. Будем решать задачу Коши

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

Эти формулы и начальное условие являются расчетными формулами метода Эйлера.

Геометрическая интерпретация одного шага метода Эйлера заключается в том, что решение на отрезке заменяется касательной, проведенной в точкек интегральной кривой, проходящей через эту точку. После выполненияшагов неизвестная интегральная кривая заменяется ломаной линией(ломаной Эйлера).

Оценка погрешности. Для оценки погрешности метода Эйлера воспользуемся следующей теоремой.

Теорема. Пусть функция удовлетворяет условиям:

.

Тогда для метода Эйлера справедлива следующая оценка погрешности: , где– длина отрезка. Мы видим, что метод Эйлера имеет первый порядок точности.

Оценка погрешности метода Эйлера часто бывает затруднительна, так как требует вычисления производных функции . Грубую оценку погрешности даетправило Рунге (правило двойного пересчета), которое используется для различных одношаговых методов, имеющих -ый порядок точности. Правило Рунге заключается в следующем. Пусть– приближения, полученные с шагом, а– приближения, полученные с шагом. Тогда справедливо приближенное равенство:

.

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

Используя правило Рунге, можно построить процедуру приближенного вычисления решения задачи Коши с заданной точностью . Для этого нужно, начав вычисления с некоторого значения шага , последовательно уменьшать это значение в два раза, каждый раз вычисляя приближенное значение,. Вычисления прекращаются тогда, когда будет выполнено условие: . Для метода Эйлера это условие примет вид:. Приближенным решением будут значения,.

Пример 1. Найдем решение на отрезке следующей задачи Коши:,. Возьмем шаг. Тогда.

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

, .

Решение представим в виде таблицы 1:

Таблица 1

Исходное уравнение есть уравнение Бернулли. Его решение можно найти в явном виде: .

Для сравнения точного и приближенного решений представим точное решение в виде таблицы 2:

Таблица 2

Из таблицы видно, что погрешность составляет

Определение дифференциального уравнения Эйлера. Рассмотрены методы его решения.

Содержание

Дифференциальное уравнение Эйлера - это уравнение вида
a 0 x n y (n) + a 1 x n-1 y (n-1) + ... + a n-1 xy′ + a n y = f(x) .

В более общем виде уравнение Эйлера имеет вид:
.
Это уравнение подстановкой t = ax+b приводится к более простому виду, которое мы и будем рассматривать.

Приведение дифференциального уравнения Эйлера к уравнению с постоянными коэффициентами.

Рассмотрим уравнение Эйлера:
(1) .
Оно сводится к линейному уравнению с постоянными коэффициентами подстановкой:
x = e t .
Действительно, тогда
;
;
;

;
;
..........................

Таким образом, множители, содержащие x m , сокращаются. Остаются члены с постоянными коэффициентами. Однако на практике, для решения уравнений Эйлера, можно применять методы решения линейных ДУ с постоянными коэффициентами без использования указанной выше подстановки.

Решение однородного уравнения Эйлера

Рассмотрим однородное уравнение Эйлера:
(2) .
Ищем решение уравнения (2) в виде
.
;
;
........................
.
Подставляем в (2) и сокращаем на x k . Получаем характеристическое уравнение:
.
Решаем его и получаем n корней, которые могут быть комплексными.

Рассмотрим действительные корни. Пусть k i - кратный корень кратности m . Этим m корням соответствуют m линейно независимых решений:
.

Рассмотрим комплексные корни. Они появляются парами вместе с комплексно сопряженными. Пусть k i - кратный корень кратности m . Выразим комплексный корень k i через действительную и мнимую части:
.
Этим m корням и m комплексно сопряженным корням соответствуют 2 m линейно независимых решений:
;
;
..............................
.

После того как получены n линейно независимых решений, получаем общее решение уравнения (2):
(3) .

Примеры

Решить уравнения:


Решение примеров > > >

Решение неоднородного уравнения Эйлера

Рассмотрим неоднородное уравнение Эйлера:
.
Метод вариации постоянных (метод Лагранжа) также применим и к уравнениям Эйлера.

Сначала мы решаем однородное уравнение (2) и получаем его общее решение (3). Затем считаем постоянные функциями от переменной x . Дифференцируем (3) n - 1 раз. Получаем выражения для n - 1 производных y по x . При каждом дифференцировании члены, содержащие производные приравниваем к нулю. Так получаем n - 1 уравнений, связывающих производные . Далее находим n -ю производную y . Подставляем полученные производные в (1) и получаем n -е уравнение, связывающее производные . Из этих уравнений определяем . После чего интегрируя, получаем общее решение уравнения (1).

Пример

Решить уравнение:

Решение > > >

Неоднородное уравнение Эйлера со специальной неоднородной частью

Если неоднородная часть имеет определенный вид, то получить общее решение проще, найдя частное решение неоднородного уравнения. К такому классу относятся уравнения вида:
(4)
,
где - многочлены от степеней и , соответственно.

В этом случае проще сделать подстановку
,
и решать

Кафедра физхимии ЮФУ (РГУ)
ЧИСЛЕННЫЕ МЕТОДЫ И ПРОГРАММИРОВАНИЕ
Материалы к лекционному курсу
Лектор – ст. преп. Щербаков И.Н.

РЕШЕНИЕ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ

Постановка задачи

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

Обыкновенным дифференциальным уравнением (ОДУ) n-го порядка называется следующее уравнение, которое содержит одну или несколько производных от искомой функции y(x):

Здесь y (n) обозначает производную порядка n некоторой функции y(x), x – это независимая переменная.

В ряде случаев дифференциальное уравнение можно преобразовать к виду, в котором старшая производная выражена в явном виде. Такая форма записи называется уравнением, разрешенным относительно старшей производной (при этом в правой части уравнения старшая производная отсутствует):

Именно такая форма записи принята в качестве стандартной при рассмотрении численных методов решения ОДУ.

Линейным дифференциальным уравнением называется уравнение, линейное относительно функции y(x) и всех ее производных.

Например, ниже приведены линейные ОДУ первого и второго порядков

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

Общее решение ОДУ n -го порядка содержит n произвольных констант C 1 , C 2 , …, C n

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

Так как для решения ДУ n -го порядка необходимо провести n интегрирований, то в общем решении появляется n констант интегрирования.

Частное решение ОДУ получается из общего, если константам интегрирования придать некоторые значения, определив некоторые дополнительные условия, количество которых позволяет вычислить все неопределенные константы интегрирования.

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

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

Например, для дифференциального уравнения первого порядка таблица решения будет представлять собой два столбца – x и y .

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

Или , i = 1, …, N

Для определения частного решения необходимо задать дополнительные условия, которые позволят вычислить константы интегрирования. Причем таких условий должно быть ровно n . Для уравнений первого порядка – одно, для второго - 2 и т.д. В зависимости от способа их задания при решении дифференциальных уравнений существуют три типа задач:

· Задача Коши (начальная задача): Необходимо найти такое частное решение дифференциального уравнения, которое удовлетворяет определенным начальными условиям, заданным в одной точке :

то есть, задано определенное значение независимой переменной (х 0) , и значение функции и всех ее производных вплоть до порядка (n-1) в этой точке. Эта точка (х 0) называется начальной . Например, если решается ДУ 1-го порядка, то начальные условия выражаются в виде пары чисел (x 0 , y 0)

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

· Краевая задача . В этом случае известны значения функции и (или) ее производных в более чем одной точке, например, в начальный и конечный момент времени, и необходимо найти частное решение дифференциального уравнения между этими точками. Сами дополнительные условия в этом случае называются краевыми (граничными ) условиями. Естественно, что краевая задача может решаться для ОДУ не ниже 2-го порядка. Ниже приведен пример ОДУ второго порядка с граничными условиями (заданы значения функции в двух различных точках):

· Задача Штурма-Лиувиля (задача на собственные значения). Задачи этого типа похожи на краевую задачу. При их решении необходимо найти, при каких значениях какого-либо параметра решение ДУ удовлетворяет краевым условиям (собственные значения) и функции, которые являются решением ДУ при каждом значении параметра (собственные функции). Например, многие задачи квантовой механики являются задачами на собственные значения.

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

Рассмотрим некоторые численные методы решения задачи Коши (начальной задачи) обыкновенных дифференциальных уравнений первого порядка. Запишем данное уравнение в общем виде, разрешенном относительно производной (правая часть уравнения не зависит от первой производной):

(6.2)

Необходимо найти значения функции y в заданных точках сетки , если известны начальные значения , где есть значение функции y(x) в начальной точке x 0 .

Преобразуем уравнение умножением на d x

И проинтегрируем левую и правую части между i -ым и i+ 1-ым узлами сетки.

(6.3)

Мы получили выражение для построения решения в i+1 узле интегрирования через значения x и y в i -ом узле сетки. Сложность, однако, заключается в том, что интеграл в правой части есть интеграл от неявно заданной функции, нахождение которого в аналитическом виде в общем случае невозможно. Численные методы решения ОДУ различным способом аппроксимируют (приближают) значение этого интеграла для построения формул численного интегрирования ОДУ.

Из множества разработанных для решения ОДУ первого порядка методов рассмотрим методы , и . Они достаточно просты и дают начальное представление о подходах к решению данной задачи в рамках численного решения .

Метод Эйлера

Исторически первым и наиболее простым способом численного решения задачи Коши для ОДУ первого порядка является метод Эйлера. В его основе лежит аппроксимация производной отношением конечных приращений зависимой (y ) и независимой (x ) переменных между узлами равномерной сетки:

где y i+1 это искомое значение функции в точке x i+1 .

Если теперь преобразовать это уравнение, и учесть равномерность сетки интегрирования, то получится итерационная формула, по которой можно вычислить y i+1 , если известно y i в точке х i :

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

Графическая интерпретация метода Эйлера также не представляет затруднений (см. рисунок ниже). Действительно, исходя из вида решаемого уравнения () следует, что значение есть значение производной функции y(x) в точке x=x i - , и, таким образом, равно тангенсу угла наклона каcательной, проведенной к графику функции y(x) в точке x=x i .

Из прямоугольного треугольника на рисунке можно найти

откуда и получается формула Эйлера. Таким образом, суть метода Эйлера заключается в замене функции y(x) на отрезке интегрирования прямой линией, касательной к графику в точке x=x i . Если искомая функция сильно отличается от линейной на отрезке интегрирования, то погрешность вычисления будет значительной. Ошибка метода Эйлера прямо пропорциональна шагу интегрирования:

Ошибка ~ h

Процесс вычислений строится следующим образом. При заданных начальных условиях x 0 и y 0 можно вычислить

Таким образом, строится таблица значений функции y(x) с определенным шагом (h ) по x на отрезке . Ошибка в определении значения y(x i) при этом будет тем меньше, чем меньше выбрана длина шага h (что определяется точностью формулы интегрирования).

При больших h метод Эйлера весьма неточен. Он дает все более точное приближение при уменьшении шага интегрирования. Если отрезок слишком велик, то каждый участок разбивается на N отрезков интегрирования и к каждому их них применяется формула Эйлера с шагом , то есть шаг интегрирования h берется меньше шага сетки, на которой определяется решение.

Пример:

Используя метод Эйлера, построить приближенное решение для следующей задачи Коши:

На сетке с шагом 0,1 в интервале (6.5)

Решение:

Данное уравнение уже записано в стандартном виде, резрешенном относительно производной искомой функции.

Поэтому, для решаемого уравнения имеем

Примем шаг интегрирования равным шагу сетки h = 0,1. При этом для каждого узла сетки будет вычислено только одно значение (N=1 ). Для первых четырех узлов сетки вычисления будут следующими:

Полные результаты (с точностью до пятого знака после запятой) приведены в в третьей колонке - h =0,1 (N =1). Во второй колонке таблицы для сравнения приведены значения, вычисленные по аналитическому решению данного уравнения .

Во второй части таблицы приведена относительная погрешность полученных решений. Видно, что при h =0,1 погрешность весьма велика, достигая 100% для первого узла x =0,1.

Таблица 1 Решение уравнения методом Эйлера (для колонок указан шаг интегрирования и число отрезков интегрирования N между узлами сетки)

x	Точное решение	0,1	0,05	0,025	0,00625	0,0015625	0,0007813	0,0001953
		1	2	4	16	64	128	512
0	0,000000	0,000000	0,000000	0,000000	0,000000	0,000000	0,000000	0,000000
0,1	0,004837	0,000000	0,002500	0,003688	0,004554	0,004767	0,004802	0,004829
0,2	0,018731	0,010000	0,014506	0,016652	0,018217	0,018603	0,018667	0,018715
0,3	0,040818	0,029000	0,035092	0,037998	0,040121	0,040644	0,040731	0,040797
0,4	0,070320	0,056100	0,063420	0,066920	0,069479	0,070110	0,070215	0,070294
0,5	0,106531	0,090490	0,098737	0,102688	0,105580	0,106294	0,106412	0,106501
0,6	0,148812	0,131441	0,140360	0,144642	0,147779	0,148554	0,148683	0,148779
0,7	0,196585	0,178297	0,187675	0,192186	0,195496	0,196314	0,196449	0,196551
0,8	0,249329	0,230467	0,240127	0,244783	0,248202	0,249048	0,249188	0,249294
0,9	0,306570	0,287420	0,297214	0,301945	0,305423	0,306284	0,306427	0,306534
1	0,367879	0,348678	0,358486	0,363232	0,366727	0,367592	0,367736	0,367844
Относительные погрешности вычисленных значений функции при различных h
x	h	0,1	0,05	0,025	0,00625	0,0015625	0,0007813	0,0001953
	N	1	2	4	16	64	128	512
0,1		100,00%	48,32%	23,76%	5,87%	1,46%	0,73%	0,18%
0,2		46,61%	22,55%	11,10%	2,74%	0,68%	0,34%	0,09%
0,3		28,95%	14,03%	6,91%	1,71%	0,43%	0,21%	0,05%
0,4		20,22%	9,81%	4,83%	1,20%	0,30%	0,15%	0,04%
0,5		15,06%	7,32%	3,61%	0,89%	0,22%	0,11%	0,03%
0,6		11,67%	5,68%	2,80%	0,69%	0,17%	0,09%	0,02%
0,7		9,30%	4,53%	2,24%	0,55%	0,14%	0,07%	0,02%
0,8		7,57%	3,69%	1,82%	0,45%	0,11%	0,06%	0,01%
0,9		6,25%	3,05%	1,51%	0,37%	0,09%	0,05%	0,01%
1		5,22%	2,55%	1,26%	0,31%	0,08%	0,04%	0,01%

Уменьшим шаг интегрирования вдвое, h = 0.05, в этом случае для каждого узла сетки вычисление будет проводиться за два шага (N =2). Так, для первого узла x =0,1 получим:

(6.6)

Данная формула оказывается неявной относительно y i+1 (это значение есть и в левой и в правой части выражения), то есть является уравением относительно y i+1 , решать которое можно, например, численно, применяя какой-либо итерационный метод (в таком виде его можно рассматривать как итерационную формула метода простой итерации). Однако, можно поступить иначи и приблизительно вычислить значение функции в узле i+1 с помощью обычной формулы :

,

которое затем использовать при вычислении по (6.6).

Таким образом получается метод Гюна или метод Эйлера с пересчетом. Для каждого узла интегрирования производится следующая цепочка вычислений

(6.7)

Благодаря более точной формуле интегрирования, погрешность метода Гюна пропорциональна уже квадрату шага интегрирования.

Ошибка ~ h 2

Подход, использованный в методе Гюна, используется для построения так называемых методов прогноза и коррекции , которые будут рассмотрены позже.

Пример:

Проведем вычисления для уравения () с помощью метода Гюна.

При шаге интегрирования h =0,1 в первом узле сетки x 1 получим:

Что намного точнее значения, полученного методом Эйлера при том же шаге интегрирования. В таблице 2 ниже приведены сравнительные результаты вычислений при h = 0,1 методов Эйлера и Гюна.

Таблица 2 Решение уравнения методами Эйлера и Гюна

x	Точное	Метод Гюна		Метод Эйлера
		y	отн. погрешность	y	отн. погрешность
0	0,000000	0,00000		0,00000
0,1	0,004837	0,00500	3,36%	0,00000	100,00%
0,2	0,018731	0,01903	1,57%	0,01000	46,61%
0,3	0,040818	0,04122	0,98%	0,02900	28,95%
0,4	0,070320	0,07080	0,69%	0,05610	20,22%
0,5	0,106531	0,10708	0,51%	0,09049	15,06%
0,6	0,148812	0,14940	0,40%	0,13144	11,67%
0,7	0,196585	0,19721	0,32%	0,17830	9,30%
0,8	0,249329	0,24998	0,26%	0,23047	7,57%
0,9	0,306570	0,30723	0,21%	0,28742	6,25%
1	0,367879	0,36854	0,18%	0,34868	5,22%

Отметим существенное увеличение точности вычислений метода Гюна по сравнению с методом Эйлера. Так, для узла x =0,1 относительное отклонение значения функции, определенного методом Гюна, оказывается в 30 (!) раз меньше. Такая же точность вычислений по формуле Эйлера достигается при числе отрезков интегрирования N примерно 30. Следовательно, при использовании метода Гюна при одинаковой точности вычислений понадобится примерно в 15 раз меньше времени ЭВМ, чем при использовании метода Эйлера.

Проверка устойчивости решения

Решение ОДУ в некоторой точке x i называется устойчивым, если найденное в этой точке значение функции y i мало изменяется при уменьшении шага интегрирования. Для проверки устойчивости, таким образом, надо провести два расчета значения (y i ) – с шагом интегрирования h и при уменьшенной (например, двое) величине шага

В качестве критерия устойчивости можно использовать малость относительного изменения полученного решения при уменьшении шага интегрирования (ε – наперед заданная малая величина)

Такая проверка может осуществляться и для всех решений на всем интервале значений x . Если условие не выполняется, то шаг снова делится пополам и находится новое решение и т.д. до получения устойчивого решения.

Методы Рунге-Кутты

Дальнейшее улучшение точности решения ОДУ первого порядка возможно за счет увеличения точности приближенного вычисления интеграла в выражении .

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

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

Достоинством многошаговых методов Адамса при решении ОДУ заключается в том, что в каждом узле рассчитывается только одно значение правой части ОДУ - функции F(x,y ). К недостаткам можно отнести невозможность старта многошагового метода из единственной начальной точки, так как для вычислений по k -шаговой формуле необходимо знание значения функции в k узлах. Поэтому приходится (k-1) решение в первых узлах x 1 , x 2 , …, x k-1 получать с помощью какого-либо одношагового метода, например метода