§1. Источники и классификация погрешностей Часто точное численное решение математических задач бывает достаточно сложным, поэтому используют приближённые вычисления. Однако при вычислении вручную некоторые выкладки могут привести к ошибкам. Типы ошибок приближенных вычислений: вычислительные ошибки; ошибки округления; проблемы устойчивости вычислительной схемы.
Определение. Численные методы это методы решения задач, сводящиеся к арифметическим и некоторым логическим действиям над ними, т.е. к действиям которые выполняет ЭВМ. При численном решении математических и прикладных задач на том или ином этапе возможно появление погрешностей следующих типов: Погрешность задачи. Она связана с приближенным характером исходной содержательной модели (невозможно учесть все факторы в процессе изучения моделируемого явления) и ее математического описания, параметрами которого служат обычно приближенные числа.
Погрешность округлений. При выполнении арифметических операций над числами, при вводе и выводе данных производится округление. Погрешность метода. При выборе способа решения поставленной математической задачи выбирают наиболее удобный приближенный способ, который не всегда является точным. Погрешности, соответствующие этим причинам называют: Неустранимая погрешность; Устранимая погрешность; Вычислительная погрешность.
§2. Погрешность численного решения задачи Введем некоторые переменные. x* точное значение вычисляемого параметра; значение этого параметра соответствующее принятому математическому описанию; решение задачи, получаемое при реализации численного метода в предположении отсутствия округлений; приближенное решение задачи, получаемое при реальных вычислениях; неустранимая погрешность; погрешность метода; вычислительная погрешность; полная погрешность;
Чтобы численный метод считался удачно выбранным, необходимо выполнение некоторых условий: А также должно выполняться условие: Рассмотрим некоторые возможные подходы к учету погрешностей действий. Пусть А и а два близких числа. А точное значение некоторой величины; a приближенное значение этой величины. его погрешность должна быть в несколько раз меньше неустранимой погрешности ( 0 < 1 ); вычислительная погрешность в несколько раз меньше погрешности метода ( 3 < 2 ).
Определение. Число а называется приближенным значением по недостатку, если оно меньше истинного значения, и по избытку, если больше. Например, число 3,14 является приближенным значением числа по недостатку, а 2,72 – приближенным значением числа е по избытку. Величина a = |А a| называется абсолютной погрешностью приближенного числа a, его относительная погрешность; А = n n 1 0, 1 2 … любое число (общий вид); Определение. Значащими цифрами числа a называют все цифры его записи, начиная с первой ненулевой цифры слева. (27,076 пять значащих цифр, 0, три значащих цифры).
Определение. Значащую цифру называют верной, если абсолютная погрешность числа не превосходит единицы разряда, соответствующего этой цифре: Пусть a = 27,01768 a = 0,003 а) a = 27,01(7)68; «7»: единица ее разряда 0,001; 0,003 0,001 цифра неверная; б) a = 27,017(6)8; «6»: единица ее разряда 0,0001; 0,003 0,0001 цифра неверная; в) a = 27,0(1)768; «1»: единица ее разряда 0,01; 0,003 0,01 цифра верная; Теорема: Погрешность, возникающая при округлении не превосходит ½ единицы (по абсолютной величине) меньшего из оставленных разрядов.
А = 3, округляют до двух знаков после запятой: а = 3,14. a = |А а| = 0, (3, – 3,14 = 0, ). Единица меньшего из оставленных разрядов 0,01 0, ,005 теорема верна. Повторное округление не рекомендуется, т.к. приводит к увеличению погрешности. Например: Пусть число А = 34,24463, тогда а) а = 34,25; Рассмотрим абсолютные погрешности пунктов а) и b): В случае а) погрешность a = |А а| = |34,24463 – 34,25| = 0,00537, больше погрешности пункта б) a = |А а| = |34,24463 – 34,24| = 0,00463, т.е. 0, , б) а = 34,24
Следовательно, приближенное значение числа А в пункте b) является более точным. При сложении и вычитании приближенных чисел с различным числом верных значащих цифр после запятой результат округляется по минимальному числу верных цифр после запятой у исходных чисел. При умножении и делении приближенных чисел с различным числом верных значащих цифр производится округление результата с числом значащих цифр, совпадающим с минимальным числом верных значащих цифр у исходных чисел.
§ 1. Общая постановка задачи Пусть дано уравнение f(x) = 0, где f(x) непрерывная функция. Требуется вычислить действительные корни этого уравнения с заданной точностью. I.Отделение корней уравнения. Отделить корни значит указать отрезки [a i, b i ], на каждом из которых содержится ровно один корень уравнения. а) Обозначим y = f(x) и построим график этой функции, корнями уравнения являются точки пересечения графика функции с осью (ох).
y x б) Если уравнение задано в виде g(x) = h(x) (или g(x) h(x) = 0) Введем обозначения y = g(x), y = h(x) и построим эти графики в одной системе координат. Абсциссы точек пересечения и являются корнями уравнения f(x) = 0.
y x 0 y = h(x) y = g(x) II. Определение отрезка. На выбранном отрезке [a, b] находится один корень уравнения f(x) = 0.
Рассмотрим несколько методов решения уравнений с одной переменной. а) Метод половинного деления (метод вилки) Вилка это любой отрезок, на концах которого функция имеет различные знаки. Пусть функция y = f(x) определена и непрерывна при всех x на отрезке [a,b] и имеет на концах этого отрезка значения разных знаков, т.е. f(а) f(b) 0, то существует по крайней мере одна точка С из этого отрезка, значение функции в которой равна нулю. (f(с) = 0)
Будем называть отрезок [a, b] промежутком существования корня, а точку с пробной точкой. Обозначим a = a 0, b = b 0. В результате возможны два случая: 1) f(с 0 ) = 0 с 0 точное значение корня; 2) f(с 0 ) 0 ( 0) данный отрезок разбивается на два отрезка [a 0, с 0 ] и [с 0, b 0 ]. Найдём середину этого отрезка: (Соответственно, если знак функции в т.С совпадает со знаком функции в т. а 0, то вместо f(а 0 ) будем использовать f(с 0 )) Найдем значение функции в этой точке f(с 0 ). Если знак функции в т.С совпадает со знаком функции в т. b 0, то в дальнейшем вместо f(b 0 ) будем использовать f(с 0 ).
Таким образом из этих двух отрезков выбирают тот, на концах которого функция имеет значения разных знаков. Получается новый отрезок [a 1, b 1 ]. Т.е. f(а 1 ) f(b 1 ) 0, и, определив точку с 1 [a 1, b 1 ] как середину этого отрезка производим те же операции. В результате: Длина отрезка [a 1, b 1 ] равна половине длины отрезка [a 0, b 0 ]: |[a 1, b 1 ]| = ½ |[a 0, b 0 ]|; Длина отрезка [a 2, b 2 ] равна ¼ длины отрезка [a 0, b 0 ]: |[a 2, b 2 ]| = ¼ |[a 0, b 0 ]| и т.д. Вывод:
В качестве приближенного значения корня берем середину отрезка [a n, b n ] x* точное значение корня c n приближенное значение корня Метод половинного деления позволяет вычислять искомый корень с любой наперед заданной точностью. Он особенно удобен для проведения вычислений на ЭВМ. число шагов алгоритма
б) метод касательных (метод Ньютона) Пусть действительный корень уравнения f(x) = 0 изолирован на отрезке [a, b], где f(x) – непрерывная функция, а значения f(a) и f(b) на концах отрезка имеют разные знаки и обе производные f (x) и f (x) сохраняют знак на всем отрезке [a, b]. Возьмем на [a, b] такое число x 0, при котором f(x 0 ) имеет тот же знак, что и f (x 0 ), т.е. такое, что f(x 0 ) f (x 0 ) 0 (в частности за x 0 можно принять тот из концов отрезка [a, b], в котором соблюдено это условие). Проведем в точке М 0 (x 0, f(x 0 )) касательную к кривой f(x). За приближенное значение корня берется абсцисса точки пересечения этой касательной с осью (ох), т.е. точка М 1 (х 1, 0).
Рассмотрим случай, когда f (x) и f (x) имеют одинаковые знаки. Пусть у = f(x). Напишем уравнение касательной в точке М 0 (x 0, f(x 0 )) y = f(x 0 ) + f (x 0 )(x – x 0 ) – общее уравнение касательной. f (x) 0 y a b x1x1 x*x* A B x f (x) < 0 y x a b x1x1 x*x* A B
Подставим точку М 1 (x 1, 0) в уравнение касательной: 0 = f(x 0 ) + f (x 0 )(x 1 – x 0 ) f (x 0 )(x 1 – x 0 ) = – f(x 0 ) Проведем теперь касательную в точке М 1 (x 1, f(x 1 )) y = f(x 1 ) + f (x 1 )(x – x 1 ) и подставим точку М 2 (х 2, 0).
0 = f(x 1 ) + f (x 1 )(x 2 – x 1 ) f (x 1 )(x 2 – x 1 ) = – f(x 1 ) Полученная таким образом последовательность x 0, x 1, x 2, … имеет своим пределом искомый корень. Вывод: где x n x n+1 Замечание: В случае, когда f (x) и f (x) имеют разные знаки, формулы для нахождения x n+1 имеют тот же вид. (Доказать самостоятельно)
Правила выбора исходной точки x 0 : За исходную точку х 0 следует выбирать тот конец отрезка [a, b], в котором знак функции совпадает со знаком f (x). Оценка погрешности (имеет место не только для метода касательных). y = f(x) – дифференцируема на [a, b] x* – точное значение корня уравнения f(x) = 0. x* [a, b],с [a, b],с x* Рассмотрим отрезок с концами x* и с и к этому отрезку применим теорему Лагранжа: Существует такая точка из отрезка с концами x* и с для которой выполняется равенство f(с) – f(x*) = f ( )(с – x*),
так как x* – точное значение корня уравнения f(x) = 0, то f(x*) = 0 f(с) = f ( )(с – x*). Т.к. точка с не является корнем этого уравнения, то можем записать f(с) 0 f ( ) 0. Выразим (с – x*): Отсюда можем написать где – условие окончания вычислений.
в) метод хорд Метод заключается в том, что дуга графика функции f(x) = 0 на отрезке [a,b] заменяется стягивающей ее хордой. В качестве приближенного значения корня берется точка пересечения хорды с осью (ох). Пусть функция f(x) определена и непрерывна при всех x [a, b] и на [a, b] меняет знак, т.е. f(а) f(b) 0. Тогда данное уравнение имеет хотя бы один корень. I случай: f (x) и f (x) имеют одинаковые знаки.
a b B A f ( x) 0 x1x1 x2x2 x* A1A1 A2A2 x y
Напишем общее уравнение прямой, проходящей через две точки М 1 (x 1, y 1 ) и М 2 (x 2, y 2 ): (М 1 М 2 ): Используя это уравнение, проведем хорду через точки А(а, f(а)) и B(b, f(b)): возьмем точку пересечения хорды с осью (ох) с координатами у = 0, x = x 1., отсюда, или
Проведем хорду через точки А(х 1, f(х 1 )) и B(b, f(b)): возьмем точку пересечения хорды с осью (ох) с координатами у = 0, x = x 2. отсюда или Полученная таким образом последовательность x 0, x 1, x 2, … имеет своим пределом искомый корень. Вывод:
II случай: f (x) и f (x) имеют разные знаки. a b x* x2x2 x1x1 x y B2B2 B1B1 B f ( x) < 0 f ( x) 0 A
, у = 0, x = x 1, отсюда, или Проведем хорду через точки А(a, f(a)) и B(x 1, f(x 1 )):, возьмем точку пересечения хорды с осью (ох) с координатами у = 0, x = x 2., отсюда, или
Полученная таким образом последовательность x 0, x 1, x 2, … имеет своим пределом искомый корень. Вывод: В качестве начального приближения x 0 надо взять тот конец [a, b], в котором знак f(x) и знак f (x) противоположны. Оценка погрешности и условие окончания вычислений такие же, что и в методе касательных.
г) комбинированное применение методов хорд и касательных Необходимо найти действительный корень уравнения f(x) = 0, изолированный на отрезке [a, b]. Предполагается, что f(а) и f(b) имеют разные знаки, а каждая из производных сохраняет определенный знак на отрезке изоляции. Возьмем на отрезке [a, b] такую точку x 0, что f(x 0 ) и f (x 0 ) имеют одинаковые знаки. Воспользуемся формулами методов хорд и касательных:
Величины x 11 и x 12 принадлежат промежутку изоляции, причем f(x 11 ) и f(x 12 ) имеют разные знаки. Построим новую пару приближений к корню: Точки x 21 и x 22 на числовой оси расположены между точками x 11 и x 12, причем f(x 21 ) и f(x 22 ) имеют разные знаки. Вычислим теперь значения: и т.д.
Каждая из последовательностей x 11, x 21, x 31, …, x n1, … и x 12, x 22, x 32, …, x n2, … стремится к искомому корню, причем одна из последовательностей монотонно возрастает, а другая монотонно убывает. Пусть – приближенное значение корня, полученное на n-ом шаге методом касательных, а – приближенное значение корня, полученное на n-м шаге методом хорд. Тогда условием окончания вычисления будет: а
§1. Постановка задачи Рассмотрим систему из n линейных уравнений с n неизвестными (1) Предположим, что система имеет единственное решение. Найти решение системы А, т.е. с заданной точностью.
I. Точные (Прямые) методы Точные (Прямые) методы – это методы, которые приводят к решению за конечное число арифметических операций. а) Формулы Крамера: Запишем систему (1) в матричном виде: А X = В, где – матрица системы;
– вектор (столбец) неизвестных; – вектор (столбец) свободных членов. Потребуем от системы выполнение определенных условий: 1. Матрица А должна быть размерностью n n; 2. det A 0.
Тогда при небольших n можем использовать формулы Крамера:, (i = 1, 2, 3, …, n) (2) где A i – матрица n n, получена из матрицы системы А заменой i-го столбца столбцом свободных членов. При использовании формулы Крамера необходимо вычислить (n + 1) определителей. Если определители вычислять разложением по строке (столбцу), то на вычисление определителя n-го порядка будет затрачиваться n! операций умножения. Факториальный рост числа операций с увеличением размерности n называется проклятием размерности (100! ), а размерность n = 100 для современных задач не велика.
б) Метод Гаусса. Последовательное исключение неизвестных. Выпишем расширенную матрицу. (3) Предположим, что коэффициент a 11, называемый ведущим элементом первой строки, не равен нулю. Разделив первую сроку на a 11, получим где
С помощью первой строки получаем в первом столбце, начиная со второй строки, нули. Для этого, сложим первую строку со следующими строками. Умножив ее на элементы а i1 (i = 2, 3, …, n) с противоположным знаком. где Допустим, что ведущий элемент второй строки, т.е. коэффициент a 22 (1), тоже отличен от нуля. Тогда, разделив на него вторую сроку, получим
С помощью второй строки получаем во втором столбце ниже единицы все нули. Получаем и т.д. В результате получим
(4) Переход от расширенной матрицы к матрице (4) называется прямой ход метода Гаусса (получение треугольной матрицы). С помощью последней строки, получаем в последнем столбце все нули выше единицы. Затем с помощью предпоследней строки в предпоследнем столбце получаем нули выше единицы и т.д.
Получим Это преобразование называется – обратный ход метода Гаусса (переход от треугольной матрицы к единичной). (5) – решение системы (1).
в) Метод Гаусса с выбором главного элемента Рассмотренный выше простейший вариант метода Гаусса, называемый схемой единственного деления, обладает следующими недостатками. Если делить на число, близкое к нулю, то получится большая ошибка округления, поэтому если на диагонали матрицы стоят числа близкие к нулю, то приблизительное решение системы получаются с большой погрешностью. Чтобы уменьшить ошибки округления, используют метод Гаусса с выбором главного элемента. Пусть дана система (1). Выпишем расширенную матрицу. В первом столбце среди всех элементов выбираем наибольший по модулю элемент и ту строку, в которой он находится, меняем местами с первой строкой.
Затем во втором столбце среди элементов, начиная со второго, выбираем наибольший по модулю элемент и строку, в которой он находится, меняем со второй строкой и т.д. Затем как в методе Гаусса получаем в первом столбце первой строки единицу, а ниже ее нули. Обратный ход тот же, что и в методе Гаусса. Прямые методы приводят к точным решениям, если все вычисления производились без округлений. Невязкой решения называется вектор который равен |AX 0 – B|.,
1 = |a 11 x a 12 x … + a 1n x n 0 – b 1 | По малости невязки решения можно с осторожностью судить о близости найденного приближенного решения x 0 и точного решения x*. Замечание. Правило Крамера в ЭВМ не применяется, т.к. оно требует значительно большего числа арифметических действий, чем метод Гаусса. Метод Гаусса используется в ЭВМ при решении систем до порядка 10 3, а итерационные методы – до порядка n = | a n1 x a n2 x … + a nn x n 0 – b n | 2 = |a 21 x a 22 x … + a 2n x n 0 – b 2 | ……………………………………….
II. Итерационные методы Итерационные методы – это методы, в которых точное решение может быть получено лишь в результате бесконечного повторения единообразных, как правило, простых действий. а) метод простых итераций. Предположим, что диагональные элементы системы (1) не равны 0 (а ij 0). Из первого уравнения выражаем х 1, из второго – х 2 и т.д. Из n-ого – х n.
Обозначим, (6) Система (6) приведена к нормальному виду. Запишем ее в матричном виде: X = + x, где,,
х 0 = – нулевое приближение. Общая формула имеет вид: Будем решать систему (6) методом последовательных приближений. х k+1 = х k + т.е. х 0 = ; х (1) = х (0) + ; х (2) = х (1) + … Теорема 1. Если максимальная сумма модулей коэффициентов при неизвестных в правой части каждого уравнения системы (6) меньше единицы, то последовательность приближений имеет предел.
Если последовательность приближений имеет предел, то – точное решение системы (6) следовательно, системы (1). Теорема 2. Необходимым и достаточным условием сходимости метода простых итераций при любом начальном векторе x 0 к решению x* системы (6) является требование, чтобы все собственные числа матрицы были по модулю меньше 1. Условие окончания вычисления. Пусть – точное решение системы.
– k-ое приближенное значение неизвестных, вычисленных по методу итераций. Тогда | x i (k+1) – x i (k) | < для любого i = 1, 2, …, n. Оценка погрешности: Определим: (Норма) || || = max {| 1 |, | 2 |, …, | n |} Тогда max {| x 1 – x 1 (k) |, | x 2 – x 2 (k) |, …, | x n – x n (k) |}
Пример. Показать, что для данной системы итерационный процесс сходится и определить, сколько итераций следует выполнить, чтобы найти неизвестные с точностью = 10 –4. || || = max {0,42; 0,55; 0,4} 1 || || = 0,55 – итерационный процесс сходится. Найдем количество итераций, которое необходимо выполнить, чтобы найти неизвестные с точностью = 10 –4.
|| || = 0,41,, (k + 1) lg 0,55 + lg 0,41 – lg 0,45 < – 4;, (k + 1) lg 0,55 < – 4 – lg 0,41 + lg 0,45;
б) метод Зейделя. Пусть дана система линейных уравнений Ax = b, (7) Где у матрицы все диагональные элементы отличны от нуля, т.е. a ii 0, i = 1,2, …, n. Если i-ое уравнение системы (7), i = 1,2,…,n, разделить на a ii, а затем все неизвестные, кроме x i, перенести вправо, то мы придем к эквивалентной системе вида x = Cx + d, (8) где d = (d 1, d 2, …, d n ),,
Метод Зейделя состоит в том, что итерации производятся по формуле (9) где произвольны, i = 1, 2, …, n, k = 1, 2, … Итерации (9) по методу Зейделя отличаются от простых итераций тем, что при нахождении i-ой компоненты k-ого приближения сразу используются уже найденные компоненты k-ого приближения с меньшими номерами. Условия сходимости метода простых итераций и метода Зейделя не совпадают, но пересекаются. В некоторых случаях метод Зейделя дает более быструю сходимость.
Сформулируем теорему о двух различных достаточных условиях сходимости метода Зейделя. Теорема 3. Для существования единственного решения системы (7) и сходимости метода Зейделя достаточно выполнения хотя бы одного из двух условий: а), i = 1, 2, …, n; в) матрица А – симметричная положительно определенная (все ее собственные значения положительны).
§1. Интерполяционные многочлены Пусть в точках x 0, x 1, …, x n заданы значения функции f(x 0 ), f(x 1 ), …, f(x n ) (a
Часто для решения этой задачи строится алгебраический многочлен L n (x) степени n, значения которого в точках x i совпадают со значениями функции в заданных точках. (1) Точки x 0, x 1, …, x n называются узлами интерполяции. Сам многочлен L n (x i ) называется интерполяционным многочленом. Для удобства под многочленом степени n будем подразумевать многочлен не выше n.
Например, если f i = 0, i = 0, 1, …, n, то интерполяционный многочлен L n (x) 0 фактически имеет нулевую степень, но его тоже будем называть интерполяционным многочленом n-ой степени. Приближенное восстановление функции f по формуле f(x) L n (x) (2) называется интерполяцией функции f (с помощью алгебраического многочлена). Если x расположен вне минимального отрезка, содержащего все узлы интерполяции x 0, x 1, …, x n, то замену функции f по формуле (2) называют также экстраполяцией.
Терема 1. Существует единственный интерполяционный многочлен n-ой степени, удовлетворяющий условиям (1). Доказательство. Запишем выражение интерполяционного многочлена. Пусть n = 1, тогда (3) При n = 2 (4)
и, наконец, в общем случае при любом натуральном n (5) где (6) (i = 0, 1, …, n.) Действительно, выражение (3) представляет собой линейную функцию, т.е. многочлен первой степени, причем, L 1 (x 0 ) = f 0, L 1 (x 1 ) = f 1. Таким образом, требования (1) при n = 1 выполнены.
Аналогично, формула (4) задает некоторый многочлен L 2 (x) второй степени, удовлетворяющий при n = 2 условиям (1). При произвольном натуральном n функции (6), описываемая дробью, в числителе которой стоит произведение n линейных множителей, а в знаменателе – некоторое отличное от нуля число, являются алгебраическими многочленами степени n. Следовательно, функция (5) тоже является алгебраическим многочленом степени n, причем, поскольку p ni (x i )=1, а p ni (x j ) = 0 при j i, 0 j n, то выполнены требования (1).
Докажем единственность интерполяционного многочлена. Допустим, что кроме интерполяционного многочлена (5) имеется еще некоторый алгебраический многочлен n-й степени, удовлетворяющий условиям (7) Тогда согласно (1) и (7) (8) Если, то эта разность, будучи алгебраическим многочленом не выше n-ой степени, в силу основной теоремы высшей алгебры имеет не более n корней, что противоречит равенствам (8), число которых равно n + 1.
Следовательно,. Теорема 1 полностью доказана. Интерполяционный многочлен, представленный в виде (5), называется интерполяционным многочленом Лагранжа, а функции (многочлены) (6) – лагранжевыми коэффициентами.
Запишем равенство f(x) = L n (x) + R n (x), где R n (x) – остаточный член, т.е. погрешность интерполяции. Возьмем некоторую точку [a, b], обозначим n (x)= (x – x 0 ) (x – x 1 ) …(x – x n ) Тогда R n (x) = n (x) (9) Следовательно, f(x) = L n (x) + n (x) (10) Из равенства (10) вытекает оценка погрешности интерполяции (в частности, экстраполяции) в текущей точке x [a, b]:
|f(x) – L n (x)| | n (x)|, (11) где М n+1 = |f (n+1) (x)|, и оценка максимальной погрешности интерполяции на всем отрезке [a, b]: |f(x) – L n (x)| | n (x)| (12)
§4. Интерполяционный многочлен Ньютона Предположим, что узлы интерполяции отстоят друг от друга на одинаковом расстоянии x 0, x 1 = x 0 + h, x 2 = x 0 + 2h, … x k = x 0 + k h, (1) h > 0, k = 0, 1, …, n (т.е. узлы интерполяции образуют арифметическую прогрессию с разностью h) Такое расположение узлов обычно имеет место при интерполировании функций, заданных в виде таблицы с постоянным шагом.
Определение. Пусть x k = x 0 + k h, где k – целое, h > 0, f k = f(x k ). Величина f k = f k+1 – f k, называется конечной разностью первого порядка функции f в точке x k (с шагом h), Т.е. f 0 = f(x 1 ) – f(x 0 ) = y 1 – y 0, f 1 = f(x 2 ) – f(x 1 ) = y 2 – y 1, …………………………… f k = f(x k+1 ) – f(x k ) = y k+1 – y k, а величину n f k = n–1 f k+1 – n–1 f k, называют конечной разностью n-ого порядка функции f в точке x k. Т.е. 2 f k = f k+1 – f k (x k ), 3 f k = 2 f k+1 – 2 f k (x k ), и т.д.
Конечные разности функции f удобно записывать в таблице x0x1x2x3x4x0x1x2x3x4 f0f1f2f3f4f0f1f2f3f4 f 0 f 1 f 2 f 3 2 f 0 2 f 1 2 f 2 3 f 0 3 f 1
Пусть x 0, x 1 = x 0 + h, x 2 = x 0 + 2h, … x k = x 0 + k h - узлы интерполяции функции f(x). Тогда интерполяционный многочлен имеет вид P n (x) = a 0 + a 1 (x – x 0 ) + a 2 (x – x 0 )(x – x 1 ) + … + a n (x – x 0 ) …(x – x n-1 ) где a 0, a 1, …, a n найдены из условия, что P n (x i ) = f(x i ), i = 0, 1, …, n. P n (x 0 ) = a 0 = y 0 P n (x 1 ) = a 0 + a 1 (x 1 – x 0 ) = y 1 y 0 + a 1 h = y 1 ; a 1 = ; ;
P n (x 2 ) = a 0 + a 1 (x 2 – x 0 ) + a 2 (x 2 – x 0 )(x 2 – x 1 ) = y 2 y 0 + 2h + a 2 2h h = y 2 2h 2 a 2 = y 2 – y 0 – y 0 2; a 2 = ; итак, a 2 =
P n (x 3 ) = a 0 + a 1 (x 3 – x 0 ) + a 2 (x 3 – x 0 )(x 3 – x 1 ) + + a 3 (x 3 – x 0 )(x 3 – x 1 )(x 3 – x 2 ) = y 3 y 0 + 3h + 3h 2h + a 3 3h 2h h = y 3 6h 3 a 3 = y 3 – y 0 – 3 y y 0 ; a 3 = = и т.д. Общий вид a n =
Таким образом, формула Ньютона для интерполирования вперед имеет вид P n (x) = y 0 + (x – x 0 ) + (x – x 0 )(x – x 1 ) + + (x – x 0 )(x – x 1 )(x – x 2 ) + … + (x – x 0 ) …(x – x n-1 ) (2) В нем начало отсчета расположено в крайнем левом узле x 0, а используемые конечные разности идут в таблице разностей от f 0 вправо вниз. Интерполяционный многочлен (2) удобно использовать в начале таблицы и для экстраполяции левее точки x 0, т.е. < 0.
Интерполяционный многочлен с узлами x 0, x –1, …, x –n, где x – k = x 0 – k h, имеет вид P n (x) = y n + (x – x n ) + (x – x n )(x – x n-1 ) + + (x – x n )(x – x n-1 )(x – x n-2 )+ … + (x – x n ) …(x – x 1 ) (3) И называется интерполяционным многочленом Ньютона для интерполирования назад. В нем начало отсчета расположено в крайнем правом узле x 0, а используемые конечные разности идут в таблице от f 0 вправо вверх:
x -4 x -3 x -2 x -1 x 0 f -4 f -3 f -2 f -1 f 0 f -4 f -3 f -2 f -1 2 f -4 2 f -3 2 f -2 3 f -4 3 f -3 Интерполяционный многочлен (3) удобно использовать при интерполяции в конце таблицы и для экстраполяции правее точки x 0, т.е. > 0.
Рассмотрим численный процесс приближения производной f(x): (1) Выберем последовательность {h k } так, что h k 0, и вычисляем ее предел: для k = 1, 2, …, n, …(2)
Будем вычислять только конечное количество членов D 1, D 2, …, D n последовательности (2). Следовательно, для ответа следует использовать D n. Причем необходимо выбирать значение h n так, чтобы D n было хорошим приближением к производной f (x). Для примера рассмотрим функцию f (x) = e x и используем длину шагов, равную h = 1, ½ и ¼, чтобы построить секущую линию, которая проходит между точками (0; 1) и (h, f (h)) соответственно. Так как h уменьшается, то секущая приближается к касательной, как показано на рисунке.
y x 00,250,5 0,75 1 y = f(x) 1 Нужно произвести вычисления при h = 0,00001, чтобы получить приемлемый численный ответ, и для этого значения h графики касательной и секущей должны быть неразличимы.
§1. Приближенные методы вычислений определенных интегралов Формула Ньютона-Лейбница для вычисления определенного интеграла имеет вид. Однако, вычисление по этой формуле не всегда возможно. В таких случаях используются приближенные методы вычисления интегралов. Наиболее употребительными среди них являются метод прямоугольников, метод трапеций и метод парабол.
Пусть дан интеграл: - c c 0 x y y = f(x) Основная идея: Заменить подынтегральную функцию f(x) на многочлен, совпадающий с этой функцией в узлах интерполяции.
1)f(x) заменим многочленом нулевого порядка y = f(0): 2)f(x) заменим многочленом первого порядка, который совпадает с функцией f(x) в точках –с и с. Т.е. y = kx + b (площадь трапеции (а + b)/2 h) y - c c 0 x y = f(x)
3)f(x) заменим многочленом второго порядка, который совпадает с функцией f(x) в точках –с, 0 и с. Т.е. y = ax 2 + bx + c. y - c c 0 x y = f(x)
1.Метод прямоугольников Пусть дан интеграл Отрезок интегрирования [a, b] разобьем на n равных частей: a = x 0 < x 1 < x 2 < x 3
Однако для удобства вычислений поступают следующим образом: Точку x 1 выбирают таким образом, чтобы она являлась серединой первого отрезка, т.е. при разбивании отрезка на части таким образом: a = x 0 < x 2 < x 4 < x 6
a x2x2 x4x4 x6x6 b=x 2n x y y= f(x) Оценка погрешности формулы прямоугольников:, где
2.Метод трапеций Пусть дан интеграл Отрезок интегрирования [a, b] разобьем на n равных частей: a = x 0 < x 1 < x 2 < x 3
Вывод: – формула трапеций. Оценка погрешности формулы трапеций., где
3) Метод парабол (Симпсона). Пусть дан интеграл Отрезок интегрирования [a, b] разобьем на 2n равных частей: a = x 0 < x 1 < x 2 < x 3
= 2h y 0 + 2h y 0 + h 2y 0 – 2y 0 h = = h (2y 0 +2(y 1 – y 0 ) + (y 2 – 2y 1 + y 0 )) == h ( y 2 + y 1 + y 0 ) = ( y 2 + 4y 1 + y 0 ) тогда на промежутке [x 0, x 2 ] имеем:
( y 2 + 4y 1 + y 0 + y 4 + 4y 3 + y 2 + y 6 + 4y 5 + y 4 + … + 4y 2n–1 + y 2n–2 ) = = {f(a) + f(b) }
Вывод: {f(a) + f(b) } Оценка погрешности формулы парабол: где
§2. Формулы Ньютона- Котеса Необходимо вычислить Делим отрезок [a, b] на n равных частей. Шаг разбиения и x 0 = a, x i = x i –1 + h (i=1,2,…,n–1), x n = b. Тогда (1) – квадратурная формула Ньютона-Котеса, где (2) – коэффициенты Котеса. (значению x = a, соответствует значение q = 0, а x = b – значение q = n и dx = hdq)
Эти формулы определяют семейство квадратурных формул. Параметром этого семейства является число n – степень интерполяционного многочлена, которым заменяется подынтегральная функция. Рассмотрим несколько простейших частных случаев, соответствующих небольшим значениям n N. При этом конкретные формулы будем получать не на основе общих формул, а используя для этой цели вместо многочлена Лагранжа Эквивалентный ему первый интерполяционный многочлен Ньютона:
P n (x 0 + qh) = y 0 + q y y 0 + … + + n y 0 1.Пусть n=1, т.е. имеется всего две точки x 0 и x 1 =x 0 + h, в которых известны значения функции (y 0 = f(x 0 ) и y 1 = f(x 1 )) Этим точкам соответствуют значения 0 и 1 переменной q. Следовательно (4) yy yh (3)
Получена простейшая квадратурная формула трапеций, к которой можно прийти и из геометрических соображений: x0x0 hx1x1 y1y1 y0y0 y = L 1 (x) y = f(x) Остаточный член этой формулы: (5) где 1 (x 0, x 1 ) – некоторая точка.
2.Положим в (3) n = 2, т.е. проинтерполируем функцию f(x) по трем точкам: x 0, x 1 = x 0 + h, x 2 = x 0 = 2h. Тогда = h [ 2y 0 + 2(y 1 – y 0 ) + (y 2 – 2y 1 + y 0 )] = = (y 0 + 4y 1 + y 2 ) (6) Полученное приближенное равенство называется простейшей формулой Симпсона. Ее остаточный член:, (x 0, x 2 ) (7)
3.Предполагая теперь n = k, мы придем к частным формулам Ньютона-Котеса: (6) где x i = x 0 + ih, а коэффициенты B k,, и остаточные члены r k (h) задаются таблицей (точка (x 0,x k ), для каждого k своя).
Параметры некоторых частных формул Ньютона- Котеса вида (8) kBkBk a 0 (k) a 1 (k) a 2 (k) a 3 (k) a 4 (k) a 5 (k) … r r (h) …………………………
Общий вид линейной квадратурной формулы – это (8) где фиксированные аргументы x i называют узлами, а коэффициенты A i – весами (весовыми коэффициентами) квадратурной формулы (определенный интеграл приближенно равен среднему взвешенному значений подынтегральной функции, вычисленных в определенных точках промежутка интегрирования).
Все рассмотренные выше квадратурные формулы характерны тем, что узла в них брались равноотстоящими с шагом h, а веса находились в результате подмены подынтегральной функции f(x) кусочно-постоянной в случае формул прямоугольников, кусочно-линейной в случае формул трапеций, кусочно- квадратичной в случае формулы Симпсона и т.д. Например, у составной формулы трапеций набор весов получился следующий:, h, h, …, h, а у составной формулы Симпсона –
Далее откажемся от равномерного распределения узлов x i на промежутке интегрирования [a, b]. В таком случае целесообразно предварительно сделать линейную замену и преобразовать исходный интеграл к интегралу со стандартным промежутком интегрирования [–1, 1]: (9) Это равенство позволяет рассматривать вычисление интеграла
т.е. строить квадратурные формулы вида (10) от которых на основе (9) легко перейти к квадратурным формулам (8). Формула (10) имеет 2n параметров: n узлов t i и n весов A i. Если считать, что мы свободны в выборе как узлов, так и весов, можно попытаться подобрать их такими, чтобы равенство (11) было точным для многочленов степени 2n – 1 или, что тоже, для 2n степенных функций (t) = 1, t, t 2, …, t 2n – 1.
Формула (11) называется квадратурной формулой Гаусса. Ее решение упирается в решение нелинейной системы: Однако, решение этой системы затруднительно, но его не сложно обойти, если знать конечный результат. Но мы рассматривать их не будем.