Метод Эйлера относиться к численным методам, дающим решение в виде таблицы приближенных значений искомой функции у(х) . Он является сравнительно грубым и применяется в основном для ориентировочных расчетов. Однако идеи, положенные в основу метода Эйлера, являются исходными для ряда других методов.
Рассмотрим дифференциальное уравнение первого порядка
с начальным условием
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++ и представление результатов ее выполнения
Определение дифференциального уравнения Эйлера. Рассмотрены методы его решения.
СодержаниеДифференциальное уравнение Эйлера - это уравнение вида
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)
,
где - многочлены от степеней и ,
соответственно.
В этом случае проще сделать подстановку
,
и решать
Метод Эйлера. Усовершенствованный метод Эйлера.
Классический метод Рунге-Кутты
Не обошла стороной вычислительная математика и дифференциальные уравнения! Сегодня на уроке мы познакомимся с основами приближённых вычислений в этом разделе математического анализа, после чего перед вами приветливо распахнутся толстые-претолстые книги по теме. Ибо вычислительная математика стороной диффуры ещё как не обошла =)
Перечисленные в заголовке методы предназначены для приближённого нахождения решений дифференциальных уравнений , систем ДУ , и краткая постановка наиболее распространённой задачи такова:
Рассмотрим дифференциальное уравнение первого порядка , для которого требуется найти частное решение , соответствующее начальному условию . Что это значит? Это значит, нам нужно найти функцию (предполагается её существование) , которая удовлетворяет данному дифф. уравнению, и график которой проходит через точку .
Но вот незадача – переменные в уравнении разделить невозможно. Никакими известными науке способами. А если и возможно, то получается неберущийся интеграл. Однако частное-то решение существует! И здесь на помощь приходят методы приближенных вычислений, которые позволяют с высокой (а зачастую с высочайшей) точностью «сымитировать» функцию на некотором промежутке.
Идея методов Эйлера и Рунге-Кутты состоит в том, чтобы заменить фрагмент графика ломаной линией , и сейчас мы узнаем, как эта идея реализуется на практике. И не только узнаем, но и непосредственно реализуем =) Начнём с исторически первого и самого простого метода. …Вы хотите иметь дело со сложным дифференциальным уравнением? Вот и я тоже не хочу:)
Задание
Найти частное решение дифференциального уравнения , соответствующее начальному условию , методом Эйлера на отрезке с шагом . Построить таблицу и график приближённого решения.
Разбираемся. Во-первых, перед нами обычное линейное уравнение , которое можно решить стандартными способами, и поэтому очень трудно устоять перед соблазном сразу же найти точное решение:
– желающие могут выполнить проверку и убедиться, что данная функция удовлетворяет начальному условию и является корнем уравнения .
Что нужно сделать? Нужно найти и построить ломаную
, которая приближает график функции 
на промежутке . Поскольку длина этого промежутка равна единице, а шаг составляет , то наша ломаная
будет состоять из 10 отрезков:
причём, точка 
уже известна – она соответствует начальному условию . Кроме того, очевидны «иксовые» координаты других точек:
Осталось найти 
. Никакого дифференцирования
и интегрирования
– только сложение и умножение! Каждое следующее «игрековое» значение получается из предыдущего по простой рекуррентной
формуле:
Представим дифференциальное уравнение в виде :
Таким образом: 
«Раскручиваемся» от начального условия :
Понеслось:
Результаты вычислений удобно заносить в таблицу:
А сами вычисления автоматизировать в Экселе – потому что в математике важен не только победный, но ещё и быстрый конец:)
По результатам 2-го и 3-го столбцов изобразим на чертеже 11 точек и 10 отрезков, соединяющих смежные точки. Для сравнения я построю график точного частного решения 
:
Существенным недостатком простого метода Эйлера является слишком большая погрешность, при этом легко заметить, что погрешность имеет тенденцию накапливаться – чем дальше мы уходим от точки , тем преимущественно
больше становится расхождение между приближением и истиной. Это объяснимо самим принципом, который Эйлер положил в основу своего метода: отрезки параллельны соответствующим
касательным
к графику функции в точках . Данный факт, кстати, тоже хорошо просматривается по чертежу.
Как можно улучшить приближение? Первая мысль – измельчить разбиение. Разделим отрезок , например, на 20 частей. Тогда шаг составит: 
, и совершенно понятно, что ломаная из 20 звеньев заметно точнее приблизит частное решение. С помощью того же Экселя не составит труда обработать 100-1000 и даже миллион (!) промежуточных отрезков, однако зададимся вопросом: а нельзя ли КАЧЕСТВЕННО улучшить метод?
Но перед тем как раскрыть этот вопрос, не могу не остановиться на неоднократно прозвучавшей сегодня фамилии. Читая биографию Леонарда Эйлера , просто поражаешься, как невероятно много может успеть сделать за свою жизнь человек! Сопоставимо вспомнился только К.Ф. Гаусс. …Вот и мы постараемся не потерять мотивацию к обучению и новым открытиям:))
Усовершенствованный метод Эйлера
Рассмотрим тот же самый пример: дифференциальное уравнение , частное решение, удовлетворяющее условию , промежуток и его разбиение на 10 частей
( – длина каждой части).
Цель усовершенствования состоит в том, чтобы приблизить «красные квадратики» ломаной к соответствующим «зелёным точкам» точного решения 
.
И идея модификации такова: отрезки должны быть параллельны касательным , которые проведены к графику функции не на левых краях , а «посерединке» интервалов разбиения. Что, естественно, улучшит качество приближения.
Алгоритм решения работает в том же русле, но формула, как нетрудно догадаться, усложняется:
, где 
Плясать вновь начинаем от частного решения и сразу же находим 1-й аргумент «внешней» функции:
Теперь находим нашего «монстра», который на поверку оказался не таким уж и страшным – обратите внимание, что это ТА ЖЕ функция 
, вычисленная в другой точке:
Умножаем результат на шаг разбиения:
Таким образом:
Алгоритм заходит на второй круг, не поленюсь, распишу его подробно:
рассматриваем пару и находим 1-й аргумент «внешней» функции:
Рассчитываем и находим её 2-й аргумент:
Вычислим значение:
и его произведение на шаг:
Вычисления разумно провести в Экселе (растиражировав формулы по той же схеме – см. видеоролик выше)
, а результаты свести в таблицу:
Числа целесообразно округлять до 4-5-6 знаков после запятой. Нередко в условии той или иной задачи есть прямое указание
, с какой точностью следует проводить округление. Я подровнял сильно «хвостатые» значения до 6 знаков.
По результатам 2-го и 3-го столбцов (слева)
построим ломаную
, и для сравнения я снова приведу график точного решения 
:
Результат существенно улучшился! – красные квадратики практически «спрятались» за зелёными точками точного решения.
Однако нет пределов совершенству. Одна голова хорошо, а две – лучше. И снова немецкие:
Классический метод Рунге-Кутты 4-го порядка
Его цель добиться ещё бОльшего приближения «красных квадратиков» к «зелёным точкам». Вы спросите, куда ещё ближе? Во многих, в частности физических, исследованиях бывает ПРИНЦИПИАЛЬНО важен 10-й, а то и 50-й точный знак после запятой. Нет, такой точности можно достичь и простым методом Эйлера, но на СКОЛЬКО частей придётся разбить промежуток ?! …Хотя с современными вычислительными мощностями это не проблема – тысячи кочегаров китайского космического корабля гарантируют!
И, как правильно подсказывает заголовок, при использовании метода Рунге-Кутты на каждом шаге
нам придётся вычислить значение функции 
4 раза
(в отличие от двукратного вычисления в предыдущем параграфе)
. Но задача эта вполне и вполне подъёмная если нанять китайцев. Каждое следующее «игрековое» значение получается из предыдущего – ловим формулы:
, где 
, где:
Готовы? Ну тогда начинаем:))
Таким образом:
Первая строка запрограммирована, и я копирую формулы по образцу:
Не думал, что так быстро разделаюсь с методом Рунге-Кутты =)
В чертеже нет смысла, поскольку он уже не показателен. Давайте лучше проведём аналитическое сравнение точности
трёх методов, ибо когда известно точное решение 
, то грех не сравнить. Значения функции в узловых точках элементарно рассчитываются в том же Экселе – один раз забиваем формулу и тиражируем её на остальные .
В нижеследующую таблицу я сведу значения (для каждого из трёх методов) и соответствующие абсолютные погрешности
приближённых вычислений:
Как видите, метод Рунге-Кутты даёт уже 4-5 верных знака после запятой по сравнению с 2 верными знаками усовершенствованного метода Эйлера! И это не случайность:
– Погрешность «обычного» метода Эйлера не превосходит шага разбиения. И в самом деле – взгляните на самый левый столбец погрешностей – там после запятых только один ноль, что и говорит нам о точности 0,1.
– Усовершенствованный метод Эйлера гарантирует точность: 
(смотрим на 2 нуля после запятой в средней колонке погрешностей)
.
– И, наконец, классический метод Рунге-Кутты обеспечивает точность 
.
Изложенные оценки погрешностей строго обосновывается в теории.
Как можно ЕЩЁ улучшить точность приближения? Ответ прямо-таки философский: качеством и/или количеством =) В частности, существует и другие, более точные модификации метода Рунге-Кутты. Количественный путь, как уже отмечалось, состоит в уменьшении шага, т.е. в разбиении отрезка на бОльшее количество промежуточных отрезков. И с увеличением этого количества ломаная 
всё больше и больше будет походить на график точного решения 
и в пределе
– совпадёт с ним.
В математике это свойство называется спрямляемостью кривой . К слову (небольшой оффтоп) , «спрямить» удаётся далеко не всё – рекомендую прочитать интереснейшую , в которых уменьшение «участка исследования» не влечёт за собой упрощение объекта исследования.
Так получилось, что я разобрал всего лишь одно дифференциальное уравнение и поэтому пара дополнительных замечаний. Что ещё нужно иметь в виду на практике? В условии задачи вам может быть предложен другой отрезок и другое разбиение, причём иногда встречается следующая формулировка: «найти методом… …на промежутке , разбив его на 5 частей». В этом случае нужно найти шаг разбиения 
, после чего придерживаться обычной схемы решения. Кстати, начальное условие должно быть такого вида: , то есть «икс нулевое», как правило, совпадает с левым концом отрезка. Образно говоря, ломаная всегда «выходит» из точки .
Безусловным достоинством рассмотренных методов, является тот факт, что они применимы к уравнениям с очень сложной правой частью. И безусловный недостаток – далеко не каждый диффур можно представить в таком виде.
Но почти всё в этой жизни поправимо! – ведь мы рассмотрели лишь малую толику темы, и моя фраза о толстых-претолстых книгах была вовсе не шуткой. Существует великое множество приближённых методов нахождения решений ДУ и их систем, в которых применяются, в том числе, принципиально другие подходы. Так, например, частное решение можно приблизить степенным рядом . Однако это уже статья другого раздела.
Надеюсь, мне удалось разнообразить скучноватую вычислительную математику, и вам было интересно!
Спасибо за внимание!
Пусть нам известна входная динамическая последовательность X (входной сигнал) и модель (способ преобразования входного сигнала в выходной сигнал). Рассматривается задача определения выходного сигнала y (t ) (см. рис. 10.1 ).
Модель динамической системы может быть представлена дифференциальным уравнением. Основное уравнение динамики:
y " = f (x (t ), y (t ), t ) .
Известны начальные условия в нулевой момент времени t 0 : y (t 0) , x (t 0) . Чтобы определить выходной сигнал, заметим, что по определению производной:
Нам известно положение системы в точке «1», требуется определить положение системы в точке «2». Точки отделены друг от друга расстоянием Δt (рис. 10.2 ). То есть расчет поведения системы производится по шагам. Из точки «1» мы скачком (дискретно) переходим в точку «2», расстояние между точками по оси t называется шагом расчета Δt .
методом Эйлера на одном шаге
Последняя формула называется формулой Эйлера .
Очевидно, чтобы узнать состояние системы в будущем y (t + Δt ) , надо к настоящему состоянию системы y (t ) прибавить изменение Δy , прошедшее за время Δt .
Рассмотрим еще раз это важное соотношение, выведя его из геометрических соображений (рис. 10.3 ).
Пусть A точка, в которой состояние системы известно. Это «настоящее» состояние системы.
В точке A к траектории движения системы проведем касательную. Касательная это производная функции f (x (t ), y (t ), t ) по переменной t . Производную в точке всегда легко вычислить, достаточно подставить известные переменные (в момент «Настоящее» они известны) в формулу y " = f (x (t ), y (t ), t ) .
Заметим, что, по определению, производная связана с углом наклона касательной: y " = tg(α ) , значит, угол α легко вычислить (α = arctg(y " ) ) и провести касательную.
Проведем касательную до пересечения с линией t + Δt . Момент t + Δt соответствует «будущему» состоянию системы. Проведем линию параллельно оси t от точки A до пересечения с линией t + Δt . Линии образуют прямоугольный треугольник ABC, один катет которого равен Δt (известен). Известен также угол α . Тогда второй катет в прямоугольном треугольнике ABC равен: a = Δt · tg(α ) . Теперь легко вычислить ординату точки B. Она состоит из двух отрезков y (t ) и a . Ордината символизирует положение системы в точке y (t + Δt ) . То есть y (t + Δt ) = y (t ) + a или далее y (t + Δt ) = y (t ) + Δt · tg(α ) или, подставляя дальше, имеем: y (t + Δt ) = y (t ) + Δt · y " и, наконец, y (t + Δt ) = y (t ) + Δt · f (x (t ), y (t ), t ) . Снова мы получили формулу Эйлера (из геометрических соображений).
Эта формула может дать точные результаты только при очень малых Δt (говорят при Δt > 0 ). При Δt ≠0 формула дает расхождение между истинным значением y и расчетным, равное ε , поэтому в ней должен стоять знак приближенного равенства, либо она должна быть записана так:
y (t + Δt ) = y (t ) + Δt · f (x (t ), y (t ), t ) + ε .
И в самом деле. Взгляните еще раз на рис. 10.3 . Будем мысленно сдвигать линию t + Δt влево (фактически, будем приближать значение Δt к нулю). Как нетрудно видеть, расстояние BB * = ε , то есть ошибка! будет сокращаться. В пределе (при Δt > 0 ) значение ошибки ε будет равно нулю.
Итак, заменяя реальную кривую прямой (касательной) на отрезке Δt , мы вносим в решение ошибку, попадая в результате не в точку «2» (см. рис. 10.2 ), а рядом, в точку «3». Очевидно, что этот численный метод на каждом шаге имеет погрешность расчета ε .
Из рисунка видно, что чем меньше взять величину Δt , тем меньше будет ошибка расчета ε . То есть для расчета поведения системы на сколько-нибудь продолжительном отрезке времени (например, от t 0 до t k ), чтобы уменьшить ошибку на каждом шаге, шаги Δt делают по возможности малыми. Для достижения точки t k отрезок (t k t 0) делится на отрезки длиной Δt ; таким образом, всего получится N = (t k t 0)/Δt шагов. В результате расчета придется формулу Эйлера применить для каждого шага, то есть N раз. Но следует иметь в виду, что ошибки ε i на каждом i -ом шаге (в простейшем случае) складываются, а общая ошибка быстро накапливается (см. рис. 10.4 ). И в этом состоит существенный недостаток данного метода. Хотя с помощью этого метода можно получить (в численном виде) решение любого дифференциального уравнения (в том числе и неразрешимого аналитически). Уменьшая шаг, мы получаем более точные решения, но при этом не следует забывать, что увеличение числа шагов ведет к вычислительным затратам и снижению быстродействия. Кроме того, при большом числе итераций в расчет вносится другая существенная погрешность из-за ограниченной точности вычислительных машин и ошибок округления.
Задача 1 . Дано дифференциальное уравнение y " = 2t y . Задано начальное положение системы: y (0) = 1 . Требуется найти y (t ) , то есть поведение системы на интервале времени t от 0 до 1.
Аналитический способ решения задачи 1
y " = 2t y .
Методом разделения переменных найдем:
y " /y = 2t
Будем интегрировать от 0 до t i , тогда согласно правилам интегрирования имеем:
Полученное аналитическое решение характеризуется тем, что оно является абсолютно точным, но если уравнение окажется сколько-нибудь сложным, то решение не будет найдено вовсе. Аналитический путь решения не универсален.
Численный способ решения задачи 1
Численный способ решения предполагает, что расчет будет вестись по формуле Эйлера на ряде последовательных шагов. На каждом шаге решение имеет свою ошибку (см. рис. 10.2 ), поскольку на каждом шаге кривая заменяется прямым отрезком.
При алгоритмической реализации расчет реализуется циклом, в котором изменяется t (счетчик t ) и y :
Блок-схема при реализации метода на компьютере показана на рис. 10.5 .
В реализации Стратум запись будет выглядеть так (наличие символа «~» при t ):
Будем искать значение y рассмотренного ранее примера в численном виде на промежутке от T = 0 до T = 1 . Возьмем число шагов n = 10 , тогда шаг приращения Δt составит: Δt = (1 0)/n = (1 0)/10 = 0.1 .
| Таблица 10.1. Численный расчет уравнения методом Эйлера и сравнение результата с точным решением на каждом шаге |
|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| i | t i | y i = y i 1 + y" i 1 · Δt | y" i = 2t i · y i | Δy i = y" i · Δt | y i + 1 = y i + Δy i | y точн. = exp(t i 2) |
| 0 | 0.0 | 1 | 0 | 0 | 1 | 1 |
| 1 | 0.1 | 1 | 0.2 | 0.02 | 1.02 | 1.0101 |
| 2 | 0.2 | 1.02 | 0.408 | 0.0408 | 1.0608 | 1.0408 |
| 3 | 0.3 | 1.061 | 0.636 | 0.0636 | 1.1246 | 1.0942 |
| 4 | 0.4 | 1.124 | 0.900 | 0.0900 | 1.2140 | 1.1735 |
| 5 | 0.5 | 1.214 | 1.214 | 0.1214 | 1.3354 | 1.2840 |
| 6 | 0.6 | 1.336 | 1.603 | 0.1603 | 1.4963 | 1.4333 |
| 7 | 0.7 | 1.496 | 2.095 | 0.2095 | 1.7055 | 1.6323 |
| 8 | 0.8 | 1.706 | 2.729 | 0.2729 | 1.9789 | 1.8965 |
| 9 | 0.9 | 1.979 | 3.561 | 0.3561 | 2.3351 | 2.2479 |
| 10 | 1.0 | 2.335 | 4.669 | 0.4669 | 2.8019 | 2.7183 |
Обратите внимание на то, что рассчитанное численно значение (y i + 1 ) отличается от точного (y точн. ), и погрешность (разница столбцов y i + 1 и y точн. ) в процессе расчета нарастает подобно тому, как было показано на рис. 10.4 .
Теперь подсчитаем относительную погрешность σ для расчетного значения y (1) , полученного численно, в сравнении с теоретическим точным y теор. по следующей формуле:
σ = (1 y расч. /y теор.) · 100%
и сравним σ при различных значениях Δt .
Если будем менять значение шага Δt , например, уменьшать шаг, то относительная погрешность расчета тоже будет уменьшаться. Вот что получится при вычислении значения y (1) с разными значениями шага (см. табл. 10.2).
| Таблица 10.2. Зависимость погрешности расчета от размера шага Δt |
|||||||||||||||
| Δt | y расч. (1) | y теор. (1) | σ |
| 1/10 | 2.3346 | 2.7183 | 14% |
| 1/20 | 2.5107 | 2.7183 | 8% |
| 1/100 | 2.6738 | 2.7183 | 2% |
Как видим, с уменьшением шага приращения Δt уменьшается величина относительной погрешности, а значит, повышается точность расчета.
Обратите внимание, что изменение шага в 10 раз (с 1/10 до 1/100) ведет к изменению величины ошибки примерно тоже в 10 раз (с 14% до 2%). При изменении шага в 100 раз ошибка примерно уменьшится тоже в 100 раз. Иными словами размер шага и ошибка для метода Эйлера связаны линейно. Хотите уменьшить в 10 раз ошибку уменьшайте в 10 раз шаг и увеличивайте соответственно в 10 раз количество вычислений. Этот факт в математике принято обозначать символом ε = O (Δt ) , а метод Эйлера называют методом первого порядка точности.
Поскольку в методе Эйлера ошибка достаточно велика и от шага к шагу накапливается, а точность пропорциональна количеству вычислений, то метод Эйлера обычно применяют для грубых расчетов, для оценки поведения системы в принципе. Для точных количественных расчетов применяют более точные методы.
Примечания
- Каждый численный метод обладает точностью, поскольку результат отличается от теоретического. Точность метода зависит от величины шага. Различные методы имеют различную точность. Порядок зависимости точности от величины шага обозначают как O (h ) . У метода Эйлера первый порядок точности, зависимость ошибки от величины шага линейна.
- Если при уменьшении шага предел y n стремится к значению y теор. , то говорят, что метод обладает сходимостью. Исследователей интересует скорость сходимости метода.
- Метод должен быть устойчив. Устойчивость связана с некоторой критической величиной шага. При проявлении неустойчивости наблюдается полное искажение качественной картины расчета, «разболтка» результата.
- При выборе метода рекомендуется сначала добиться устойчивости, а внутри области устойчивости сходимости результата. Устойчивость обеспечивает качественную картину. Сходимость обеспечивает количественный результат (см. также рис. 10.10 ).
Изложенное в пп. 1-4 поясним на примере.
Пример . Пусть
Качественно это уравнения описывают процесс теплообмена между двумя телами, температуры которых в некоторый момент времени обозначим как A и B . Вообще A и B переменные, меняющиеся во времени t . Найти поведение системы означает, что надо найти, как будут меняться температуры A (t ) и B (t ) .
Интуитивно ясно, что при начальной разнице температур A = 8 и B = 5 температуры тел постепенно со временем должны выровняться, так как более горячее тело будет отдавать энергию более холодному, и его температура будет уменьшаться, а более холодное тело будет принимать энергию от более горячего, и его температура будет увеличиваться. Процесс теплообмена закончится (то есть изменения прекратятся) тогда, когда температуры двух тел станут одинаковыми.
Проведем несколько расчетов поведения A (t ) и B (t ) с разной величиной шага Δt .
Будем брать различную величину шага Δt и находить соответствующие значения A и B во времени по следующим формулам Эйлера:
A
нов. = A
пред. + (B
пред. A
пред.) · Δt
,
B
нов. = B
пред. + (A
пред. B
пред.) · Δt
.
Расчет при Δt = 2 (табл. 10.3).
Наблюдается явление «разболтки» (см. рис. 10.6 ). Неустойчивое решение. Из физических соображений очевидно, что так вести себя два тела в процессе теплообмена не могут.
 |
неверно. Решение неустойчиво
Расчет при Δt = 1 (табл. 10.4).
| Таблица 10.4. Изменение температур тел при численном расчете с шагом 1 |
|||||||||||||||
| № шага |
t | A | B |
| 0 | 0 | 8 | 5 |
| 1 | 1 | 5 | 8 |
| 2 | 2 | 8 | 5 |
Наблюдается поведение решения системы на границе устойчивости (см. рис. 10.7 ).
 |
неверно. Решение находится на грани устойчивости
Расчет при Δt = 0.5 (табл. 10.5).
| Таблица 10.5. Изменение температур тел при численном расчете с шагом 0.5 |
|||||||||||||||
| № шага |
t | A | B |
| 0 | 0 | 8 | 5 |
| 1 | 0.5 | 6.5 | 6.5 |
| 2 | 1.0 | 6.5 | 6.5 |
Решение устойчиво, соответствует правильной качественной картине (см. рис. 10.8 ). Температуры тел постепенно сближаются, становятся со временем одинаковыми. Но решение пока имеет большую погрешность.
Решение (поведение системы) имеет большую погрешность
Расчет при Δt = 0.1 (табл. 10.6).
| Таблица 10.6. Изменение температур тел при численном расчете с шагом 0.1 |
|||||||||||||||||||||||||||
| № шага |
t | A | B |
| 0 | 0 | 8 | 5 |
| 1 | 0.1 | 7.7 | 5.3 |
| 2 | 0.2 | 7.46 | 5.54 |
| 3 | 0.3 | 7.27 | 5.73 |
| 4 | 0.4 | 7.12 | 5.88 |
| 5 | 0.5 | 7.00 | 6.00 |
Решение устойчиво. Решение более точно (см. рис. 10.9 ).
Количественно решение более точно
Роль изменения величины шага иллюстрирует рис. 10.10 .
Системой дифференциальных уравнений называется система вида
где x - независимый аргумент,
y i - зависимая функция, ,
y i | x=x0 =y i0 - начальные условия.
Функции y i (x), при подстановке которой система уравнений обращается в тождество, называется решением системой дифференциальных уравнений .
Численные методы решения систем дифференциальных уравнений.
Дифференциальным уравнением второго порядка называется уравнение вида
Функция y(x), при подстановке которой уравнение обращается в тождество, называется решением дифференциального уравнения .
Численно ищется частное решение уравнения (2), которое удовлетворяет заданным начальным условиям, то есть решается задача Коши.
Для численного решения дифференциальное уравнение второго порядка преобразуется в систему двух дифференциальных уравнений первого порядка и приводится к машинному виду (3). Для этого вводится новая неизвестная функция , слева в каждом уравнении системы оставляют только первые производные неизвестных функций, а в правых частях производных быть не должно
 .
| (3) |
Функция f 2 (x, y 1 , y) в систему (3) введена формально для того, чтобы методы, которые будут показаны ниже, могли быть использованы для решения произвольной системы дифференциальных уравнений первого порядка. Рассмотрим несколько численных методов решения системы (3). Расчетные зависимости для i+1 шага интегрирования имеют следующий вид. Для решения системы из n уравнений расчетные формулы приведены выше. Для решения системы из двух уравнений расчетные формулы удобно записать без двойных индексов в следующем виде:
- Метод Эйлера
.
у 1,i+1 =у 1,i +hf 1 (x i , y 1,i , y i),
у i+1 =у i +hf 2 (x i , y 1,i , y i),
- Метод
Рунге-Кутта четвертого порядка
.
у 1,i+1 =у 1,i +(m 1 +2m 2 +2m 3 +m 4)/6,
у i+1 =у i +(k 1 +2k 2 +2k 3 +k 4)/6,
m 1 =hf 1 (x i , y 1,i , y i),
k 1 =hf 2 (x i , y 1,i , y i),
m 2 =hf 1 (x i +h/2, y 1,i +m 1 /2, y i +k 1 /2),
k 2 =hf 2 (x i +h/2, y 1,i +m 1 /2, y i +k 1 /2),
m 3 =hf 1 (x i +h/2, y 1,i +m 2 /2, y i +k 2 /2),
k 3 =hf 2 (x i +h/2, y 1,i +m 2 /2, y i +k 2 /2),
m 4 =hf 1 (x i +h, y 1,i +m 3 , y i +k 3),
k 4 =hf 2 (x i +h, y 1,i +m 3 , y i +k 3),
где h - шаг интегрирования. Начальные условия при численном интегрировании учитываются на нулевом шаге: i=0, x=x 0 , y 1 =y 10 , y=y 0 .
Контрольное задание по зачетной работе.
Колебания с одной степенью свободы
Цель. Изучение численных методов решения дифференциальных уравнений второго порядка и систем дифференциальных уравнений первого порядка.
Задание. Численно и аналитически найти:
- закон движения материальной точки на пружинке х(t),
- закон изменения силы тока I(t) в колебательном контуре (RLC - цепи) для заданных в табл.1,2 режимов. Построить графики искомых функций.
Варианты заданий.
Таблица режимов
Варианты заданий и номера режимов:
- движение точки
- RLC - цепь
Рассмотрим более подробно порядок составления дифференциальных уравнений и приведения их к машинному виду для описания движения тела на пружинке и RLC-цепи.
- Название, цель работы и задание.
- Математическое описание, алгоритм (структограмма) и текст программы.
- Шесть графиков зависимости (три точные и три приближенные) x(t) или I(t), выводы по работе.
