В.А. Антонюк, А.П. Иванов - Программирование и информатика (Краткий конспект лекций) (1109543), страница 7
Текст из файла (страница 7)
Требуется найти (приближённо), чему равна функцияв других точках , = 1, . . . , заданного отрезка (т.е., на заданной сетке): 0 < 1 < . . . < ≤ . Довольно часто сетка выбирается равномерной с шагом ℎ, т.е. +1 − = ℎ, =0, . . . , − 1. Далее рассматриваются некоторые методы приближённого решения (иногдаговорят – интегрирования) ОДУ первого порядка: метод Эйлера (метод ломаных), методсредней точки (модифицированный метод Эйлера), метод «предиктор-корректор»(метод Эйлера с пересчётом), классический метод Рунге-Кутты.7.3.Ошибки приближённых методовназывается ошибка на одном шаге метода.
Глобальная ошибка– суммарная ошибка за все шаги. Если глобальная ошибка при уменьшении шага уменьшается примерно пропорционально уменьшению шага, то говорят, что это – метод первогопорядка. Методы второго порядка имеют глобальную ошибку, уменьшающуюся пропорционально второй степени уменьшения шага и т.д.Локальной ошибкойПро рассматриваемые методы приближённого нахождения решений дифференциальныхуравнений можно сказать следующее: метод Эйлера – первого порядка, метод средней точки и метод «предиктор-корректор» – методы второго порядка, классический метод РунгеКутты – метод четвёртого порядка.7.4.Устойчивость приближённого решенияЧисленный метод интегрирования ОДУ называется устойчивым (при заданном шаге ℎи для определённого ), если результат численного интегрирования уравнения ′ () = ()остаётся ограниченным при → ∞.Устойчивость зависит от шага ℎ и значения .
На практике устойчивость зависит от произведения ℎ. Например, для метода Эйлера (как выяснится далее) +1 = (1 + ℎ) , а потому,соответственно, +1 = (1 + ℎ) 0 , значит, метод будет устойчив при |1 + ℎ| ≤ 1, т.е., когда−2 ≤ ℎ ≤ 0.317.5.Метод Эйлера (метод ломаных)Подстановка начального условия в исходное ОДУ даёт значение производной функции = () в начальной точке0′ = (0 , 0 )Значение в следующей точке сетки можно оценитьпо значению этой производной и шагу сетки1 ≃ 0 + (1 − 0 ) (0 , 0 )Поэтому общая формула метода Эйлера (на -ом шаге) для равномерной сетки, очевидно, будет такова:+1 ≃ + ℎ ( , ), = 0, .
. . , − 1.Метод Эйлера является одношаговым, поскольку используется значение функции толькоиз предыдущего шага. Локальная ошибка метода ∼ ℎ2 , глобальная (суммарная) погрешность ∼ ℎ, поэтому говорят, что этот метод имеет первый порядок точности.Если бы правая часть исходного ОДУ не зависела от , то точное значение искомой функциив точке +1 = + ℎ определялось бы интегралом∫︁ +ℎ+1 = + () а это значит, что общая формула метода Эйлера – это аналог формулы численного интегрирования методом левых прямоугольников.Этот метод из-за быстрого накопления глобальной ошибки на практике не используется,но удобен для теоретических оценок.7.6.Метод средней точки (модифицированный метод Эйлера)Оценивается значение в «половинной» точке + 1 по2методу Эйлера:ℎ+ 1 ≃ + ( , )22Затем в этой половинной точке вычисляется угловой коэффициент касательной, примерно равный (+ 1 , + 1 ), и в этом направлении совершается пе22реход из точки в точку +1 , чтобы получить новоеприближённое значение искомой функции:+1 ≃ + ℎ (+ 1 , + 1 )2Этот метод является методом второго порядка точности.3227.7.Метод «предиктор-корректор» (метод Эйлера с пересчётом)На первом этапе (прогноз) предсказывается ˜+1 пометоду Эйлера:˜+1 ≃ + ℎ ( , )На втором этапе (коррекция) значение корректируется путём усреднения угловых коэффициентов вточках ( , ) и (+1 , ˜+1 ):+1 ≃ + ℎ ( , ) + (+1 , ˜+1 )2Т.е., сначала мы «грубо» вычисляем значение функции ˜+1 и наклон интегральной кривой (+1 , ˜+1 ) вновой точке.
Затем находим средний наклон на этомшаге и по нему корректируем значение +1 в новойточке. За счёт коррекции точность метода повышается по сравнению с методом Эйлера,так что этот метод тоже является методом второго порядка точности.7.8.Классический метод Рунге-КуттыОписывается системой соотношений для 0, . . . , − 1:=+1 = + ℎ6 (k1 + 2k2 + 2k3 + k4 ),k1 = ( , ),k2 = ( + ℎ2 , + ℎ2 k1 ),k3 = ( + ℎ2 , + ℎ2 k2 ),k4 = ( + ℎ, + ℎk3 ).Последовательно вычисляются приближающие угловые коэффициенты: k1 – в исходной точке, k2 – наполовинном шаге (как в методе средней точки), k3– тоже на половинном шаге, но по уточнённому значению углового коэффициента k2 вместо k1 , k4 – нацелом шаге по предыдущему значению k3 .
Стоит отметить, что получаемые здесь на каждом этапе k1 , k2 , k3 , k4 – угловые коэффициенты четырёх разных интегральных кривых в трёх точках: , + ℎ2 , +1 = + ℎ. Для получениянового значения искомой функции на полном шаге используется взвешенное среднее этихугловых коэффициентов.Классический метод Рунге-Кутты имеет четвёртый порядок точности.
Он является явными допускает расчёт с переменным шагом.В случае, когда правая часть исходного ОДУ не зависит от , легко заметить, что точноезначение искомой функции в точке +1 = + ℎ, определяемое интегралом∫︁ +ℎ+1 = + () аппроксимируется в соответствии с формулой Симпсона, а это значит, что формула классического метода Рунге-Кутты – это аналог формулы численного интегрирования Симпсона.33На случай систем дифференциальных уравнений различные вариации метода Рунге-Куттыпереносятся при помощи формальной замены скалярных величин на нужное число векторных компонент.7.9.Алгоритм ВерлеАлгоритм численного интегрирования Верле (переоткрыт L.Verlet, 1967) часто используется в молекулярной динамике для вычисления траекторий частиц.
Он более устойчив, чемпростой метод Эйлера. В основной форме алгоритма Верле вычисляется следующееместоположение частицы по текущему и прошлому положениям, причём без использованияскорости. Для получения формулы запишем разложение в ряд Тейлора вектора местоположения ⃗() частицы в моменты времени + ∆ и − ∆.⃗( + ∆) = ⃗() + ⃗ ()∆ + ⃗()(∆)2 /2 + .
. .⃗( − ∆) = ⃗() − ⃗ ()∆ + ⃗()(∆)2 /2 − . . .Здесь ⃗ – положение частицы, ⃗ – её скорость, ⃗ – ускорение. Сложив эти два равенства ивыражая ⃗( + ∆), получим:⃗( + ∆) ≈ 2⃗() − ⃗( − ∆) + ⃗()(∆)2Ускорение частицы известно, поскольку пропорционально силе, действующей на частицу, иобратно пропорционально её массе. Таким образом, действительно, положение вычисляетсябез знания скорости. Однако, алгоритм не является, как говорят, самостартующим : сего помощью нельзя в самом начале найти текущее положение, зная только прошлое, здесьнужно воспользоваться другим способом.Если скорость всё-таки нужна (например, для определения кинетической энергии частицы), можно воспользоваться скоростной формой алгоритма Верле :1. Вычисляется скорость на «половинном» шаге: ⃗ ( + 12 ∆) = ⃗ () + 12 ⃗()∆2.
Вычисляется положение на следующем шаге: ⃗( + ∆) = ⃗() + ⃗ ( + 21 ∆)∆3. Вычисляется ⃗( + ∆) в положении ⃗( + ∆) – по действующей силе4. Вычисляется скорость на следующем шаге: ⃗ ( + ∆) = ⃗ () + 12 (⃗() + ⃗( + ∆))∆Или, приводя всё к двум формулам:1⃗( + ∆) = ⃗() + ⃗ ()∆ + ⃗()(∆)221⃗ ( + ∆) = ⃗ () + (⃗() + ⃗( + ∆))∆2348.ИнтерполяцияИнтерполяцией называют такую разновидность аппроксимации, при которой кривая построенной функции проходит точно через имеющиеся точки данных. Часто нас интересуетаппроксимация сложной функции другой, более простой в вычислительном смысле. Например, полиномы характерны тем, что легко вычисляются и дифференцируются.Рассмотрим в качестве примера важный частный случай – линейную интерполяцию:необходимо для функции (), заданной в точках 1 и 2 , построить приближение с помощью линейной функции + .Для любой точки (, ) прямой, проходящей через точки (1 , (1 )) и (2 , (2 )) справедливо соотношение: − (1 ) (2 ) − (1 )=, − 1 2 − 1откуда можно выразить зависимость от :(1 )· ( − 1 ) = (1 ) + (22)−−1 (2 )− (1 )= (1 ) + 2 −1 · − ( (2 ) − (1 )) ·(1 )1 (2 )= (22)−· + 2 (12)−.−1−112 −1Таким образом,где = (2 )− (1 )2 −1и= () ≈ + ,2 (1 )−1 (2 ).2 −1Можно также записать эту формулу и по-другому, перегруппировывая слагаемые так, чтобы образовалась «взвешенная» сумма значений функции в обеих точках с соответствующими сомножителями: () ≈===(1 )· ( − 1 ) (1 ) + (22)−−111 (1 ) −[︁ (1 ) · −+ (2 ) · −−12 ]︁2 −111 (1 ) · 1 − −+ (2 ) · −2 −12 −1−1+ (2 ) · −.
(1 ) · 22−12 −1Обратите внимание на поведение сомножителей при значениях функции (1 ), (2 ) : ониравны единице для той точки, в которой вычисляется функция, и нулю – для другой точки.8.1.Общая постановка задачи интерполяцииПусть даны точки , = 0, 1, ..., (их называют узлами интерполяции), в которых известны значения функции () : = ( ). Надо найти интерполирующую функцию (),такую, что ( ) = для = 0, 1, ..., , а сама – из заданного класса функций (чащевсего применяется интерполяция полиномами).8.2.Интерполяционные многочлены ЛагранжаМногочлены минимальной степени, приближающие данные значения на заданном набореточек.35Пусть даны пары чисел (0 , 0 ), (1 , 1 ), ..., ( , ), где все – различны.
Тогда существуетединственный многочлен () степени не более , для которого ( ) = , = 0, 1, ..., .Попытаемся построить базисные функции 0 (), 1 (), ..., (), такие, что{︃ ( ) =1, = 0, ≠ Т.е., -ая базисная функция равна нулю во всех узлах интерполяции, кроме -ого, где онаравна единице.Имея такие базисные функции, мы можем сконструировать интерполирующую функциютак:() =∑︁ ()=0Рассмотрим для иллюстрации простой пример (см.
рис.7): известны значения функции в точках 1, 2, 4; они равны 1, 1/3 и 1/5 соответственно. Требуется найти многочлен, проходящийчерез эти точки.Таким образом, известна такая таблица значений функции:11241315Из общих соображений понятно, что будет достаточно многочлена второго порядка. Попытаемся сконструировать его на основе трёх вспомогательных функций 0 (), 1 (), 2 (), удовлетворяющих условиям: 0 () равна нулю в точках 1 = 2, 2 = 4 и единице – в точке 0 , 1 ()равна нулю в точках 0 = 1, 2 = 4 и единице – в точке 1 , 2 () равна нулю в точках 0 = 1,1 = 2 и единице – в точке 2 .Тогда искомый многочлен можно выразить как сумму() = 0 ()0 + 1 ()1 + 2 ()2и он будет, очевидно, принимать в заданных точках нужные значения.(a) Базисные функции 0 (), 1 (),для заданных точек 1, 2, 42 ()(b) Интерполирующий многочлен()и его составляющиеРис.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
