Численные методы

Министерство образования и науки Российской Федерации В.Г. Пименов ЧИСЛЕННЫЕ МЕТОДЫ Учебное электронное текстовое издание Пособие предназначено для учебно-методического обеспечения унифицированных модулей бакалавриата всех направлений подготовки института математики и компьютерных наук. Подготовлено: кафедрой вычислительной математики ИМКН. Екатеринбург 2012 Пособие содержит лекции первой части двухсеместрового курса "Численные методы", который читается в институте математики и компьютерных наук кафедрой вычислительной математики для всех направлений подготовки бакалавриата. В пособие вошли разделы "Теория погрешностей", "Численные методы решения нелинейных уравнений", "Численные методы решения линейных и нелинейных систем", "Теория интерполяции", "Численное дифференцирование", "Численное интегрирование". 2 ОГЛАВЛЕНИЕ 1. Введение 2. Теория погрешностей 6 10 2.1. Числа и характеристики их точности . . . . . . . . . . . . . . . 10 2.2. Погрешность арифметических действий . . . . . . . . . . . . . . 11 2.3. Эффект вычитания близких чисел . . . . . . . . . . . . . . . . . 12 2.4. Погрешность функции нескольких переменных . . . . . . . . . . 13 2.5. Погрешность представления вещественных чисел в вычислительном устройстве . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.6. Эффект сложения чисел с разным порядком . . . . . . . . . . . 16 2.7. Ускорение сходимости ряда . . . . . . . . . . . . . . . . . . . . . 18 3. 21 Численное решение нелинейных уравнений 3.1. Постановка задачи . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.2. Описание численных методов решения нелинейных уравнений . 21 3.3. Погрешность методов . . . . . . . . . . . . . . . . . . . . . . . . 25 3.4. Метод простой итерации . . . . . . . . . . . . . . . . . . . . . . 26 3.5. Сходимость метода Ньютона. Квадратичная сходимость метода 4. Ньютона . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 Численные методы решения задач линейной алгебры 29 4.1. Численные методы решения линейных систем. Классификация методов . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 4.2. Компактная схема Гаусса . . . . . . . . . . . . . . . . . . . . . . 31 4.3. Трехдиагональная прогонка . . . . . . . . . . . . . . . . . . . . 37 3 4.4. Нормы матриц . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 4.5. Метод простой итерации для решения систем линейных уравнений. Критерий сходимости . . . . . . . . . . . . . . . . . . . . . 43 4.6. Достаточные условия метода простой итерации для решения систем линейных уравнений . . . . . . . . . . . . . . . . . . . . . . 48 4.7. Метод Якоби . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 4.8. Метод Гаусса-Зейделя . . . . . . . . . . . . . . . . . . . . . . . . 52 4.9. Неустранимая погрешность при решении линейных систем. Обу- 5. словленность матриц . . . . . . . . . . . . . . . . . . . . . . . . . 54 Численное решение систем нелинейных уравнений 58 5.1. Постановка задачи и предварительные сведения . . . . . . . . . 58 5.2. Метод Ньютона для решения систем нелинейных уравнений . . 59 5.3. Метод простой итерации для решения систем нелинейных урав- 6. нений . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 Интерполяция 63 6.1. Постановка задачи . . . . . . . . . . . . . . . . . . . . . . . . . . 63 6.2. Интерполяционный многочлен в форме Лагранжа . . . . . . . . 65 6.3. Погрешность интерполяционного многочлена Лагранжа . . . . 65 6.4. Разделенные разности и интерполяционный многочлен в форме Ньютона . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 6.5. Интерполяция с кратными узлами . . . . . . . . . . . . . . . . . 71 6.6. Разделённые разности с кратными узлами . . . . . . . . . . . . 72 6.7. Интерполяционный многочлен Эрмита . . . . . . . . . . . . . . 73 6.8. Дополнительные свойства разделенных разностей . . . . . . . . 74 7. 76 Численное дифференцирование 7.1. Общий подход к численному дифференцированию . . . . . . . . 4 76 7.2. Численное дифференцирование по двум узлам . . . . . . . . . . 76 7.3. Численное дифференцирование по трем узлам . . . . . . . . . . 77 7.4. Метод неопределенных коэффициентов . . . . . . . . . . . . . . 80 7.5. Неустранимая погрешность при численном дифференцировании 81 7.6. Выбор оптимального шага при численном дифференцировании 82 8. 83 Численное интегрирование 8.1. Интерполяционные квадратурные формулы . . . . . . . . . . . 83 8.2. Погрешность интерполяционных квадратурных формул . . . . 85 8.3. Элементарные квадратурные формулы (Формулы Ньютона-Котеса) . . . . . . . . . . . . . . . . . . . . 86 8.4. Погрешность элементарных интерполяционных квадратурных формул . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 8.5. Составные квадратурные формулы . . . . . . . . . . . . . . . . 93 8.6. Погрешность составных квадратурных формул . . . . . . . . . . 95 8.7. Метод Рунге практической оценки погрешности . . . . . . . . . 96 8.8. Формулы наивысшей алгебраической степенью точности . . . . 98 8.9. Существование и единственность квадратуры Гаусса . . . . . . 101 8.10. Алгоритм построения квадратуры Гаусса . . . . . . . . . . . . . 104 8.11. Погрешность квадратуры Гаусса . . . . . . . . . . . . . . . . . . 105 8.12. Вычисление интегралов с особенностями . . . . . . . . . . . . . 107 5 1. ВВЕДЕНИЕ Рассмотрим место изучаемой дисциплины в ряду предметов, связанных с математическим моделированием. Предположим, что имеется объект любой природы. Его требуется изучить, выявив закономерности и характеристики. С помощью конкретной науки, относящейся к объекту, составляется математическая модель, в которой отражаются самые существенные характеристики объекта. Модель содержит связи между характеристиками, которые, как правило, выражаются математическими уравнениями. С целью установить закономерности необходимо исследовать математическую модель, например решить уравнения. Но чем адекватнее составлена математическая модель, тем она сложнее, и напрямую аналитические методы решения зачастую неприменимы. Кроме аналитических методов в настоящее время все активнее применяют численные методы, алгоритмы, основанные на аппроксимации исходной математической модели другими, реализуемыми в виде вычислительных операций, зачастую громоздких. Следует отметить,что для правильной картины изучения математической модели, численные методы должны быть дополнены качественными методами, которые составляют значительную часть математики. Громоздкость численных расчетов приводит к необходимости использования вычислительных устройств (в широком смысле этого слова), и именно развитие компьютеров и их программного обеспечения привело к большинству современных достижений в науке и технике. Следуя словам академика А.А.Самарского отметим сказанное тезисом: основу математического моделирования составляет триада: модель – алгоритм – программа. В результате на выходе получаем некоторую характеристику изучаемого объекта, которую условно будем считать числом. Но эта характери6 стика в результате обработки каждым элементом триады в той или иной степени огрубляется, вносится погрешность. Задача состоит в её изучении на каждой стадии. Aнеустр. - погрешность, которая вносится на стадии математического моделирования, называется неустранимой погрешностью. Пример: Физический маятник, возможно, с трением, описывается дифференциальным уравнением второго порядка     φ̈ + a sin φ + bφ̇ = f (t),    φ(0) = α,      φ̇(0) = β; где φ - угол отклонения маятника от положения равновесия. Начальные условия α, β, а также коэффициенты уравнения a и b должны быть измерены некоторыми приборами, которые имеют определенную точность, которая и дает неустранимую погрешность. Aметода - погрешность, которая вносится на стадии составления численного алгоритма, называется погрешностью метода. Например, приведенное уравнение физического маятника аналитически не решается. Можно, конечно, заменить sin φ на φ и получим линейное уравнение второго порядка (математический маятник), решение которого выписывается в элементарных функциях, но такое упрощение модели можно сделать только при малых φ. Между тем, в теории численных методов разработаны эффективные методы решения дифференциальных уравнений, подобных уравнению физического маятника. К числу них относятся метод Эйлера, методы Рунге-Кутты, методы Адамса и другие. Подбирая методы и их параметры, прежде всего шаг дискретизации, можно получить решение с требуемой погрешностью метода. Изучение разнообразия 7 таких методов и способов оценки их погрешностей составляет основное содержание курса. Aвычисл. - погрешность, которую вносит вычислительное устройство, называется вычислительной погрешностью. Так подавляющее большинство численных алгоритмов, в том числе и упомянутые методы, требует громадных вычислений для достижения хорошей точности, поэтому необходимо эффективно работающее вычислительное устройство. Но компьютер представляет информацию дискретно, поэтому любые действия с вещественными числами, наиболее распространенными в моделировании, вносят вычислительную погрешность. По этой причине число операций нужно по возможности оптимизировать, оптимизировать нужно также и вычислительные средства. Все сказанное можно изобразить в виде схемы Aнеуст. ↓ Aметода ↓ Aвыч. ↓ обьект→ модель → метод → выч.устройство →число Формула полной погрешности Aпол. = Анеустр. + Aмет. + Aвыч. Считается, что идеальная ситуация наблюдается тогда, когда все три вида погрешности примерно равны или, по крайней мере, имеют один порядок. Например, нет смысла применять очень точные алгоритмы и считать на вычислительном устройстве с большой степенью точности, если имеется большая неустранимая погрешность. Основной задачей данной дисциплины является средняя часть триады – изучение численных алгоритмов и погрешностей метода. Однако, при 8 этом, обязательно нужно уметь учитывать влияние неустранимой погрешности, а также нужно уметь выбирать третью часть триады (программу) с соответствующим влиянием вычислительной погрешности на конечный результат. 9 2. ТЕОРИЯ ПОГРЕШНОСТЕЙ 2.1. Числа и характеристики их точности По умолчанию будем считать,что числа рассматриваются вещественные. Пусть x - идеальное вещественное число, а x∗ - его приближение. Абсолютной погрешностью числа x будет являться величина Ax∗ , удовлетворяющая соотношению |x − x∗ | 6 Ax∗ . Относительной погрешностью числа x называется величина ∆x∗ = Ax ∗ |x∗ | (|x∗ | = 6 0). Рассмотрим число в десятичном позиционном представлении. Цифра, начиная с первой ненулевой, называется значащей. Например, 20.02 - все цифры значащие, 002.02 - третья цифра является значащей. Цифра называется верной, если абсолютная погрешность числа не превосходит половины единицы соответствующего разряда. Все остальные цифры сомнительные (то есть Ax∗ 6 0.5 · 10−m ) Например, 1) число 21.145, и пусть Ax∗ = 0.01. Тогда 211 - верные цифры, 45 сомнительные. 2) число 3.1415926, Ax∗ = 0.007. Верные цифры 31, в самом деле, посмотрим на 3.14, единица разряда последней цифры 1/100, ее половина 5/1000, а погрешность больше 5/1000 < 0.007, следовательно 4 – цифра сомнительная при данной оценке погрешности. Понятие верной цифры после запятой связано с понятием абсолютной погрешности, а понятие верной значащей с относительной погрешностью. 10 2.2. Погрешность арифметических действий Cложение и вычитание. Пусть точные значения некоторых чисел – x1 , x2 , . . . , xn , а их приближения x∗1 , x∗2 , . . . , x∗n . Известны погрешности Ax∗1 , Ax∗2 , . . . , Ax∗n . Пусть треn P буется сложить числа, т.е. посчитать y = xi . В качестве приближения суммы возьмем y ∗ = n P i=1 x∗i . Оценим погрешность i=1 n n n n X X X X ∗ ∗ |y − y | = | xi − xi | 6 |xi − xi | 6 Ax∗i = Ay∗ ∗ i=1 i=1 i=1 i=1 Получили утверждение Теорема 1. При сложении и вычитании абсолютные погрешности складываются. Рассмотрим относительную погрешность сложения. n P ∆y∗ = Ax∗i n P Ay∗ = i=1 = n ∗ P |y | | x∗i | ∆x∗i |x∗i | i=1 | i=1 n P x∗i | i=1 Предположим, что складываются числа одного знака, тогда n P ∆y∗ = ∆x∗i |x∗i | i=1 n P max ∆x∗i 6 n P |x∗i | i=1 n P |x∗i | i=1 |x∗i | i=1 = max ∆x∗i Аналогичным образом доказывается, что при сложений чисел одного знака выполняется min ∆x∗i ≤ ∆y∗ . Однако при сложении чисел разных знаков (при вычитании) может наблюдаться эффект возрастания относительной погрешности. 11 2.3. Эффект вычитания близких чисел Пример: Требуется посчитать √ если √ 11 − 11 ≈ 3, 32; √ √ 10 10 ≈ 3, 16 (все цифры верные). Тогда √ 11 − √ 10 ≈ 0, 16. √ √ Оценим погрешность. Обозначим x1 = 11, x2 = 10, x∗1 = 3, 32, x∗2 = √ √ 3, 16, y = 11 − 10. Тогда Ax∗1 = Ax∗2 = 0, 005; по теореме сложения √ √ абсолютных погрешностей Ay∗ = 0, 01 и в соотношении 11 − 10 ≈ 0, 16 только один знак после запятой верный. Относительная погрешность, как можно посчитать, резко возросла. Посчитаем искомую величину по другому, используя домножение на сопряженное: √ 11 − √ 1 √ ≈ 0, 154321 11 + 10 √ √ 11 + 10 ≈ 6, 48 10 = √ Оценим погрешность, используя, кроме ранее введенных обозначений, z = √ √ 11 + 10. Тогда Az ∗ = 0, 01; ∆z ∗ = 0,01 6,48 = 0, 0016. Для оценки погрешности деления y = 1/z используем утверждение Теорема 2. При умножении и делении складываются относительные погрешности. Доказательство этого утверждения будет дано позже. 12 Тогда, так как 1 в числителе для y дана точно, ∆y∗ = ∆z ∗ = 0, 0016. Из относительной погрешности можно получить абсолютную Ay∗ = ∆y∗ |y ∗ | = 0, 0016 · 0, 16 = 0, 000256. Поэтому в числе y ∗ три верных знака после запятой. Таким образом √ 11 − √ 10 ≈ 0, 154 При этом все цифры верные. 2.4. Погрешность функции нескольких переменных Исследуем погрешность для произвольной функции многих переменных y = f (x1 , x2 , . . . , xn ) Здесь и в дальнейшем будем предполагать, что функция обладает свойствами, которые обеспечивают дальнейшие действия, например, имеет необходимую гладкость. Предположим, что заданы приближенные значения аргументов функции x∗1 , x∗2 , . . . , x∗n с абсолютными погрешностями Ax∗1 , Ax∗2 , . . . , Ax∗n . Тогда y ∗ = f (x∗1 , x∗2 , . . . , x∗n ). Оценим абсолютную погрешность функции |y − y ∗ | = |f (x1 , x2 , . . . , xn ) − f (x∗1 , x∗2 , . . . , x∗n )| Сведем разность значений функции векторного аргумента к разности значений функции скалярного аргумента. Введем векторные обозначения X = (x1 , . . . , xn ) ∈ Rn , X ∗ = (x∗1 , . . . , x∗n ) ∈ Rn 13 Введем также вспомогательную векторную функцию скалярного аргумента: X(t) = X + t(X ∗ − X). Тогда X(0) = X, X(1) = X ∗ , |f (x1 , x2 , . . . , xn ) − f (x∗1 , x∗2 , . . . , x∗n )| = |f (X(0)) − f (X(1))| = = |φ(0) − φ(1)| Здесь φ(t) = f (X(t)) – сложная функция. Применим к функции φ(t) формулу конечных приращений Лагранжа |φ(0) − φ(1)| = |φ0 (c)| · |0 − 1|, с ∈ [0; 1] Вычисляя производную функции φ(t), и производя оценку, получаем n n X X ∂f ∂f ∗ |φ(0) − φ(1)| = | (X(c))||xi − xi | 6 (X(c))|Ax∗i | ∂x ∂x i i i=1 i=1 Таким образом n X ∂f Ay∗ = | (X(c))|Ax∗i ∂x i i=1 При малых погрешностях аргументов абсолютная погрешность функции приближенно описывается формулой n X ∂f Ay ∗ ≈ | (X ∗ )|Ax∗i ∂xi i=1 (обычно пишут знак "=" вместо "≈". Рассмотрим в качестве примера вывод теоремы о погрешности умножения. Пусть y = x1 x2 14 Вычисляя частные производные этой функции двух переменных и применяя формулу для погрешностей, получаем Ay∗ = |x∗2 |Ax∗1 + |x∗1 |Ax∗2 ∆y∗ = Ax∗1 Ax∗2 + = ∆x∗1 + ∆x∗2 . |x∗1 | |x∗2 | Аналогично выводится теорема о погрешности деления. 2.5. Погрешность представления вещественных чисел в вычислительном устройстве Информация о вещественных числах в вычислительном устройстве хранится в виде знака, порядка и мантиссы. Ограниченность разрядов, отводимых для хранения мантиссы приводит к погрешности представления вещественного числа. Формы представления в разных программных средствах различные, но принцип одинаковый, поэтому проиллюстрируем на одной из распространенных, когда на хранение вещественного числа отводится 48 ячеек (бит), из которых одна хранит информацию о знаке числа, восемь ячеек хранят информацию о порядке числа, остальные 39 – о мантиссе числа. Такая форма принята, например, в стандартном представлении вещественного числа в алгоритмическом языке ПАСКАЛЬ. Пусть точное число представимо в виде x = (−1)s 2g−127 0.a1 a2 . . . a39 a40 . . . , где s – двоичная цифра, g – двоичный порядок, ai – двоичные цифры нормализованной ( a1 = 1) мантиссы. Приближенное же число имеет в силу ограниченности мантиссы форму x∗ = (−1)s 2g−127 0.a1 a2 . . . a39 15 Оценим абсолютную погрешность такого представления |x − x∗ | = 2g−127 0. . . . a40 . . . ≤ 2g−127 0. . . . 111 . . . ≤ ≤ 2g−127 2−39 = 2g−166 , таким образом Ax∗ = 2g−166 . Недостаток этой оценки состоит в том, что она зависит от порядка числа. Оценим относительную погрешность ∆x∗ = Ax∗ |x∗ | В силу нормализованности мантиссы справедлива оценка 2g−128 ≤ |x∗ | ≤ 2g−127 , откуда следует 2−40 ≤ ∆x∗ ≤ 2−39 . Воспользуемся правилом двоичной тысячи: 210 = 1024 ≈ 1000 = 103 , тогда 2−40 ≈ 10−12 , 2−39 ≈ 5 · 10−12 и получаем оценку 10−12 ≤ ∆x∗ ≤ 0.5 · 10−11 , таким образом одинадцать верных десятичных цифр после запятой. Погрешность представления вещественных чисел ведет к различным эффектам, которые нужно учитывать при проведении вычислений. 2.6. Эффект сложения чисел с разным порядком Пусть в вычислительном устройстве складываются два вещественных числа, имеющих представление x∗ = 2p · m1 и y ∗ = 2q · m2 16 Пусть порядки у них разные, скажем p >> q, при сложении происходит выравнивание порядков за счет денормализации мантиссы, поэтому, если разница в порядках очень большая, то может оказаться, что x∗ + y ∗ = x∗ Этот эффект приводит к парадоксам, которые "опровергают" известные факты классической математики. Упомянем пару примеров. Пример 1. При счете по определению гармонический ряд ∞ X 1 n=1 n сходится. В самом деле, частичная сумма ряда SN = N X 1 n=1 n считается по формуле 1 , N и, так как второе слагаемое стремится к нулю, а первое не стремится, то SN = SN −1 + начиная с некоторого номера N при реализации счета на некотором вычислительном устройстве будет наблюдаться SN = SN −1 , что означает стабилизацию частичных сумм, т.е. "сходимость" ряда. Пример 2. При счете по определению производные всех функций f (x) "равны нулю", если только x 6= 0. В самом деле, по определению f (x + h) − f (x) h→0 h f 0 (x) = lim Можно организовать цикл для счета предела, взяв, например, h = 10−1 , 10−2 , . . . При уменьшении найдется такое h, что с учетом сложения чисел с разным порядком выполняется x+h = x и числитель в определении предела равен нулю, а знаменатель хоть и мал, но не ноль. 17 2.7. Ускорение сходимости ряда В качестве иллюстрации влияния вычислительной погрешности рассмотрим задачу о суммирования ряда. Эта задача содержит много однотипных операций, которые легко реализовать с помощью любых программных средств, но в силу того, что погрешности могут накапливаться, результат может сильно искажаться. Рассмотрим пример: требуется посчитать сумму ряда S= ∞ X n=1 1 n2 + 0, 4 с точностью ε = 0, 5 · 10−6 Оценим погрешность метода, который состоит в том, что для счета ряда нужно посчитать его частичную сумму SN = N X n=1 1 : n2 + 0, 4 ∞ X |S − SN | = n=N +1 1 n2 + 0, 4 Для оценки остатка ряда используем интегральный признак и оценим интеграл ∞ X n=N +1 1 ≤ n2 + 0, 4 Z∞ dx ≤ x2 + 0, 4 N Z∞ 1 dx = x2 N N Потребуем, что погрешность метода не превосходила половины полной погрешности, т.е. |S − SN | ≤ ε 1 ≤ = 0, 25 · 10−6 , N 2 для этого нужно, чтобы было вычислено N ≥ 4 · 106 слагаемых. Так как при сложении абсолютные погрешность складываются, на подсчет каждого слагаемого (потребуем, чтобы вычислительная погрешность всей суммы 18 также была равна половине полной) отводится ε = 0, 25 · 10−12 , 2N т.е. 12 верных знаков после запятой, что невозможно сделать например, если вещественное число имеет форму представления, описанную выше. Рассмотрим метод, который позволяет ускорить сходимость ряда, резко уменьшить число слагаемых, необходимых для обеспечения заданной точности и, тем самым, уменьшить вычислительную погрешность. Пусть нужно вычислить исходный ряд ∞ X A= an n=1 Ряд B= ∞ X bn n=1 называется эталонным для исходного, если 1. Известна сумма ряда B. 2. an = λ 6= 0 n→∞ bn lim Если известен эталонный ряд, то исходный ряд можно представить в виде A= ∞ X (an − λbn ) + λB, n=1 при этом ряд C= ∞ X (an − λbn ) n=1 будет иметь более высокую скорость сходимости чем исходный ряд. В качестве эталонного ряда могут выступать ряды ∞ X 1 π2 = = 1, 6449340668 . . . 2 n 6 n=1 19 ∞ X 1 = 1, 2020569032 . . . 3 n n=1 ∞ X π4 1 = = 1, 0823232337 . . . 4 n 90 n=1 Применим этот прием к ряду, который был рассмотрен в качестве ∞ P 1 примера. Эталонный ряд для него n2 , λ = 1. n=1 Новый ряд ∞ X ∞ X 1 1 0, 4 C= ( 2 − 2) = − n + 0, 4 n n2 (n2 + 0, 4) n=1 n=1 Оценим для него погрешность метода ∞ X n=N +1 0, 4 ≤ 0, 4 n2 (n2 + 0, 4) Z∞ dx ≤ 0, 4 x2 (x2 + 0, 4) N Z∞ dx 0, 4 = x4 3N 3 N Чтобы теперь обеспечить требуемую погрешность метода 0, 4 ε ≤ = 0, 25 · 10−6 , 3 3N 2 достаточно взять порядка N = 100 слагаемых и каждое слагаемое нужно считать с точностью до 8 знаков после запятой. Замечание. При необходимости процедура ускорения сходимости ряда может быть применена еще раз, к уже полученному ряду. Так в примере, ∞ P 1 после применения эталонного ряда n4 число слагаемых во вновь полуn=1 ченном ряде сокращается до десятка. 20 3. ЧИСЛЕННОЕ РЕШЕНИЕ НЕЛИНЕЙНЫХ УРАВНЕНИЙ 3.1. Постановка задачи Рассмотрим задачу о численном решении нелинейных уравнений f (x) = 0 Несмотря на внешнюю простоту, нелинейные уравнения крайне редко решаются аналитически: даже среди алгебраических эффективно аналитически решаются только квадратные; уравнения третьей и четвертой степени хотя и имеют формулы для аналитического решения (формулы Кардано и Феррари), но в силу их громоздкости в настоящее время такие уравнения предпочитают решать численно; уравнения пятой и выше степени хотя и имеют действительный корень (в силу основной теоремы алгебры), но аналитических формул нет; трансцендентные же уравнения вообще сложны для аналитического решения. Задач решения нелинейных уравнений состоит из двух этапов: отделение корней и уточнения корней. На первом этапе находится такой отрезок [a, b], на котором существует корень (в дальнейшем будем обозначать его ξ) и он единственен. Эта работа, как правило, проводится аналитически. На втором этапе строится последовательность xn , сходящаяся к корню. По способу построения этой последовательности и различаются численные методы. Опишем простейшие из них. 3.2. Описание численных методов решения нелинейных уравнений Метод половинного деления. Пусть на отрезке [a, b] функция f (x) непрерывна и концах отрезка 21 принимает разные знаки: f (a)f (b) < 0, тогда по теореме Больцано-Коши на этом отрезке существует корень ξ. На этом утверждении и основан метод: находим середину отрезка c = a+b c и сужаем отрезок так, чтобы на его концах функция принимала разные знаки: если f (a)f (с) < 0, то в качестве нового значения правого конца отрезка нужно взять b = c, иначе a = c. Далее деление отрезка повторяется до тех пор, пока длина отрезка не станет меньше наперед заданной точности ε. Алгоритм всегда сходится, но его недостатком является большое число итераций: число итераций не зависит от функции, а только от длины отрезка. Замечание. Геометрическая интерпретация этого и других рассмотренных ниже методов обычно рассматривается на практике, поэтому рисунки с иллюстрациями методов приводятся в методической разработке по практическим и лабораторным работам дисциплины "Численные методы". Метод касательных (Ньютона). Возьмем начальное приближение x0 и проведем касательную к графику функции f (x) в этой точке. Точку пересечения касательной с осью абсцисс обозначим за x1 . Далее проводим касательную к графику функции f (x) в точке x1 и точку пересечения касательной с осью абсцисс обозначим за x2 . Далее процесс повторяется. Из уравнения касательной выводится формула xn+1 = xn − f (xn ) f 0 (xn ) В качестве начального приближения x0 берут один из концов отрезка, рекомендации по выбору начального приближения этого и других методов будут обсуждаться на практических занятиях. Как будет доказано позже, метод быстро сходится. Но его недостатком является необходимость использовать производную, а эта аналитическая работа не всегда удобна. Поэтому используют модификации метода. 22 Упрощенный (модифицированный) метод Ньютона. Возьмем начальное приближение x0 и проведем касательную к графику функции f (x) в этой точке. Точку пересечения касательной с осью абсцисс обозначим за x1 . Далее проводим через точку с координатами (x1 , f (x1 )) проводим прямую, параллельную первой касательной. Точку пересечения этой прямой с осью абсцисс обозначим за x2 . Далее процесс повторяется по параллельным прямым. Формула следующая xn+1 = xn − f (xn ) f 0 (x0 ) Этот метод требует вычисления всего одного значения производной, но проигрывает методу Ньютона в скорости сходимости. Метод хорд (неподвижных хорд). Зафиксируем начальное приближение x0 и возьмем первое приближение x1 . Через точки с координатами (x0 , f (x0 )) и (x1 , f (x1 )) на графике функции проведем хорду. Точку пересечения хорды с осью абсцисс обозначим за x2 . Далее проводим хорду через точки с координатами (x0 , f (x0 )) и (x2 , f (x2 )), пересечения этой хорды с осью абсцисс обозначим за x3 и так далее. Используя уравнение прямой линии, проходящей через две заданные точки, получаем формулу xn+1 = xn − f (xn )(xn − x0 ) f (xn ) − f (x0 ) Метод подвижных хорд. В отличие от ранее рассмотренных методов метод подвижных хорд двухшаговый – по двум меняющимся значениям считается третье, а затем происходит переадресация. Возьмем два приближения x0 и x1 . Через точки с координатами (x0 , f (x0 )) и (x1 , f (x1 )) на графике функции проведем хорду. Точку пересечения хорды с осью абсцисс обозначим за x2 . Далее проводим хорду через точки с координатами (x1 , f (x1 )) и (x2 , f (x2 )), пе23 ресечения этой хорды с осью абсцисс обозначим за x3 , далее через точки с координатами (x2 , f (x2 )) и (x3 , f (x3 )) и так далее. Используя уравнение прямой линии, проходящей через две заданные точки, получаем формулу xn+1 = xn − f (xn )(xn − xn−1 ) f (xn ) − f (xn−1 ) Из более сложных методов рассмотрим Метод секущих парабол (Майера). Пусть вторая производная функции на отрезке ограничена |f 00 (x)| ≤ M для всех x ∈ [a, b] и константа M известна. По значению функции в очередном приближении, по первой производной в этой точке и по M построим производную так, чтобы она пересекала ось абсцисс. Так для случая f (xn ) > 0 уравнение этой параболы (ветвями вниз) будет y = f (xn ) + f 0 (xn )(x − xn ) − M (x − xn )2 2 Найдем точки пересечения этой параболы с осью абсцисс, их две, договоримся брать правую, чтобы процесс дал монотонно возрастающую последовательность xn+1 = xn + f 0 (xn ) + p [f 0 (xn )]2 + 2M f (xn ) M Докажем, что между xn и xn+1 корней нет, для этого покажем, что график построенной параболы лежит ниже графика функции f (x). Разложим f (x) по формуле Тейлора в окрестности точки xn : f 00 (c) f (x) = f (xn ) + f (xn )(x − xn ) + (x − xn )2 , 2 где c некоторая точка между x и xn . Вычитая из этого уравнения уравнение параболы, получаем M + f 00 (c) y − f (x) = (x − xn )2 ≥ 0. 2 24 Таким образом метод дает последовательность, монотонно приближающуюся к корню. Как только нужная точность ε будет достигнута, можно добавить ε к очередному приближению, окажемся справа от корня и пойдем к другому корню (формулы в случае f (xn ) < 0 несколько меняются). Достоинство метода состоит в том, что он перебирает все корни на отрезке, если их число конечно, т.е. является глобальным. Замечание. Если xn близко к корню, то в формуле метода содержится вычитание близких чисел и, поэтому, нужно её преобразовать, домножив на сопряженное. 3.3. Погрешность методов Если xn – очередное приближение, ξ – корень уравнения, то требуемая точность ε достигнута, если выполняется условие |xn − ξ| ≤ ε. Однако это условие нельзя использовать для окончания итерационного процесса, т.к. корень неизвестен. Выведем другую оценку погрешности. Разложим функцию по формуле Тейлора в окрестности приближения xn и поставим в разложение корень ξ, получим 0 = f (ξ) = f (xn ) + f 0 (c)(ξ − xn ), c ∈ [ξ, xn ], откуда ξ − xn = − f (xn ) f 0 (c) Предположим, |f 0 (x)| ≥ m > 0 на отрезке [a, b], тогда |xn − ξ| ≤ f (xn ) , m и, следовательно, оценку f (xn ) ≤ ε, m можно взять в качестве условия, обеспечивающего заданную точность. 25 Замечание. Если оценка производной снизу m неизвестна, то используют эвристическое условие |xn − xn+1 | ≤ ε, которое, однако, не всегда обеспечивает заданную точность. 3.4. Метод простой итерации Чтобы исследовать главный вопрос – сходимость методов, рассмотрим еще один метод, который играет роль некоторого унифицированного, к которому сводятся все остальные. Перейдем от исходного уравнения f (x) = 0 к эквивалентному уравнению в форме x = ϕ(x) Методом простой итерации для нелинейных уравнений назовем алгоритм xn+1 = ϕ(xn ) Теорема 3. Пусть функция ϕ(x) непрерывно дифференцируема в окрестности корня ξ, причем выполняется |ϕ(ξ)| < 1. Тогда существует такая окрестность U = (ξ − r, ξ + r) корня ξ, что взяв начальное приближение x0 из этой окрестности, мы получим сходящийся к ξ метод простой итерации. Доказательство. В силу непрерывности производной выполняется kϕ(x)k ≤ q < 1 в некоторой окрестности U = (ξ − r, ξ + r) корня ξ. 26 Докажем, что если некоторое приближение xk ∈ U , то xk+1 ∈ U . В самом деле, по теореме Лагранжа |xk+1 − ξ| = |ϕ(xk ) − ϕ(ξ)| = |ϕ0 (c)||xk+1 − ξ| ≤ |xk+1 − ξ| < r. Таким образом, если взять начальное приближение x0 из этой окрестности, то все последующие приближения также будут в этой окрестности. Поэтому имеем |xk − ξ| ≤ q|xk−1 − ξ| ≤ · · · ≤ q k |x0 − ξ| откуда следует сходимость xk → ξ при k → ∞. Замечание. Можно доказать также, что если функция ϕ(x) непрерывно дифференцируема в окрестности корня ξ, причем выполняется |ϕ(ξ)| > 1, тогда существует такая окрестность U = (ξ − r, ξ + r) корня ξ, что взяв начальное приближение x0 из этой окрестности, мы получим расходящийся метод простой итерации. 3.5. Сходимость метода Ньютона. Квадратичная сходимость метода Ньютона В качестве примера доказанного утверждения о сходимости метода простой итерации, докажем сходимость метода Ньютона xn+1 = xn − f (xn ) f 0 (xn ) Метод является частным случаем метода простой итерации, в котором ϕ(x) = x − f (x) f 0 (x) Тогда, f 0 (x)f 0 (x) − f (x)f 00 (x) ϕ (x) = 1 − [f 0 (x)]2 27 Подставляя корень, в котором f (ξ) = 0, получаем ϕ0 (ξ) = 0. Таким образом, метод Ньютона сходится, если "удачно", т.е. в нужной окрестности корня взять начальное приближение. Более того, метод Ньютона обладает квадратичной скоростью сходимости, т.е. расстояние от последующего приближения до корня убывает пропорционально квадрату расстояния от последующего приближения до корня. В самом деле, ξ − xn+1 f (xn ) f (xn ) + f 0 (xn )(ξ − xn ) = ξ − xn + 0 = f (xn ) f 0 (xn ) Но используя разложение по формуле Тейлора, получаем 1 0 = f (ξ) = f (xn ) + f 0 (xn )(ξ − xn ) + f 00 (c)(ξ − xn )2 , 2 тогда 1 f 00 (c)(ξ − xn )2 2 f 0 (xn ) Если |f 0 (x)| ≥ m > 0, |f 00 (x)| ≤ M на отрезке [a, b], то ξ − xn+1 = − 1 M |ξ − xn |2 |ξ − xn+1 | ≤ 2 m Введем величину rn = 1 M |ξ − xn | , 2 m тогда rn+1 = 1 M |ξ − xn+1 | M 1 M |ξ − xn |2 ≤ = rn2 , 2 m 2m 2 m n+1 4 rn+1 ≤ rn2 ≤ rn−1 ≤ · · · ≤ r02 , что и означает квадратичную скорость сходимости. 28 4. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ ЗАДАЧ ЛИНЕЙНОЙ АЛГЕБРЫ Среди многочисленных задач линейной алгебры наибольшее значение имеет задача о решении линейных систем, к которой сводятся другие задачи. Рассмотрим эту задачу. 4.1. Численные методы решения линейных систем. Классификация методов Пусть дана линейная неоднородная система    a11 x1 + a12 x2 + · · · + a1n xn = b1 ,       a21 x1 + a22 x2 + · · · + a2n xn = b2 ,   ·································       an1 x1 + an2 x2 + · · · + ann xn = bn . Введем обозначения  a ... a1n  11  A =  . ... .  an1 ... ann      b x  1    1       , x =  . , b =  . ,      bn xn тогда система запишется в виде Ax = b Будем предполагать, что определитель матрицы A отличен от нуля, т.е. система имеет единственное решение. Все численные методы делятся на две группы: точные и итерационные. 29 В точных методах за конечное число шагов достигается точное решение системы. В итерационных методах за конечное число шагов достигается лишь некоторое приближение к точному решению. Рассмотрим сначала группу точных методов. Основной критерий при выборе точных методов – число операций, необходимых для решения системы. При этом под количеством операций обычно понимают количество операций умножения, т.к. на современных вычислительных устройствах операции сложения и вычитания выполняются настолько быстро по сравнению с операцией умножения, что ими пренебрегают в расчетах затрат времени, а операция деления встречается крайне редко, её стараются избегать. Известный из курса алгебры метод Крамера для решения линейных систем, относящийся к точным методам, в этом смысле крайне нерационален, число операций растет как факториал относительно размерности n системы. Даже на современных компьютерах для решения систем размерности 20 потребуется 108 лет. А приходится решать системы и гораздо большей размерности. Совсем другая ситуация с методом исключения неизвестных Гаусса, который также относится к точным методам. Там число операций имеет порядок n3 и для решения той же задачи потребуется менее секунды. Классическая схема исключения неизвестных обладает недостатком: многочисленные переходы от одной матрицы к другой. Она предназначена для аналитического счета, но не для программирования алгоритма. Рассмотрим вариант этого метода, который позволяет осуществить переход от исходной матрицы к последней за один шаг. 30 4.2. Компактная схема Гаусса Для того, чтобы изложение алгоритма было понятнее, распишем его для системы размерности четыре. Пусть решается система    a11 x1 + a12 x2 + a13 x3 + a14 x4 = a15 ,       a21 x1 + a22 x2 + a23 x3 + a24 x4 = a25 ,   a31 x1 + a32 x2 + a33 x3 + a34 x4 = a35 ,       a41 x1 + a42 x2 + a43 x3 + a44 x4 = a45 . Таким образом  a  11   a21 A=   a31  a41 a12 a13 a14 a22 a23 a24 a32 a33 a34 a42 a43 a44         Рассмотрим вспомогательную задачу факторизации матрицы A, т.е. представление матрицы в виде A = BC, где матрица B – левая треугольная, C – правая треугольная, по главной диагонали которой стоят единицы:    b 0 0 0 1 c12 c13 c14  11       b21 b22 0 0   0 1 c23 c24 , C =  B=     b31 b32 b33 0   0 0 1 c34    b41 b42 b43 b44 0 0 0 1         Запись этих двух матриц можно объединить в одну (подрузамевает- 31 ся, но не пишется, что на главной диагонали матрицы C стоят единицы)   b c c c  11 12 13 14      b b c c 21 22 23 24  A∗ =     b31 b32 b33 c34    b41 b42 b43 b44 Расписывая элементы первого столбца произведения матриц B и C, последовательно получаем a11 = b11 c11 = b11 , a21 = b21 , a31 = b31 , a41 = b41 , т.е. первый столбец "новой" матрицы A∗ совпадает с первым столбцом "старой" матрицы A. Расписывая элементы первой строки произведения матриц B и C, получаем a12 = b11 c12 , a13 = b11 c13 , a14 = b11 c14 , откуда c12 = a12 a13 a14 , c13 = , c14 = , b11 b11 b11 т.е. элементы первой строки "новой" матрицы A∗ получаются из соответствующих элементов "старой" матрицы A делением на диагональный элемент "новой" матрицы. Расписывая элементы второго столбца произведения матриц B и C, получаем a22 = b21 c12 + b22 , a32 = b31 c12 + b32 , a42 = b41 c12 + b42 , откуда b22 = a22 − b21 c12 , b32 = a32 − b31 c12 , b42 = a42 − b41 c12 , 32 т.е. элементы второго столбца матрицы B получаются из соответствующих элементов "старой" матрицы A, из которых нужно вычесть произведение уже подсчитанных соответствующих (стоящих в той же строке и колонке) элементов "новой" матрицы. Расписывая элементы второй строки произведения матриц B и C, получаем a23 = b21 c13 + b22 c23 , a24 = b21 c14 + b22 c24 , откуда c23 = (a23 − b21 c13 ) (a24 − b21 c14 ) , c24 = , b22 b22 т.е. элементы второй строки матрицы C получаются из соответствующих элементов "старой" матрицы A, из которых нужно вычесть произведение уже подсчитанных соответствующих (стоящих в той же строке и колонке) элементов "новой" матрицы и затем поделить на диагональный элемент "новой" матрицы. Аналогичным образом последовательно получаются оставшиеся элементы новой матрицы: b33 = a33 − b31 c13 − b32 c23 , b43 = a43 − b41 c13 − b42 c23 , c34 = (a34 − b31 c14 − b32 c24 ) , b33 b44 = a44 − b41 c14 − b42 c24 − b43 c34 Таким образом, можем сформулировать общие правила. Последовательно считаются столбец, затем строка, снова столбец, снова строка и т.д. "новой" матрицы. Элементы матрицы B, т.е.элементы "новой" матрицы, стоящие ниже и на главной диагонали, получаются из соответствующих элементов "старой" матрицы A, из которых нужно вычесть произведение уже подсчитанных соответствующих (стоящих в той же строке и колонке) элементов "новой" матрицы. Элементы матрицы C, 33 т.е.элементы "новой" матрицы, стоящие выше главной диагонали, получаются из соответствующих элементов "старой" матрицы A, из которых нужно вычесть произведение уже подсчитанных соответствующих (стоящих в той же строке и колонке) элементов "новой" матрицы и всё поделить на диагональный элемент новой матрицы. Общие формулы следующие bij = aij − j−1 X bik ckj , j = 1, . . . , n, i = j, . . . , n, k=1 aij − cij = i−1 P bik ckj k=1 bii , i = 1, . . . , n − 1, j = i + 1, . . . , n Вернемся к решению линейной системы Ax = b. Предположим, что решена задача факторизации A = BC, тогда получаем BCx = b и, обозначая Cx = y, сводим исходную задачу к задаче решения двух систем с треугольными матрицами   By = b  Cx = y Распишем для нашего случая n = 4 первую из этих систем    b11 y1 = a15       b21 y1 + b22 y2 = a25   b31 y1 + b32 y2 + b33 y3 = a35       b41 y1 + b42 y2 + b43 y3 + b44 y4 = a45 34 Отсюда    y1 =       y2 =   y3 =       y4 = a15 b11 a25 −b21 y1 b22 (a35 −b31 y1 −b32 y2 ) b33 (a45 −b41 y1 −b42 y2 −b43 y3 ) b44 т.е. формулы для yi такие же, что и для ci5 . Тогда введя обозначения ci5 = yi , мы можем унифицировать действия с основной матрицей и столбцом правых частей исходной системы. Введем "старую" расширенную матрицу Ā и "новую" расширенную матрицу Ā∗ .  a a a  11 12 13   a21 a22 a23 Ā =    a31 a32 a33  a41 a42 a43  b  11   b21 Ā∗ =    b31  b41 a14 a15 a24 a25 a34 a35 a44 a45 c12 c13 c14 c15 b22 c23 c24 c25 b32 b33 c34 c35 b42 b43 b44 b45                 Переход от элементов "старой" матрицы к элементам "новой" матрицы осуществляется по указанным выше правилам. Это прямой ход в схеме исключения неизвестных, представляющий собой тройной цикл с числом операций, пропорциональным n3 . Обратный ход алгоритма состоит в решении системы Cx = y и представляет собой двойной цикл. 35 Замечание 1. Как показывают формулы алгоритма, компактная схема Гаусса осуществима тогда и только тогда, когда все диагональные элементы bii отличны от нуля. В терминах исходной матрицы это выполняется тогда и только тогда, когда все главные миноры матрицы A отличны от нуля. Замечание 2. Если элементы bii (главные элементы) отличны от нуля, но близки к нулю по абсолютной величине, то при делении на малую величину может возникнуть большая погрешность. Чтобы избежать этого явления, применяют модификации компактной схемы Гаусса: схему с выбором главного элемента по колонке, схему с выбором главного элемента по строке, схему с выбором главного элемента по всей матрице. Замечание 3. С помощью компактной схемы можно считать определители квадратных матриц. В силу того, что определитель произведения матриц равен произведению определителей матриц, а матрицы B и C треугольные, получаем det(A) = det(BC) = det(B)det(C) = det(B) = n Y bii i=1 Компактная схема Гаусса является одним из самых эффективных алгоритмов вычисления определителей и этот алгоритм содержится во многих пакетах прикладных программ. Замечание 4. Компактной схемой Гаусса можно эффективно решать несколько систем с одной и той же основной матрицей A и с разными правыми частями. В этом случае в расширенную матрицу Ā дописываются новые столбцы правых частей, а матрица B будет одна и та же, т.е. число операций возрастет не сильно. Замечание 5. С помощью компактной схемы можно считать обратные матрицы. В самом деле, пусть X = A−1 , тогда матрица X является 36 решением матричного уравнения AX = E, E – единичная матрица. Обозначим столбцы матрицы X через xi , а через ei – столбец, на i месте у которого стоит 1, а остальные координаты равны нулю. Тогда получаем n систем линейных уравнений Axi = ei , с одной и той же основной матрицей. Эти системы эффективно, согласно замечанию 4, одновременно решаются компактной схемой Гаусса. 4.3. Трехдиагональная прогонка Кроме компактной схемы Гаусса, существует много других точных методов. Некоторые из них используют специфику системы. Так метод квадратного корня, который использует симметричность матрицы A, в два раза сокращает число операций. Но особо эффективным и потому часто применяемым является алгоритм прогонки, который использует трехдиагональность матрицы A. Матрица называется трехдиагональной, если у неё отличны от нуля только элементы, стоящие на главной диагонали и на двух соседних. Рассмотрим систему с трехдиагональной матрицей    b1 x1 + c1 x2 = d1       a2 x1 + b2 x2 + c2 x3 = d2      ·········   ai xi−1 + bi xi + ci xi+1 = di       ·········      an xn−1 + bn xn = dn 37 Если выразить из первого уравнения x1 через x2 и подставить во второе, то x2 выразится линейно через x3 и т.д. Введем линейную связь предыдущего неизвестного через последующее (обратный ход прогонки) xi−1 = λi xi + µi , λi и µi называются прогоночными коэффициентами. Чтобы их определить, подставим это выражение в i-ое уравнение ai (λi xi + µi ) + bi xi + ci xi+1 = di , откуда ci d i − ai µ i xi+1 + ai λi + bi ai λi + bi Получили связь между предыдущим и последующим неизвестным и, соxi = − гласно обозначениям, λi+1 = − ci di − ai µi , µi+1 = ai λi + bi ai λi + bi Для того, чтобы начинать счет прогоночных коэффициентов по этим формулам, выразим x1 из первого уравнения системы: c1 d1 x1 = − x2 + , b1 b1 откуда c1 d1 λ2 = − , µ2 = b1 b1 Сравнивая эти формулы с общими формулами для прогоночных коэффициентов, получаем, что цикл нужно начинать с λ1 = 0, µ1 = 0. Для того, чтобы вывести формулу, с которой нужно начинать обратный ход прогонки, подставим xn−1 в последнее уравнение системы: an (λn xn + µn ) + bn xn = dn , 38 откуда xn = dn − an µn = µn+1 , an λn + bn т.е. обратный ход нужно начинать с формулы xn = µn+1 . Алгоритм прогонки очень эффективен, т.к. число операций линейно относительно размерности системы n, этот факт дает возможность решать системы очень большой размерности. Установим условия, при которых прогонка осуществима. Будем говорить, что матрица обладает диагональным преобладанием, если для каждой строки модуль диагонального элемента больше суммы модулей всех других элементов этой строки. В нашем случае диагональное преобладание матрицы A означает, что ∆i = |bi | − |ai | − |ci | > 0 для всех i = 1, 2, . . . , n. Лемма 1. Пусть для основной матрицы системы выполнены условия диагонального преобладания. Тогда выполняется |λi | < 1 для всех i = 1, 2, . . . , n. Доказательство. Проверим утверждение леммы индукцией по i. База индукции выполняется, так как λ1 = 0. Предположим, что |λi | < 1 докажем, что |λi+1 | < 1. В силу диагонального преобладания имеем оценку |ai λi + bi | ≥ |bi | − |ai ||λi | ≥ |bi | − |ai | = ∆i + |ci | 39 По формулам прогоночных коэффициентов |λi+1 | = |ci | |ci | ≤ <1 |ai λi + bi | ∆i + |ci || Доказательство. В условиях диагонального преобладания прогонка осуществима. Это следует из того, что в процессе доказательства леммы мы показали, что модуль знаменателя в формулах для прогоночных коэффициентов оценивается снизу положительной величиной. Докажем также важное утверждение об оценке решений системы с трехдиагональной матрицей. Это утверждение часто используется для доказательства устойчивости алгоритмов, в которых используется прогонка. Лемма 2. Пусть для основной матрицы системы выполнены условия диагонального преобладания. Тогда выполняется оценка координат решения системы |di | 1≤i≤n 1≤i≤n ∆i Доказательство. Пусть j – тот номер, на котором достигается max |xi | ≤ max max |xi | = |xj |. 1≤i≤n Рассмотрим j-ое уравнение системы aj xj−1 + bj xj + cj xj+1 = dj , откуда bj xj = dj − aj xj−1 − cj xj+1 , |bj ||xj | ≤ |dj | + |aj ||xj−1 | + |cj ||xj+1 | ≤ |dj | + (|aj | + |cj |)|xj |. Перенося последнее слагаемое в левую часть неравенства, получаем ∆j |xj | ≤ |dj |, или max |xi | = |xj | ≤ 1≤i≤n 40 |dj | |di | ≤ max , 1≤i≤n ∆i ∆j что и требовалось доказать. 4.4. Нормы матриц Прежде чем рассматривать итерационные методы решения линейных систем, рассмотрим вспомогательный аппарат: понятие нормы квадратной матрицы. Это понятие вводится по аналогии к норме вектора. Вспомним это определение. Пусть x – n-мерный вектор с координатами xi . Нормой вектора называется поставленное в соответствие x число kxk, удовлетворяющее свойствам: 1. Для всякого вектора x выполняется kxk ≥ 0, причем kxk = 0 тогда и только тогда, когда вектор x нулевой. 2. Для всякого вектора x и скаляра λ выполняется kλxk = |λ|kxk. 3. Для всяких векторов x и y выполняется kx + yk ≤ kxk + kyk. Среди различных видов норм наиболее распространенными являются kxk1 = n X |xi | i=1 n X kxk2 = ( x2i )1/2 i=1 kxk∞ = max |xi | 1≤i≤n Пусть A – квадратная матрица размерности n × n с элементами aij . Нормой матрицы называется поставленное в соответствии A число kAk, удовлетворяющее свойствам: 1. Для всякой матрицы A выполняется kAk ≥ 0, причем kAk = 0 тогда и только тогда, когда матрица A нулевая. 41 2. Для всякой матрицы A и скаляра λ выполняется kλAk = |λ|kAk. 3. Для всяких матриц A и B выполняется kA + Bk ≤ kAk + kBk. Это общее определение нормы матрицы. С другой стороны, при выбранном базисе в векторном пространстве матрице соответствует линейное отображение x → Ax, поэтому на матрицу можно смотреть как на линейный оператор и ввести его норму: kAk = sup x6=0 kAxk kxk Это определение зависит от способа введения нормы вектора x и называется подчиненной нормой матрицы. Не всякая норма матрицы является подчиненной. Кроме указанных трех свойств, подчиненная норма матрицы удовлетворяет еще следующим свойствам: 4. Для всякой матрицы A и вектора x выполняется kAxk ≤ kAkkxk. 5. Для всяких матриц A и B выполняется kABk ≤ kAkkBk. Укажем без доказательства формулы, которыми описываются нормы матриц, подчиненные наиболее распространенным векторным нормам. kAk1 = max 1≤j≤n X |aij | 1≤i≤n (поколонная норма); kAk∞ = max 1≤i≤n X |aij | 1≤j≤n (построчная норма); q kAk2 = max λi (AT A) i (евклидова норма). 42 В определении евклидовой нормы AT – означает транспонированную матрицу, λi (B) – собственные числа матрицы. Определение евклидовой нормы корректно, т.к. AT A – симметричная, неотрицательно определённая матрица, поэтому её собственные числа вещественны и неотрицательны и под корнем стоит неотрицательное число. В отличие от поколонной и построчной норм, евклидова норма определяется не очень конструктивно из-за операции вычисления собственных чисел, поэтому часто используют оценку sX kAk2 ≤ (aij )1/2 i,j Кроме того, если матрица A – симметричная, то kAk2 = max |λi (A)|. i В дальнейшем будем рассматривать только подчиненные нормы. 4.5. Метод простой итерации для решения систем линейных уравнений. Критерий сходимости Преобразуем исходную систему Ax = b к виду x = Bx + C Методом простой итерации назовем алгоритм x(k+1) = Bx(k) + C, где x(k) – k-ое приближение вектора x к решению. 43 Фактически все итерационные методы можно записать в таком виде и главный вопрос в исследовании метода – условия его сходимости. Распишем несколько итераций, взяв начальное приближение x(0) . Тогда x(1) = Bx(0) + C, x(2) = Bx(1) + C = B 2 x(0) + (E + B)C, x(m) = B m x(0) + (E + B + B 2 + · · · + B m−1 )C. Вопрос сходимости упирается в условия сходимости матричного ряда в этой формуле. Справедливо следующее фундаментальное утверждение о сходимости метода простой итерации для линейных систем: Теорема 4. Метод простой итерации сходится при любом начальном приближении тогда и только тогда, когда все собственные числа матрицы B удовлетворяют условию |λi (B)| < 1 Доказательство разобьем на ряд вспомогательных утверждений. Сначала дадим определение. ε-жордановой формой матрицы называется квадратная матрица, у которой диагональные клетки являются ε-жордановыми, а остальные клетки нулевые. ε-жордановой клеткой называется клетка, у которой на главной диагонали стоит собственное число матрицы, ниже по соседней диагонали стоят числа ε, а все остальные элементы нулевые:   λ 0 ··· 0      ε λ ··· 0       ··· ··· ··· ···    0 ··· λ 44 Лемма 1. Всякая квадратная матрица B для всякого ε > 0 может быть приведена к ε-жордановой форме. Доказательство. Возьмем квадратную матрицу B и любое ε > 0, обозначим A = 1ε B. По одному из основных утверждений линейной алгебры всякую квадратную матрицу A можно привести к жордановой форме (которая отличается от ε-жордановой формы тем, что в клетках на диагонали ниже главной стоят единицы), т.е. найдется невырожденная матрица T , такая, что матрица G = T −1 AT – жорданова. Умножим это соотношение на ε, получим εG = T −1 (εA)T = T −1 BT. Из определения собственного числа Bh = λ(B)h, (h – собственный вектор) вытекает, что введенная матрица A имеет собственные числа λ(A) = 1ε λ(B), такие же собственные числа имеет жорданова матрица G. Но тогда у матрицы Gε = εG такие же собственные числа, что и у матрицы B, т.е. матрица Gε является ε-жордановой формой матрицы B. Лемма 2. Для того, чтобы B m → 0 при m → ∞, необходимо и достаточно, чтобы для всех собственных чисел матрицы B выполнялось условие |λi (B)| < 1. Доказательство. −1 m Так как G2ε = T −1 BT T −1 BT = T −1 B 2 T, то Gm B T, и обратно ε =T −1 B m = T Gm . Из свойств подчиненной нормы вытекает, что ε T −1 kGm kkB m kkT k ε k ≤ kT и −1 kB m k ≤ kT kkGm | ε kkT матрицы T и T −1 – невырожденные, поэтому сходимость к нулевой матрице матрицы B m эквивалентна сходимости к нулевой матрице матрицы Gm ε . 45 Пусть выполняется |λi (B)| < 1 для всех собственных чисел матрицы B. Тогда найдется ε > 0, такое, что max |λi (B)| + ε < 1 i Рассмотрим построчную норму матрицы Gε : kGε k∞ = max |λi (Gε )| + ε = max |λi (B)| + ε = q < 1 i i m Так как kGm ε k ≤ kGε k , то это означает сходимость к нулевой матm рице матрицы Gm ε и, следовательно матрицы B . Обратно, пусть Gm ε → 0 при m → ∞. Матрица Gm ε состоит из клеток, которые имеет следующую структуру: на главной диагонали стоят числа λi (B)m . Так как она сходится к нулевой, то λi (B)m → 0, что означает выполнимость условия |λi (B)| < 1 для всех собственных чисел матрицы B. Лемма 3. Матричный ряд E +B +B 2 +· · ·+B m +· · · сходится тогда и только тогда, когда для всех собственных чисел матрицы B выполняется условие |λi (B)| < 1. При этом ∞ X B m = (E − B)−1 m=0 Доказательство. Пусть матричный ряд сходится, тогда общий член ряда стремится к нулю, т.е. B m → 0 при m → ∞, а тогда по лемме 2 выполняется |λi (B)| < 1 для всех собственных чисел матрицы B. Обратно, пусть выполняется |λi (B)| < 1 для всех собственных чисел матрицы B. Рассмотрим тождество (E + B + B 2 + · · · + B m−1 )(E − B) = E − B m , 46 которое проверяется раскрытием скобок. Покажем, что матрица E − B обратима. В самом деле, если бы это было не так, то система линейных однородных уравнений (E − B)x = 0 имела бы нетривиальное решение x, а это означало бы, что вектор x собственный для матрицы B с собственным числом λ(B) = 1, что противоречит условию |λi (B)| < 1. Итак матрица E − B обратима и указанное тождество можно умножить справа на (E − B)−1 . В полученном соотношении E + B + B 2 + · · · + B m−1 = (E − B m )(E − B)−1 перейдем к пределу при m → ∞. По лемме 2 выполняется B m → 0 и правая часть соотношения стремится к пределу (E − B)−1 . Тогда сумма матричного ряда в левой части соотношения тоже имеет предел и ряд сходится к этому пределу: ∞ X B m = (E − B)−1 m=0 Доказательство теоремы. Вернемся к соотношению x(m) = B m x(0) + (E + B + B 2 + · · · + B m−1 )C. Пусть выполняется |λi (B)| < 1 для всех собственных чисел матрицы B, тогда по лемме 2 B m → 0 и первое слагаемое в правой части этого соотношения стремится к нулю, а по лемме 3 частичная сумма матричного ряда стремится к (E − B)−1 . Тогда правая и левая части соотношения имеют предел x = (E − B)−1 C. 47 Умножая на (E − B) и преобразуя, получаем, что вектор x удовлетворяет соотношению x = Bx + C Предположим теперь, что метод простой итерации сходится при любом начальном приближении, в частности при x(0) . Тогда последовательность x(m) = (E + B + B 2 + · · · + B m−1 )C сходится, поэтому матричный ряд сходится, и по лемме 3 выполняется |λi (B)| < 1 для всех собственных чисел матрицы B. 4.6. Достаточные условия метода простой итерации для решения систем линейных уравнений Приведенные необходимые и достаточные условия сходимости метода простой итерации не совсем удобны: они связаны со сложной задачей нахождения всех собственных значений матрицы. Следующее вспомогательное утверждение дает путь получения многих простых достаточных условий сходимости. Лемма. Модуль собственного числа матрицы не превосходит нормы матрицы: |λ(B)| ≤ kBk Доказательство. Пусть λ – собственное число матрицы B, а x – соответствующий собственный вектор, тогда по определению Bx = λx. Из свойств подчиненной нормы вытекает kBkkxk ≥ kBxk = |λx| = |λ|kxk. 48 Так как собственный вектор x ненулевой, его норма не ноль, поэтому можно в неравенстве сократить на kxk, получим требуемое. Как следствие из этого утверждения и критерия сходимости метода простой итерации для линейных систем вытекает Теорема 5. Если kBk < 1, то метод простой итерации сходится при любом начальном приближении. Доказательство получается из неравенства |λ(B)| ≤ kBk < 1 В частности, метод простой итерации сходится, если выполняется одно из следующих условий: 1). Для любого j = 1, . . . , n n X |bij | < 1 i=1 2). n X b2ij < 1 i,j=1 3). Для любого i = 1, . . . , n n X |bij | < 1 j=1 4.7. Метод Якоби Конкретные итерационные методы различаются по способу перехода от исходной линейной системы к виду, удобному для организации итерационного процесса, сходимость которого определяется матрицей B. Рассмотрим простейший способ такого перехода. 49 Пусть дана система Ax = b или в координатной форме    a11 x1 + a12 x2 + · · · + a1n xn = b1 ,       a21 x1 + a22 x2 + · · · + a2n xn = b2 ,   ·································       an1 x1 + an2 x2 + · · · + ann xn = bn . Выразим x1 из первого уравнения, x2 из второго уравнения и т.д., xn из последнего уравнения (предполагается, что диагональные элементы матрицы A отличны от нуля), получим систему вида   b1 12  x1 = − aa11 x2 − · · · − aa1n x + ,  n a  11 11     x2 = − a21 x1 − · · · − a2n xn + b2 , a22 a22 a22   ·································       xn = − an1 x1 − an2 x2 − · · · + bn ann ann ann или в векторной форме x = Bx + C, где матрица B имеет вид  − aa12 11   a21  −a 22 B=   . .  n1 n2 − aann − aann 50 ... − aa1n 11 − aa2n 22 ... . ... ...         Методом Якоби назовем итерационный процесс   (k+1) a12 (k) a1n (k)  x = − x − · · · − xn + ab111 ,  1 2 a a  11 11     b2 x(k+1) = − a21 x(k) − · · · − a2n x(k) 2 a22 1 a22 n + a22 ,   ·································       x(k+1) n1 (k) n2 (k) = − aann x1 − aann x2 − · · · + n bn ann Для того, чтобы изучить условия его сходимости, представим исходную матрицу A в виде A = L + D + R, где L – содержит элементы исходной матрицы ниже главной диагонали, все остальные элементы нулевые, D – диагональная матрица, R – содержит элементы исходной матрицы выше главной диагонали. Тогда сделанные преобразования запишутся в матричном виде следующим образом: (L + D + R)x = b, Dx = −(L + R)x + b, x = −D−1 (L + R)x + D−1 b, и матрица B в методе простой итерации определяется соотношением B = −D−1 (L + R) Теорема 6. Пусть aii 6= 0 для всех i = 1, . . . , n. Тогда для сходимости метода Якоби необходимо и достаточно, чтобы все корни уравнения det(L + λD + R) = 0 удовлетворяли соотношению |λ| < 1. Доказательство. Согласно доказанному критерию сходимости метода простой итерации, метод сходится тогда и только тогда, когда все собственные числа матрицы B по модулю меньше единицы. Выпишем характеристическое уравнение, учитывая конкретный вид матрицы в методе Якоби и проделаем 51 преобразования 0 = det(λE − B) = det(λE + D−1 (L + R)) = = det(D−1 )det(L + λD + R), откуда следует заключение теоремы. Теорема 7. Если исходная матрица A имеет диагональное преобладание, то метод Якоби сходится. Доказательство. В условиях диагонального преобладания матрицы A выпишем построчную норму матрицы B: X kBk∞ = max |bij | = max 1≤i≤n 1≤i≤n 1≤j≤n X 1≤j≤n,j6=i |aij | < 1, |aii | поэтому, в силу достаточного условия, метод сходится. 4.8. Метод Гаусса-Зейделя Посмотрим на покоординатные формулы метода Якоби. Когда счита(k+1) ется вторая координата k +1-го приближения, уже известно x1 , поэтому в формулах естественно использовать уже посчитанные приближения координат. Получаем метод, называемый методом Гаусса-Зейделя   (k+1) (k) 12 (k)  x1 = − aa11 x2 − · · · − aa1n xn + ab111 ,   11     b2 x(k+1) = − a21 x(k+1) − · · · − a2n x(k) 2 a22 1 a22 n + a22 ,   ·································       x(k+1) n1 (k+1) n2 (k+1) = − aann x1 − aann x2 − ··· + n bn ann Перепишем его в матричной форме, используя те же обозначения, что и в методе Якоби: x(k+1) = D−1 (−Lx(k+1) − Rx(k) + b). 52 Чтобы свести к методу простой итерации, приведем подобные при x(k+1) : Dx(k+1) + Lx(k+1) = −Rx(k) + b, x(k+1) = = −(L + D)−1 Rx(k) + (L + D)−1 b, откуда матрица B в методе простой итерации определяется соотношением B = −(L + D)−1 R Теорема 8. Пусть aii 6= 0 для всех i = 1, . . . , n. Тогда для сходимости метода Гаусса-Зейделя необходимо и достаточно, чтобы все корни уравнения det(λL + λD + R) = 0 удовлетворяли соотношению |λ| < 1. Доказательство. Согласно доказанному критерию сходимости метода простой итерации, метод сходится тогда и только тогда, когда все собственные числа матрицы B по модулю меньше единицы. Выпишем характеристическое уравнение, учитывая конкретный вид матрицы в методе Гаусса-Зейделя и проделаем преобразования 0 = det(λE − B) = det(λE + (L + D)−1 R) = = det((D + L)−1 )det(λ(L + D) + R), откуда следует заключение теоремы. Без доказательства приведем два достаточных условия сходимости метода Гаусса-Зейделя. Теорема 9. Если исходная матрица A имеет диагональное преобладание, то метод Гаусса-Зейделя сходится. 53 Теорема 10. Если исходная матрица A симметричная и положительно определена, то метод Гаусса-Зейделя сходится. 4.9. Неустранимая погрешность при решении линейных систем. Обусловленность матриц Рассмотрим примеры. Пример 1. Пусть решается двумерная система   x1 + x2 = 1  −x1 + x2 = 1 у которой точное решение x1 = 0, x2 = 1. Однако в коэффициентах системы и у чисел в правой части может быть (и реально практически всегда бывает) неустранимая погрешность. Пусть для простоты такая погрешность есть только в неоднородности, т.е. решается система   x1 + x2 = 1 + ∆b1  −x1 + x2 = 1 + ∆b2 , где |∆b1 | ≤ ε, |∆b2 | ≤ ε. Возникает вопрос: насколько точное решение близко к приближенному, обусловленному такой погрешностью, если величина ε мала. Геометрическая иллюстрация приведена на рисунке, который показывает, что если ε мало, то приближенное решение близко к точному. Однако это не всегда так. Пример 2. Пусть решается двумерная система   0.2x1 + x2 = 1  x2 = 1 54 у которой точное решение x1 = 0, x2 = 1 и рассмотрим систему с погрешностью в неоднородности   0.2x1 + x2 = 1 + ∆b1  x2 = 1 + ∆b2 , где |∆b1 | ≤ ε, |∆b2 | ≤ ε. Геометрическая иллюстрация показывает, что даже если ε мало, то приближенное решение может значительно отличаться от к точного. Это пример "плохой", плохо обусловленой системы. Перейдем к определениям. Рассмотрим систему с точными коэффициентами и неоднородностью Ax = b и систему с погрешностью в неоднородности A(x + ∆x) = b + ∆b Ставится вопрос об оценке абсолютной и относительной погрешности решения в зависимости от абсолютной и относительной погрешности в неоднородности. Вычитая из второй системы первую, получаем A∆x = ∆b, откуда ∆x = A−1 ∆b, k∆xk ≤ kA−1 kk∆bk С другой стороны, из неравенства kAkkxk ≥ kAxk = kbk получаем 1 kAk ≤ kxk kbk 55 Таким образом, k∆xk k∆bk ≤ kAkkA−1 k . kxk kbk Это неравенство показывает, во сколько раз может возрасти неустранимая погрешность (относительная) при решении линейных систем. Числом обусловлености матрицы назовем condA = kAkkA−1 k С учетом этого определения оценка перепишется так: k∆xk k∆bk ≤ condA . kxk kbk В случае, когда имеется погрешность не только в неоднородности, но и в коэффициентах системы, т.е. когда решается система (A + ∆A)(x + ∆x) = b + ∆b, имеет место оценка k∆bk condA( k∆Ak k∆xk kAk + kbk ) ≤ . kxk 1 − condA k∆Ak kAk Эту оценку приводим без доказательства. Отметим некоторые свойства числа обусловленности. 1). condA ≥ 1 В самом деле, так как A−1 A = E, то kA−1 kkAk ≥ kA−1 Ak = kEk = 1. Это неравенство показывает, что хорошо обусловленная матрица имеет число обусловленности близкое к единице. Так в примере 1 в поколонной норме cond1 A = 2, а в примере 2 cond1 A = 12 (посчитать), что приводит к большому влиянию неустранимой погрешности в примере 2. 56 2). cond(λA) = λcondA Проверяется по определению. 3). Если матрица A симметричная, то cond2 A = max |λ(A)| min |λ(A) Доказательство вытекает из того, что для симметричной матрицы евклидова норма равна максимальному по модулю собственному числу. 57 5. ЧИСЛЕННОЕ РЕШЕНИЕ СИСТЕМ НЕЛИНЕЙНЫХ УРАВНЕНИЙ 5.1. Постановка задачи и предварительные сведения Рассмотрим систему уравнений    f1 (x1 , x2 , · · · , xn ) = 0,       f2 (x1 , x2 , · · · , xn ) = 0,   ·····················       fn (x1 , x2 , · · · , xn ) = 0. Эту систему можно записать в векторном виде F (X) = 0, где вектор X = (x1 , x2 , · · · , xn ) ∈ Rn , F (X) – векторная функция векторного аргумента. Как и в случае одного уравнения, задача состоит из двух этапов: отделение корней - нахождение такой области в Rn , где существует и единственен корень ξ ∈ Rn , и уточнение корней – организация последовательности, сходящейся к корню X k → ξ. Теория численного решения системы уравнений во многом подобна теории численного решения одного уравнения, главное отличие состоит в определении производной векторной функции по векторному аргументу. Существуют разные определении производной векторной функции по векторному аргументу. Самыми распространенными являются производные по Фреше и по Гато. Функция F (X) называется дифференцируемой по Фреше в точке X, если найдется матрица P, такая, что для любого Y ∈ Rn выполняется kF (Y ) − F (X) − P (Y − X)k = o(k(Y − X)k), 58 где o(k(Y − X)k) – величина более высокого порядка, чем k(Y − X)k. Функция F (X) называется дифференцируемой по направлению в точке X, если для любого Z ∈ Rn найдется матрица P, такая, что выполняется F (X + tZ) − F (X) = PZ t→∞ t lim Матрица P называется производной функции F (X) в точке X. Если матрица P не зависит от направления Z, то функцию назовем дифференцируемой по Гато. Заметим, что из дифференцируемости по Фреше следует дифференцируемость по Гато. В случае дифференцируемости по Гато матрицу P можно сформировать, взяв в качестве направлений Zj = (0, 0, . . . , 0, 1, 0, . . . , 0) – j-ый единичный орт, тогда в качестве производной функции F (X) в точке X выступает матрица Якоби   ∂f1 (X) ∂f1 (X) ... ∂xn  ∂x1    F (X) =  . ... .    ∂fn (X) ∂fn (X) ... ∂xn ∂x1 5.2. Метод Ньютона для решения систем нелинейных уравнений Напомним, что для решения одномерного уравнения f (x) = 0, метод Ньютона записывался так x (k+1) (k) =x f (x(k) ) − 0 (k) f (x ) Многомерный аналог можно получить, взяв линейную часть в тейлоровском разложении 0 = F (X) = F (X (k) ) + F 0 (X (k) )(X − X (k) ) + · · · тогда X (k+1) = X (k) − (F 0 (X (k) ))−1 F (X (k) ) 59 Эти формулы содержат трудоемкую операцию взятия обратной матрицы, поэтому преобразуем их, перейдя к приращениям ∆X (k) = X (k+1) − X (k) , и умножив равенство слева на F 0 (X (k) ). Тогда метод перепишется в виде F 0 (X (k) )∆X (k) = −F (X (k) ). Это метод Ньютона-Рафсона, называемый также методом НьютонаКанторовича. Часто его также называют методом лианеризации, так как на каждом шаге он представляет собой задачу решения линейных систем относительно ∆X (k) , сводя, таким образом, задачу о решении нелинейной системы к решению последовательности линейных систем. 5.3. Метод простой итерации для решения систем нелинейных уравнений Перейдем от исходной системы к системе вида    x1 = ϕ1 (x1 , x2 , · · · , xn ),       x2 = ϕ2 (x1 , x2 , · · · , xn ),   ·····················       x3 = ϕn (x1 , x2 , · · · , xn ), или в векторном виде X = Φ(X). Методом простой итерации для нелинейных систем назовем алгоритм X (k+1) = Φ(X (k) ). 60 Как и в одномерном случае, метод простой итерации является базовым, т.к. все другие итерационные методы можно представить как метод простой итерации и тем самым исследовать вопрос о сходимости. Справедливо утверждение Теорема 11. Пусть функция Φ(X) непрерывно дифференцируема в окрестности решения ξ системы, причем выполняется kΦ0 (ξ)k < 1 для подчиненной нормы. Тогда существует такая окрестность U (ξ) вектора ξ, что начальное приближение X (0) принадлежит этой окрестности, то метод простой итерации сходится к ξ. Доказательство. Рассмотрим разность скалярной функции от векторного аргумента ϕi (Y ) − ϕi (X). Сведем разность значений функции векторного аргумента к разности значений функции скалярного аргумента. Введем вспомогательную векторную функцию скалярного аргумента: X(t) = X + t(Y − X), тогда X(0) = X, X(1) = Y, ϕi (Y ) − ϕi (X) = ϕi (X(1)) − ϕi (X(0)) = ψi (1) − ψi (0). Здесь ψi (t) = ϕi (X(t)) – сложная функция. Применим к функции ψi (t) формулу конечных приращений Лагранжа ψi (1) − ψi (0) = ψi0 (ci )(1 − 0), сi ∈ [0; 1] Вычисляя производную функции ψi (t), получаем ψi (1) − ψi (0) = n X ∂ϕi ∂xj j=1 61 (X(ci ))(yj − xj ) Возвращаясь к векторной функции векторного аргумента, получаем Φ(Y ) − Φ(X) = B(c1 , c2 , · · · , cn )(Y − X), где    B = B(c1 , c2 , · · · , cn ) =   ∂ϕ1 (X(c1 )) ∂x1 ∂ϕ1 (X(c1 )) ∂xn ... . ... . ∂ϕn (X(cn )) ∂x1 ... ∂ϕn (X(cn )) ∂xn      Пусть теперь ξ – решение системы и в силу непрерывности производной выполняется kΦ0 (X)k ≤ q < 1 в некоторой окрестности U (ξ) вектора ξ. В силу непрерывности частных производных выполняется kBk ≤ q < 1, если X и Y принадлежат этой окрестности. Тогда kΦ(Y ) − Φ(X)k ≤ qkY − Xk Сузим U (ξ) до его подмножества – некоторого шара в той метрике, в которой рассматривается подчиненная норма, с центром в ξ, обозначать будем также. Докажем, что если некоторое приближение X k ∈ U (ξ), то X k+1 ∈ U (ξ). В самом деле kX k+1 − ξk ≤ qkX k − ξk ≤ kX k − ξk Таким образом, если взять начальное приближение X 0 из этой окрестности, то все последующие приближения также будут в этой окрестности. Имеем kX k − ξk ≤ qkX k−1 − ξk ≤ · · · ≤ q k kX 0 − ξk откуда следует сходимость X k → ξ при k → ∞. 62 6. ИНТЕРПОЛЯЦИЯ Теория интерполяции в численных методах важна как базовый аппарат для решения таких задач, как численное дифференцирование, численное интегрирование, численное решение дифференциальных уравнений. 6.1. Постановка задачи Даны числа x0 , x1 , . . . , xn , называемые узлами интерполяции (пока для простоты попарно различные), а также числа f0 , f1 , . . . , fn , называемые значениями в узлах. Требуется найти функцию f (x), удовлетворяющую интерполяционным условиям f (xi ) = fi , i = 0, 1, . . . , n. Конечно, если не оговорить класс функций, то задача интерполяции может иметь неединственное решение. Будем выбирать функцию f (x) в классе линейных комбинаций некоторых функций ϕi (x), i = 0, 1, . . . , n : f (x) = C0 ϕ0 (x) + C1 ϕ1 (x) + · · · + Cn ϕn (x). Тогда задача интерполяции сводится к решению линейной системы относительно коэффициентов C0 , C1 , . . . , Cn :    С0 ϕ0 (x0 ) + C1 ϕ1 (x0 ) + · · · + Cn ϕn (x0 ) = f0       С0 ϕ0 (x1 ) + C1 ϕ1 (x1 ) + · · · + Cn ϕn (x1 ) = f1   ··········································       С0 ϕ0 (xn ) + C1 ϕ1 (xn ) + · · · + Cn ϕn (xn ) = fn Система однозначно разрешима тогда и только тогда, когда её главный определитель отличен от нуля. Это свойство называется линейной незави63 симостью на системе узлов. Отметим, что одна и та же система функций может быть на одной системе узлов зависимой, а на другой – независимой. Пример: Пусть ϕ0 ≡ 1, ϕ1 = x2 . Эта система функций независима на системе узлов x0 = 0, x1 = 1, но зависима на системе узлов x0 = −1, x1 = 1. Особую роль играют системы функций, линейно независимые на любой системе узлов. Такие системы функций называются it чебышевскими. Докажем, что система функций ϕ0 ≡ 1, ϕ1 = x, ϕ2 = x2 , . . . , ϕn = xn чебышевская. В самом деле, для любой системы попарно различных узлов x0 , x1 , . . . , xn , нужно проверить отличие от нуля определителя 1 x0 x20 ... xn0 1 x1 x21 ... xn1 . . . ... . 1 xn x2n ... xnn Но это определитель Ван-Дер-Монда, который равен Q (xi − xj ), и, так i6=j как узлы попарно различны, то определитель отличен от нуля. Система чебышевская. В дальнейшем мы будем рассматривать только эту систему, рассматривая таким образом задачу интерполяции в классе многочленов n−ой степени. В силу чебышевости этот многочлен, f (x) = Ln (x) = C0 + C1 x + C2 x2 + · · · + Cn xn , обозначаемый в честь Лагранжа, единственным образом решает задачу интерполяции. Существуют две формы его конструктивного построения, применяемые в различных ситуациях. 64 6.2. Интерполяционный многочлен в форме Лагранжа Рассмотрим сначала вспомогательную задачу интерполяции, когда все значения интерполяционного многочлена в узлах равны 0, за исключением i-го узла, где значение равно 1. Решение этой задачи дается формулой n Y li (x) = j=0,j6=i x − xj xi − xj В самом деле при всяком i = 0, 1, · · · , n, функция li (x) является многочленом n−ой степени; простой подстановкой получаем, что li (xi ) = 1, li (xj ) = 0, j 6= i, а в силу единственности других многочленов этой степени, решающих вспомогательную задачу нет. Теперь рассмотрим исходную задачу интерполяции, её решение выражается формулой Ln (x) = n X fi li (x), i=0 откуда, подставляя выражение для li (x), получаем формулу Лагранжа Ln (x) = n X i=0 n Y fi j=0,j6=i x − xj xi − xj Замечание. Интерполяционный многочлен в форме Лагранжа удобно применять, когда узлы интерполяции зафиксированы, а меняются значения функции (или функций). 6.3. Погрешность интерполяционного многочлена Лагранжа Будем предполагать, как и всюду в курсе, что выполняются требуемые для действий условия гладкости, в данном случае, что интерполируемая функция f (x) n + 1-раз непрерывно дифференцируема. 65 Обозначим погрешность интерполяции Rn (x) = f (x) − Ln (x). По постановке задачи погрешность интерполяции в узлах x0 , x1 , . . . , xn равна нулю. Но такие же корни имеет функция ωn (x) = (x − x0 )(x − x1 ) · · · (x − xn ). Представим погрешность в виде Rn (x) = K(x)ωn (x) и попробуем найти функцию K(x). Зафиксируем число x̂ 6= xi , i = 0, . . . , n и введем вспомогательную функцию F (x) = f (x) − Ln (x) − K(x̂)ωn (x). Функция F (x) имеет по крайней мере n+2 различных корня – узлы интерполяции x0 , x1 , . . . , xn , и число x̂. По теореме Ролля, её первая производная F 0 (x) имеет по крайней мере n+1 различный корень и т.д. n+1-производная имеет по крайней мере один корень, обозначим его ξ. Таким образом F (n+1) (ξ) = 0. Но F (n+1) (x) = f (n+1) (x) − K(x̂)(n + 1)!, поэтому f (n+1) (ξ) − K(x̂)(n + 1)! = 0, откуда f (n+1) (ξ) . (n + 1)! Подставим K(x̂) в формулу для погрешности интерполяции, расфиксироK(x̂) = вав x̂, получаем f (n+1) (ξ) Rn (x) = ωn (x), ξ ∈ [x0 , x1 , · · · , xn , x]. (n + 1)! 66 Здесь и в дальнейшем [x0 , x1 , · · · , xn , x] означает наименьший отрезок, содержащий узлы интерполяции x0 , x1 , · · · , xn и точку x, в которой происходит интерполяция. 6.4. Разделенные разности и интерполяционный многочлен в форме Ньютона Для того, чтобы построить интерполяционный многочлен в другой форме, более удобной в некоторых ситуациях, необходимо познакомиться с элементами теории разделенных разностей, созданной Ньютоном. Разделённые разности вводятся реккурентно. Пусть даны попарно различные узлы x0 , x1 , . . . , xn , и значения в узлах f0 , f1 , . . . , fn . Разделенными разностями нулевого порядка назовем значения в узлах f (xi ) = fi . Разделенной разностью первого порядка по узлам xi , xj , назовем величину f (xi , xj ) = f (xi ) − f (xj ) xi − xj Далее разделенные разности k-го порядка вводятся через разделенные разности k − 1-го порядка. Разделенной разностью k-го порядка по узлам x0 , x1 , . . . , xk , назовем величину f (x0 , x1 , . . . , xk ) = f (x1 , x2 , . . . , xk ) − f (x0 , x1 , . . . , xk−1 ) xk − x0 Среди свойств разделенных разностей отметим следующее представление Лемма. f (x0 , x1 , . . . , xn ) = n X k=0 f (xk ) n Q i=0,i6=k 67 (xk − xi ) Доказательство проведем индукцией по n. База индукции при n = 1: f (x0 , x1 ) = f (x0 ) − f (x1 ) f (x0 ) f (x1 ) = + x0 − x1 x0 − x1 x1 − x0 Шаг индукции. Предположим, что представление справедливо при n − 1, докажем представление при n. По определению разделенной разности n-го порядка и по индуктивному предположению получаем, приводя подобные f (x1 , x2 , . . . , xn ) − f (x0 , x1 , . . . , xn−1 ) = xn − x0 f (x0 , x1 , . . . , xn ) = = n−1 X k=0 (x0 − xn ) f (xk ) n−1 Q − n X (xk − xi ) (x0 − xn ) k=1 f (x0 ) + (x0 − xi ) + k=1 (x0 − xn ) (xk − xi ) f (xn ) + n−1 Q (xn − xi ) i=1 n−1 X = i=1,i6=k i=0,i6=k = Q n f (xk ) n Q i=0 f (xk ) n−1 Q [ (xk − xi ) 1 1 − ]= xk − x0 xk − xn i=1,i6=k = Q n f (x0 ) + (x0 − xi ) i=1 f (xn ) n−1 Q + (xn − xi ) n−1 X = k=0 n Q k=1 i=0 n X f (xk ) = (xk − xi ) i=0,i6=k f (xk ) n Q (xk − xi ) i=0,i6=k Лемма доказана. Замечание. Из леммы, в частности, следует, что разделенная разность является симметричной функцией своих аргументов, т.е. порядок аргументов можно менять. 68 Вернемся к решению задачи интерполяции. Справедливо утверждение. Теорема 12. А. Ln (x) = f (x0 ) + f (x0 , x1 )(x − x0 ) + f (x0 , x1 , x2 )(x − x0 )(x − x1 ) + · · · + +f (x0 , · · · , xn )(x − x0 ) · · · (x − xn−1 ) Б. Rn (x) = f (x, x0 , · · · , xn )(x − x0 ) · · · (x − xn ) Доказательство. Докажем сначала часть Б. Используя определение погрешности интерполяции и формулу интерполяционного многочлена в форме Лагранжа, получаем Rn (x) = f (x) − Ln (x) = f (x) − n X f (xi ) i=0 =[Q n f (x) (x − xj ) n Y j=0,j6=i x − xj = xi − xj n n X f (xi ) Y ] (x − xj ) − x − x i j=0 i=0 j=0 Сравнивая выражение в квадратных скобках с утверждением леммы (считая что x дополнительный аргумент разделенной разности), получаем Rn (x) = f (x, x0 , · · · , xn ) n Y (x − xj ). j=0 Докажем часть А. Возьмем тождество Ln (x) = L0 (x) + (L1 (x) − L0 (x)) + (L2 (x) − L1 (x)) + · · · + +(Ln (x) − Ln−1 (x)), и для любого m = 1, . . . , n, рассмотрим разности Lm (x)−Lm−1 (x). Эти разности являются многочленом степени m, имеют m корней x0 , x1 , . . . , xm−1 69 (так как Lm (xi ) = Lm−1 (xi ) = fi для i = 0, 1, . . . , m − 1), следовательно, имеет место представление Lm (x) − Lm−1 (x) = Am (x − x0 )(x − x1 ) · · · (x − xm−1 . Константу Am найдем, подставляя в разности число xm : Lm (xm ) − Lm−1 (xm ) = f (xm ) − Lm−1 (xm ) = = Am (xm − x0 )(xm − x1 ) · · · (xm − xm−1 ), с другой стороны из части Б теоремы вытекает, что f (xm ) − Lm−1 (xm ) = f (xm , x0 , · · · , xm−1 )(xm − x0 )(xm − x1 ) · · · · · · (xm − xm−1 ), откуда Am = f (xm , x0 , · · · , xm−1 ). Подставляя найденные выражения в тождество, получаем часть А теоремы. Замечание 1. Интерполяционный многочлен в форме Ньютона удобно применять, когда интерполируемая функция зафиксирована, а узлы интерполяции добавляются для обеспечения большей точности. В этом случае коэффициенты многочлена в форме Ньютона непересчитываются. Вычисления удобно оформлять в виде таблицы разделенных разностей. Замечание 2. Сравнивая погрешности интерполяционного многочлена в форме Лагранжа и в форме Ньютона, получаем связь между разделенными разностями и производными. Пусть узлы x0 , x1 , . . . , xk , попарно различны. Тогда f (k) (ξ) f (x0 , x1 , · · · , xk ) = , ξ ∈ [x0 , x1 , · · · , xk ] k! 70 6.5. Интерполяция с кратными узлами Кратный узел означает, что в этом узле задано не только значение интерполируемой функции, но и производные до какого-либо порядка. Пусть даны узел x1 кратности m1 , узел x2 кратности m2 , и т.д. узел xs кратности ms , т.е. следующие интерполяционные данные x1 f 0 (x1 ) ... f (m1 −1) (x1 ) x2 f 0 (x2 ) ... f (m2 −1) (x2 ) . . ... . xs f 0 (xs ) ... f (ms −1) (xs ) Всего интерполяционных данных m1 +m2 +· · ·+ms = n−1. Требуется построить многочлен Hn (x) (обозначаемый в честь Эрмита), степени n = m1 + m2 + · · · + ms − 1, удовлетворяющий условиям Hn(j) (xi ) = f (j) (xi ), i = 1, . . . , s, j = 0, . . . , mi − 1. Докажем, что эта задача имеет единственное решение. В самом деле, задача представляет собой неоднородную линейную систему относительно коэффициентов многочлена Hn (x). А линейная неоднородная система имеет единственное решение тогда и только тогда, когда соответствующая однородная система Hn(j) (xi ) = 0, i = 1, . . . , s, j = 0, . . . , mi − 1 имеет лишь тривиальное решение. Но последняя система означает, что многочлен n-ой степени имеет n + 1 корень (с учетом кратности). А это возможно, если только многочлен нулевой. Для того, чтобы построить искомый многочлен Эрмита Hn (x), рассмотрим вспомогательные конструкции. 71 6.6. Разделённые разности с кратными узлами Пусть даны узлы xi , xi+1 , . . . , xi+n , среди которых могут быть и совпадающие. Для того, чтобы определить разделенную разность по этим узлам, заменим набор узлов набором узлов xεi , xεi+1 , . . . , xεi+n , который во-первых, в отличие от исходного, состоит из попарно различных узлов, а во вторых xεi+k → xi+k при ε → 0. Такие наборы можно всегда организовать, например, если узлы xi , xi+1 , . . . , xi+k совпадают, то можно при малом ε взять xεi+j = xi + jε при j = 0, 1, · · · , k. Разделенной разностью по возможно совпадающим узлам назовем f (xi , xi+1 , . . . , xi+n ) = lim f (xεi , xεi+1 , . . . , xεi+n ), ε→0 если такой предел существует. Отметим два частных случая для конструктивного вычисления такой разделенной разности. А. Пусть среди узлов есть хотя бы одна пара различных. Пусть это первый узел xi и последний xi+n . Тогда при фиксированном ε можно воспользоваться реккурентным определением разделенной разности (узлы различные) и сохранить при переходе к пределу реккурентное определение. f (xi , xi+1 , . . . , xi+n ) = lim f (xεi , xεi+1 , . . . , xεi+n ) = ε→0 f (xεi+1 , . . . , xεi+n ) − f (xεi , . . . , xεi+n−1 ) = lim = ε→0 xεi+n − xεi = f (xi+1 , . . . , xi+n ) − f (xi , . . . , xi+n−1 ) xi+n − xi Б. Пусть все узлы совпадают xi = xi+1 = . . . = xi+n . Тогда при фиксированном ε можно воспользоваться связью разделенных разностей и производных и сохранить эту связь при переходе к пределу. f (xi , xi+1 , . . . , xi+n ) = f (xi , xi , . . . , xi ) = 72 = lim f (xεi , xεi+1 , . . . , xεi+n ) ε→0 f (k) (ξ ε ) f (k) (xi ) = = lim ε→0 n! n! 6.7. Интерполяционный многочлен Эрмита Рассмотрим исходную задачу интерполяции с кратными узлами и возьмем узлы с учетом кратности, введя двойную нумерацию узлов: x1,1 = x1,2 = · · · = x1,m1 , x2,1 = x2,2 = · · · = x2,m2 , · · · , xs,1 = xs,2 = · · · = xs,ms Заменим этот набор узлов набором попарно различных узлов xε1,1 = xε1,2 = · · · = xε1,m1 , xε2,1 = xε2,2 = · · · = xε2,m2 , · · · , xεs,1 = xεs,2 = · · · = xεs,ms и построим по ним интерполяционный многочлен степени n в форме Ньютона. Переходя к пределу при ε → 0, получаем Hn (x) = f (x1 ) + f 0 (x1 )(x − x1 ) + · · · + f (m1 −1) (x1 ) (x − x1 )m1 −1 + (m1 − 1)! +f (x1 , x1 , · · · , x1 , x2 )(x − x1 )m1 + · · · + +f (x1 , x1 , · · · , x1 , · · · , xs )(x − x1 )m1 · · · (x − xs )ms −1 При этом погрешность выражается формулой Rn (x) = f (x, x1 , x1 , · · · , x1 , · · · , xs )(x − x1 )m1 · · · (x − xs )ms Заметим, что формула Эрмита с одной стороны обобщает формулу Ньютона с её аппаратом разделенных разностей, а с другой стороны формулу Тейлора с её аппаратом производных. 73 6.8. Дополнительные свойства разделенных разностей Установим связь между разделенными разностями и производными в общем случае. Выше такая связь была установлена в двух частных случаях: когда узлы попарно различны и когда все узлы совпадают. По определению f (xi , xi+1 , . . . , xi+n ) = lim f (xεi , xεi+1 , . . . , xεi+n ). ε→0 Так как узлы xεi , xεi+1 , . . . , xεi+n попарно различные, то для разделенной разности по этим узлам можно воспользоваться связью с производной и получаем f (n) (ξ ε ) ε , ξ ∈ [xεi , xεi+1 , . . . , xεi+n ]. ε→0 n! f (xi , xi+1 , . . . , xi+n ) = lim Последовательность ξ ε ограничена, из неё можно выделить сходящуюся подпоследовательность ξ ε → ξ, ξ ∈ [xi , xi+1 , . . . , xi+n ]. Таким образом f (n) (ξ) , ξ ∈ [xi , xi+1 , . . . , xi+n ] f (xi , xi+1 , . . . , xi+n ) = n! в случае произвольных узлов. Для дальнейшего нам потребуется уметь дифференцировать разделенную разность по одному из узлов. Справедлива формула d f (x, xi , xi+1 , . . . , xi+n ) = f (x, x, xi , xi+1 , . . . , xi+n ). dx В самом деле, по определению производной и индуктивному определению разделенной разности d f (x, xi , xi+1 , . . . , xi+n ) = dx f (x + ∆x, xi , xi+1 , . . . , xi+n ) − f (x, xi , xi+1 , . . . , xi+n ) = ∆x→0 ∆x = lim = lim f (x + ∆x, x, xi , xi+1 , . . . , xi+n ) = f (x, x, xi , xi+1 , . . . , xi+n ). ∆x→0 74 Аналогичным образом проверяется формула d2 f (x, xi , xi+1 , . . . , xi+n ) = 2f (x, x, x, xi , xi+1 , . . . , xi+n ). dx2 75 7. ЧИСЛЕННОЕ ДИФФЕРЕНЦИРОВАНИЕ Даны узлы x0 , x1 , . . . , xn , и значения функции в узлах f0 , f1 , . . . , fn , требуется найти значения производных функции в некоторых точках. 7.1. Общий подход к численному дифференцированию Общий подход состоит в том, что строится интерполяционный многочлен Ln (x) и он дифференцируется, при этом погрешности численного дифференцирования Rд равна производной погрешности интерполяции: f 0 (x) = L0n (x) + Rn0 (x), f 0 (x) ≈ L0n (x), Rд = Rn0 (x). 7.2. Численное дифференцирование по двум узлам Заданы узлы x0 , x1 , (в дальнейшем будем считать, что h = x1 − x0 > 0) и значения f0 , f1 в узлах. Выпишем интерполяционный многочлен, например, в форме Ньютона, и его продифференцируем. L1 (x) = f (x0 ) + f (x0 , x1 )(x − x0 ), L01 (x) = f (x0 , x1 ) = f1 − f0 , h тогда формула численного дифференцирования по двум узлам будет f 0 (x) ≈ f1 − f0 h Погрешность интерполяции R1 (x) = f (x, x0 , x1 )(x − x0 )(x − x1 ), дифференцируя погрешность, в том числе используя правило дифференцирования разделенной разности по узлу, получаем Rд = R10 (x) = f (x, x, x0 , x1 )(x − x0 )(x − x1 )+ 76 +f (x, x0 , x1 )(2x − x0 − x1 ). Погрешность дифференцирования зависит от x, в отличие от самой формулы. Рассмотрим наиболее важные частные случаи. А. Дифференцирование на левый край x = x0 . Тогда, используя связь разделенной разности и производной, получаем f 00 (ξ) Rд = −f (x0 , x0 , x1 )h = − h, ξ ∈ [x0 , x1 ] 2 Аналогичным образом, при дифференцировании на правый край x = x1 , получаем f 00 (ξ) Rд = −f (x0 , x0 , x1 )h = h, ξ ∈ [x0 , x1 ] 2 Обычно при дифференцировании на левый или правый край применяют оценку погрешности |Rд | ≤ M2 h, M2 = max |f 00 (x)| 2 x∈[x0 ,x1 ] Б. Дифференцирование на середину x = x0 +x1 2 . Используя связь раз- деленной разности и производной, получаем Rд = −f ( x0 + x1 x0 + x1 h2 h2 , , x0 , x1 ) = −f 000 (ξ) , ξ ∈ [x0 , x1 ], 2 2 4 24 и оценку |Rд | ≤ M3 2 h , M3 = max |f 000 (x)| 24 x∈[x0 ,x1 ] Сравнивая погрешности численного дифференцирования на середину и край, заметим, что порядок малости по h на середину второй, а на край только первый. 7.3. Численное дифференцирование по трем узлам Пусть даны три равноотстоящих узла x0 , x1 , x2 ( x1 − x0 = x2 − x1 = h > 0) и значения f0 , f1 , f2 в узлах. Требуется построить формулы чис77 ленного дифференцирования для первой и второй производной и оценить погрешности. Выпишем интерполяционный многочлен, например, в форме Ньютона, и его продифференцируем два раза L2 (x) = f (x0 ) + f (x0 , x1 )(x − x0 ) + f (x0 , x1 , x2 )(x − x0 )(x − x1 ), L02 (x) = f (x0 , x1 ) + f (x0 , x1 , x2 )(2x − x0 − x1 ), L002 (x) = 2f (x0 , x1 , x2 ). Распишем разделенную разность второго порядка f (x1 , x2 ) − f (x0 , x1 ) f (x0 , x1 , x2 ) = = 2h = f2 −f1 h − f1 −f h = 2h f0 − 2f1 + f2 2h2 тогда формулы численного дифференцирования по трем узлам будут f 0 (x) ≈ f1 − f0 f0 − 2f1 + f2 (2x − x0 − x1 ), + h 2h2 f 00 (x) ≈ f0 − 2f1 + f2 h2 Подставляя различные x в формулу для первой производной мы можем получить несколько вариантов, например, формула численного дифференцирования на левый край (x = x0 ) по трем равноотстоящим узлам выглядит так: f 0 (x0 ) ≈ −3f0 + 4f1 − f2 2h Но наибольшее значение имеет формула численного дифференцирования для второй производной по трем равноотстоящим узлам, которая не зависит от x, это вообще одна из самых распространенных формул в современной математике. 78 Изучим вопрос о погрешности, особенно формулы для второй производной. Погрешность интерполяции R2 (x) = f (x, x0 , x1 , x2 )(x − x0 )(x − x1 )(x − x2 ) и вычисление производных становится очень громоздким. Чтобы немного сократить записи, введем обозначения αi = x − xi , i = 0, 1, 2, тогда R2 (x) = f (x, x0 , x1 , x2 )α0 α1 α2 R20 (x) = f (x, x, x0 , x1 , x2 )α0 α1 α2 + +f (x, x0 , x1 , x2 )[α1 α2 + α0 α2 + α0 α1 ] Отсюда можно получать оценки погрешности при разных x для формул численного дифференцирования для первой производной. Но мы это делать не будем, а найдем вторую производную R200 (x) = 2f (x, x, x, x0 , x1 , x2 )α0 α1 α2 + +2f (x, x, x0 , x1 , x2 )[α1 α2 + α0 α2 + α0 α1 ]+ +2f (x, x0 , x1 , x2 )(α0 + α1 + α2 ) Рассмотрим наиболее важный случай – формулу дифференцирования на середину (x = x1 ) для второй производной. Тогда α1 = 0, α0 = h, α2 = −h, и в формуле для погрешности остается только одно слагаемое Rд = R200 (x1 ) = 2f (x1 , x1 , x0 , x1 , x2 )(−h2 ) Выражая разделенную разность четвертого порядка через соответствующую производную, получаем погрешность Rд == −f 0000 (ξ) h2 , ξ ∈ [x0 , x2 ] 12 и оценку погрешности |Rд | ≤ M4 2 h , M4 = max |f 0000 (x)| 12 x∈[x0 ,x2 ] 79 7.4. Метод неопределенных коэффициентов Приведем простой способ вывода формул численного дифференцирования, этот способ можно успешно применять и в других разделах численных методов. Так как интерполяционный многочлен линеен относительно значений интерполируемой функции, а операция дифференцирования также линейна, то и формула численного дифференцирования линейна относительно заданных значений функции. Это дает возможность искать формулу в виде линейной комбинации заданных значений функции. Коэффициенты линейной комбинации можно получить, подставляя в формулу многочлены, начиная с низших степеней, так как интерполяционный многочлен совпадает с интерполируемой функцией, которая сама является многочленом не более высокой степени, а дифференцирование сохраняет равенство. Проиллюстрируем метод примером. Пусть в трех равноотстоящих узлах x0 , x1 , x2 ( x1 −x0 = x2 −x1 = h > 0) заданы значения f0 , f1 , f2 функции f (x). Требуется построить формулу для f 0 (x0 ). Запишем искомую формулу в виде f 0 (x0 ) ≈ Af0 + Bf1 + Cf2 . Приближенное равенство становится точным, если f (x) является многочленом не выше второй степени. Последовательно подставим в равенство f (x) ≡ 1, f (x) = x − x0 , f (x) = (x − x0 )2 , получим систему уравнений     0=A+B+C    1 = Bh + 2CH      0 = BH 2 + 4Ch2 , Решая эту систему получаем C = −1/(2h), B = 2/h, A = −3/(2h). 80 Подставляя коэффициенты, получаем формулу f 0 (x0 ) ≈ −3f0 + 4f1 − f2 , 2h уже выведенную нами ранее путем построения интерполяционного многочлена и его последующего дифференцирования. Отметим, что метод неопределенных коэффициентов очень удобен для вывода формул, но не дает оценку погрешностей. 7.5. Неустранимая погрешность при численном дифференцировании На примере формулы численного дифференцирования по двум узлам на левый край рассмотрим погрешность с учетом неустранимой. Пусть вместо точных значений функции f0 , f1 заданы их приближенные значения f0∗ , f1∗ с погрешностью Af . Оценим погрешность формулы f1∗ − f0∗ f (x0 ) ≈ h С учетом погрешности метода и влияния неустранимой погрешности имеем f1 − f0 f1 − f0 f1∗ − f0∗ f1∗ − f0∗ |f (x0 ) − | ≤ |f (x0 ) − |+| − |≤ h h h h ≤ M2 2Af h+ 2 h Таким образом, полная погрешность этой формулы численного дифференцирования оценивается (без учета влияния вычислительной погрешности) величиной Aпол = M2 2Af h+ 2 h Эта оценка показывает, что при уменьшении величины шага погрешность сначала уменьшается, а затем увеличивается. Причина этого явления состоит в том, что операция численного дифференцирования относится к 81 числу некорректных задач, т.е. задач, в которых малые изменения в исходных данных могут дать большие изменения в результате. Это явление характерно не только для данной формулы, но и для других формул численного дифференцирования. 7.6. Выбор оптимального шага при численном дифференцировании Как показывает формула полной погрешности численного дифференцирования по двум узлам на левый край, существует оптимальный шаг hопт при котором полная погрешность минимальна. Решая средствами математического анализа задачу на экстремум, получаем r M 2A Af 2 f A0пол (hопт ) = − 2 = 0, hопт = 2 , 2 h M2 p Aпол (hопт ) = 2 Af M2 Рассмотрим пример. Пусть даны четырехзначные таблицы функции y = sin x (четырехзначные означает, что все четыре знака в таблице верные и что неустранимая погрешность Af = 0, 5 · 10−4 ). Требуется вычислить в одном из узлов функцию y = cos x по формуле численного дифференцирования на левый край и оценить количество верных знаков в результате. Каким следует взять оптимальный шаг h = hопт в формуле численного дифференцирования? Табличный шаг аргумента x предполагается равным 0,001. В задаче, прежде всего, нужно определить при данной информации hопт . По выведенной выше формуле, взяв в качестве M2 = 1, получаем hопт ≈ 0, 014, т.е. 14 табличных шагов. При этом Aпол (hопт ) = 0, 014, т.е. всего один верный знак после запятой. И этот результат при данной информации улучшить за счет выбора шага нельзя. 82 8. ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ В отличие от численного дифференцирования, численное интегрирование – корректная задача. Эта одна из самых старых практических задач математики, её история отразилась в терминологии (квадратура). 8.1. Интерполяционные квадратурные формулы Основная задача – посчитать приближенно величину определенного интеграла, используя вычисление значений подинтегральной функции. Формула вида Zb f (x)dx ≈ n X Ak f (xk ) k=0 a называется квадратурной. Один из основных подходов состоит в замене подинтегральной функции интерполяционным многочленом и его последующем интегрировании. Пусть f (x) ≈ Ln (x), запишем интерполяционный многочлен в форме Лагранжа Ln (x) = n X n Y f (xk ) k=0 j=0,j6=k x − xj , xk − xj подставим его в интеграл и поменяем местами интеграл и сумму. Получим Zb f (x)dx ≈ a n X f (xk ) k=0 Обозначим Ak = Zb Y n a Zb Y n a j=0,j6=k j=0,j6=k x − xj dx xk − xj x − xj dx xk − xj и получим квадратурную формулу, которая называется интерполяционной квадратурной. 83 Интерполяционные квадратуры обладают свойством, которое выделяет их из всех квадратур. Теорема 13. Для того, чтобы квадратурная формула была интерполяционной, необходимо и достаточно, чтобы она была точна для любого многочлена степени n и ниже. Необходимость. Пусть квадратура интерполяционная, т.е. получилась заменой подинтегральной функции её интерполяционным многочленом Ln (x). Но если подинтегральная функция является многочленом степени n и ниже, то она совпадает с интерполяционным многочленом и квадратура становится точной. Достаточность. Пусть квадратура точна для любого многочлена степени n и ниже. Возьмем в качестве такого многочлена Qi (x) = n Y j=0,j6=i x − xj xi − xj и получим формулу Zb Qi (x)dx = n X Ak Qi (xk ) k=0 a Но Qi (xk ) = 0, k 6= i, Qi (xi ) = 0, и последняя сумма вырождается в одно слагаемое Ai . Подставляя в эту формулу Qi (x), получаем Zb Y n a j=0,j6=i x − xj dx = Ai , xi − xj что совпадает с определением интерполяционной квадратуры. Замечание. Введем важное определение алгебраической степени точности формулы. Число N называется алгебраической степенью точности формулы, 84 если 1) она точна для всех многочленов степени N и ниже, 2) среди многочленов степени N + 1 найдется хотя бы один, для которого формула неточна. С учетом этого определения и обозначения N теорему можно переформулировать так. Для того, чтобы квадратурная формула была интерполяционной, необходимо и достаточно, чтобы N ≥ n. 8.2. Погрешность интерполяционных квадратурных формул Для краткости будем применять следующие обозначения Zb f (x)dx, S[f ] = I[f ] = n X Ak f (xk ), k=0 a тогда по определению погрешность интегрирования Rинт [f ] = I[f ] − S[f ], а для интерполяционных квадратур, у которых S[f ] получено интегрированием интерполяционного многочлена Rинт [f ] = I[f ] − I[Ln ] = I[f − Ln ] = I[Rn ]. Подставляя формулу для погрешности интерполяции получаем Zb Rинт [f ] = f (n+1) (ξ) ωn (x)dx (n + 1)! a Заметим, что ξ зависит от x, поэтому оценить интеграл в общем случае затруднительно. Однако, если мы перейдем к модулю погрешности чис- 85 ленного интегрирования, то получим оценку Mn+1 |Rинт [f ]| ≤ (n + 1)! Zb |ωn (x)|dx a 8.3. Элементарные квадратурные формулы (Формулы Ньютона-Котеса) Элементарные квадратурные формулы (Формулы Ньютона-Котеса) – это широко распространенные интерполяционные квадратуры при небольшом количестве узлов и при их стандартном расположении. Их изучение начнем с n = 0 (один узел). Формулы прямоугольников. Рассмотрим три случая. А. x0 = a. Тогда L0 (x) = f (x0 ) = f (a), Zb f (a)dx = f (a)(b − a) a и формула имеет вид Zb f (x)dx ≈ f (a)(b − a), a эта формула называется элементарной формулой левых прямоугольников. Алгебраическая степень точности N = 0, как можно проверить по определению. Замечание. Геометрическая интерпретация этого и других рассмотренных ниже методов обычно рассматривается на практике, поэтому рисунки с иллюстрациями методов приводятся в методической разработке по практическим и лабораторным работам дисциплины "Численные методы". 86 Б. x0 = b. Тогда L0 (x) = f (b), Zb f (b)dx = f (b)(b − a) a и формула имеет вид Zb f (x)dx ≈ f (b)(b − a), a эта формула называется элементарной формулой правых прямоугольников. Алгебраическая степень точности N = 0. В. x0 = a+b 2 . Тогда L0 (x) = f ( a+b 2 ), Zb f (a)dx = f ( a+b )(b − a) 2 a и формула имеет вид Zb f (x)dx ≈ f ( a+b )(b − a), 2 a эта формула называется элементарной формулой средних прямоугольников. Алгебраическая степень точности N = 1. Это можно проверить, подставляя в формулу последовательно 1, x, x2 . Формула трапеций. Пусть теперь n = 1 (два узла), x0 = a, x1 = b. Тогда L1 (x) = f (a) + f (a, b)(x − a) = f (a) + Zb (f (a) + f (b)−f (a) b−a (x − a), f (b) − f (a) (x − a))dx = f (a)(b − a)+ b−a a + f (b) − f (a) (b − a)2 f (a) + f (b) = (b − a), b−a 2 2 87 и формула имеет вид Zb f (x)dx ≈ f (a) + f (b) (b − a). 2 a Эта формула называется элементарной формулой трапеций. Алгебраическая степень точности N = 1. Формула Симпсона. Пусть n = 2 (три узла), x0 = a, x1 = a+b x2 = b. Тогда L2 (x) = f (a) + f (a, a+b 2 )(x − a) + f (a, 2 , b)(x − a)(x − a+b 2 , a+b 2 ). Чтобы не проделывать дальнейшее громоздкое расписывание разделенных разностей и последующее интегрирование, воспользуемся методом неопределенных коэффициентов. Запишем формулу Zb f (x)dx ≈ Af (a) + Bf ( a+b ) + Cf (b). 2 a Подставим f (x) ≡ 1, получим равенство b−a=A+B+C Подставим f (x) = x − a, получим равенство (b − a)2 (b − a) =B + C(b − a) 2 2 Подставим f (x) = (x − a)2 , получим равенство (b − a)3 (b − a)2 =B + C(b − a)2 3 4 Решая эту систему, получим A = 16 (b − a), B = 64 (b − a), C = 16 (b − a), т.е. формула приобретает вид Zb f (x)dx ≈ b−a a+b (f (a) + 4f ( ) + f (b)), 6 2 a 88 эта формула называется элементарной интерполяционной формулой парабол или элементарной формулой Симпсона. Алгебраическая степень точности N = 3. Это можно проверить, подставляя в формулу последовательно 1, x, x2 , x3 , x4 . Формула "3/8". Пусть n = 3, четыре равноотстоящих узла, расстояние между которыми b−a 3 , тогда x0 = a, x1 = 2a+b 3 , x2 = a+2b 3 , x3 = b. Запишем формулу с неопределёнными коэффициентами Zb f (x)dx ≈ Af (a) + Bf ( 2a + b a + 2b ) + Cf ( ) + Df (b). 3 3 a Коэффициенты определяются, подставляя многочлены низших степеней, в результате A = D = b−a 8 , B=C= 3(b−a ) 8. Полученная формула Zb f (x)dx ≈ 2a + b a + 2b b−a (f (a) + 3f ( ) + 3f ( ) + f (b)) 8 3 3 a называется формулой "3/8". Заметим, что с увеличением числа узлов элементарные интерполяционные квадратурные формулы всё усложняются, что создает препятствия для их практического применения. Для повышения точности получил распространение другой подход – составные квадратурные формулы. 8.4. Погрешность элементарных интерполяционных квадратурных формул Рассмотрим погрешность элементарной формулы левых прямоугольников. Как уже отмечалось, для интерполяционных квадратурных формул Zb Rинт [f ] = I[Rn ] = Rn (x)dx, a 89 в частности, для элементарной формулы левых прямоугольников Zb Zb Rинт [f ] = I[Rn ] = f (x, a)(x − a)dx. R0 (x)dx = a a Воспользуемся известной из математического анализа теоремой о среднем, которая утверждает, что Zb Zb h(x)dx, c ∈ [a, b], g(x)h(x)dx = g(c) a a если g(x) – непрерывная функция, h(x) – непрерывная функция, сохраняющая знак на отрезке [a, b]. Тогда, т.к. (x − a) ≥ 0 на отрезке [a, b], то по теореме о среднем Zb Zb f (x, a)(x − a)dx = f (c, a) a (b − a)2 , c ∈ [a, b], (x − a)dx = f (c, a) 2 a а используя связь разделенной разности и производной, получаем формулу для погрешности (b − a)2 Rинт [f ] = f (η) , η ∈ [a, b]. 2 и её оценку |Rинт [f ]| ≤ M1 (b − a)2 . 2 Аналогичным образом выводится погрешность элементарной формулы правых прямоугольников (b − a)2 Rинт [f ] = −f (η) , η ∈ [a, b]. 2 и её оценка |Rинт [f ]| ≤ M1 (b − a)2 . 2 Рассмотрим погрешность элементарной формулы трапеций. 90 Zb Rинт [f ] = I[Rn ] = Zb f (x, a, b)(x − a)(x − b)dx. R1 (x)dx = a a Используя теорему о среднем (т.к. (x − a)(x − b) ≤ 0 на отрезке [a, b]) и считая интеграл, получаем Zb Zb f (x, a, b)(x − a)(x − b)dx = f (c, a, b) a (x − a)(x − b)dx = a = −f (c, a, b) (b − a)3 , c ∈ [a, b] 6 Используя связь разделенной разности второго порядка и соответствующей производной, получаем формулу для погрешности (b − a)3 Rинт [f ] = −f (η) , η ∈ [a, b]. 12 00 и её оценку |Rинт [f ]| ≤ M2 (b − a)3 . 12 Рассмотрим погрешность элементарной формулы средних прямоугольников. Прием, основанный на прямом интегрировании погрешности интерполяции не проходит, т.к. теорема о среднем неприменима, ибо x − a+b 2 не сохраняет знак на отрезке. Оценка модуля погрешности на этом пути может быть получена, но она слишком завышена. Применим прием, основанный на повышенной алгебраической степени точности формулы N = 1. Пусть H1 (x) – интерполяционный многочлен Эрмита, построенный по одному двухкратному узлу x0 = a+b 2 , т.е. построенный по интерполяци- 0 a+b онным данным f ( a+b 2 ), f ( 2 ). Так как в узле x0 = a+b 2 , функции f (x) и H1 (x) совпадают, то для формулы средних прямоугольников S[f ] = S[H1 ]. 91 Так как H1 (x) – многочлен первой степени, а алгебраическая степень точности формулы N = 1, то I[H1 ] = S[H1 ]. Получаем цепочку соотношений Rинт [f ] = I[f ] − S[f ] = I[f ] − S[H1 ] = I[f ] − I[H1 ] = = I[f − H1 ] = I[R1 ]. Посчитаем интеграл от погрешности интерполяции многочленом H1 (x), используя теорему о среднем Zb Rинт [f ] = I[R1 ] = f (x, a+b a+b a+b 2 , )(x − ) dx = 2 2 2 a a+b a+b = f (c, , ) 2 2 Zb a+b 2 f 00 (η) (b − a)3 (x − ) dx = , η ∈ [a, b] 2 2 12 a Таким образом получаем формулу для погрешности (b − a)3 , η ∈ [a, b]. Rинт [f ] = f (η) 24 00 и её оценку |Rинт [f ]| ≤ M2 (b − a)3 . 24 Аналогичный прием применим для вывода погрешности элементарной формулы Симпсона, которая обладает повышенной алгебраической степени точности формулы N = 3. Пусть H3 (x) – интерполяционный многочлен Эрмита, построенный по двухкратному узлу x1 = a+b 2 , и однократным узлам x0 = a и x2 = b 0 a+b т.е. построенный по интерполяционным данным f (a), f ( a+b 2 ), f ( 2 ), f (b). Так как в узлах f (x) и H3 (x) совпадают, то для формулы Симпсона S[f ] = S[H3 ]. Так как H3 (x) – многочлен третьей степени, а алгебраическая степень точности формулы N = 3, то I[H3 ] = S[H3 ]. 92 Получаем цепочку соотношений Rинт [f ] = I[f ] − S[f ] = I[f ] − S[H3 ] = I[f ] − I[H3 ] = = I[f − H3 ] = I[R3 ]. Посчитаем интеграл от погрешности интерполяции многочленом H1 (x), используя теорему о среднем Zb Rинт [f ] = f (x, a, a+b a+b a+b 2 , , b)(x − a)(x − b)(x − ) dx = 2 2 2 a a+b a+b = f (c, a, , , b) 2 2 Zb (x − a)(x − b)(x − a+b 2 ) dx = 2 a f 0000 (η) (b − a)5 , η ∈ [a, b] 24 120 Последний интеграл лучше считать по частям. =− Таким образом получаем формулу для погрешности (b − a)5 , η ∈ [a, b]. Rинт [f ] = −f (η) 2880 0000 и её оценку M4 (b − a)5 . 2880 |Rинт [f ]| ≤ 8.5. Составные квадратурные формулы Составные квадратурные формулы получаются путем разбиения отрезка интегрирования на равные части и применения на каждом отрезке одинаковых элементарных квадратурных формул. Рассмотрим этот процесс на примере составной формулы левых прямоугольников. Пусть требуется посчитать Zb I[f ] = f (x)dx. a 93 Разобьем отрезок интегрирования на m частей с шагом h = b−a m , введем точки xi = a+ih, i = 0, 1, . . . , m разобьем интеграл на сумму m интегралов и для каждого применим элементарную формулу левых прямоугольников: Zb f (x)dx = m−1 X xi+1 Zxi+1 m−1 m−1 XZ X f (x)dx ≈ f (xi )dx = h f (xi ) i=0 x i i=0 x i a т.е. I[f ] ≈ h m−1 X i=0 f (xi ). i=0 Аналогичным образом выводится составная формула правых прямоугольников I[f ] ≈ h m X f (xi ), i=1 составная формула средних прямоугольников I[f ] ≈ h m X i=0 h f (xi + ). 2 и составная формула трапеций m−1 f (a) + f (b) X I[f ] ≈ h( + f (xi )). 2 i=1 Выведем составную формулу Симпсона m−1 X Zxi+1 m−1 h hX f (x)dx ≈ (f (xi ) + 4f (xi + ) + f (xi+1 )) 6 i=0 2 i=0 x i Приводя подобные, получаем m−1 m−1 X X h h f (xi ) + 4 f (xi + )) I[f ] ≈ (f (a) + f (b) + 2 6 2 i=1 i=0 94 8.6. Погрешность составных квадратурных формул Для того, чтобы получить погрешность составной квадратурной формулы, нужно просуммировать погрешности элементарных квадратурных формул на каждом отрезке. Так для составной формулы левых прямоугольников Rинт [f ] = m−1 X i Rинт [f ] i=0 = m−1 X i=0 h2 f (ηi ) , ηi ∈ [xi , xi+1 ]. 2 Отсюда, используя соотношение mh = b − a получаем оценку погрешности |Rинт [f ]| ≤ M1 (b − a) h, M1 = max |f 0 (x)|. 2 x∈[a,b] Эта оценка показывает, что метод сходится, т.е. погрешность стремится к нулю при стремлении шага h к нулю, если подинтегральная функция непрерывно дифференцируема. Та же оценка погрешности справедлива для составной формулы правых прямоугольников. Для составной формулы средних прямоугольников справедлива оценка погрешности |Rинт [f ]| ≤ M2 (b − a) 2 h , M2 = max |f 00 (x)|, 24 x∈[a,b] а для составной формулы трапеций |Rинт [f ]| ≤ M2 (b − a) 2 h , M2 = max |f 00 (x)|, 12 x∈[a,b] эти две формулы имеют второй порядок сходимости. Составная формула Симпсона имеет оценку погрешности |Rинт [f ]| ≤ M4 (b − a) 4 h , M4 = max |f 0000 (x)|. 2880 x∈[a,b] 95 8.7. Метод Рунге практической оценки погрешности Выведенные оценки погрешностей на практике применять затруднительно, т.к. необходимо делать оценки соответствующих производных. Изложим другой метод оценки погрешностей, который позволяет по двум вычисленным значениям искомой величины (в данном случае интеграла) оценивать погрешность и уточнять результат. Пусть требуется вычислить некоторую идеальную величину I, которая приближенно считается с помощью величины Sh , зависящей от параметра h. Обозначим погрешность через Rh , таким образом I = Sh + Rh Предположим, что погрешность представляется в виде асимптотического разложения Rh = Chp + o(hp ), где известен порядок p. Вычислим Sh при двух значения параметра h, например, Sh и Sh/2 . Из системы уравнений   I = Sh + Rh  I = Sh/2 + Rh/2 исключим неизвестную величину I, получим связь Sh/2 − Sh + Rh/2 − Rh = 0. С другой стороны, отбрасывая слагаемые более высокого порядка чем hp , имеем представление Rh ≈ Chp , Rh/2 ≈ C(h/2)p , 96 т.е. Rh ≈ 2p Rh/2 . Отсюда получаем Rh/2 ≈ Sh/2 − Sh , 2p − 1 2p Rh ≈ p (Sh/2 − Sh ), 2 −1 эти формулы называются формулами Рунге практической оценки погрешности. Они не только позволяют оценивать погрешность, но и уточнять результат: Sh/2 − Sh 2p − 1 Рассмотрим применение этого метода к одной из составных формул, I ≈ Sh/2 + скажем к составной формуле левых прямоугольников. Сначала нужно доказать, что справедливо асимптотическое разложение погрешности. Имеем Rh = Rинт [f ] = m−1 X i=0 m−1 h2 h X 0 f (ηi ) = ( f (ηi )h) = 2 2 i=0 Zb h = ( f 0 (x)dx + O(h)), 2 a т.е. для составной формуле левых прямоугольников имеет место асимптотическое разложение с p = 1. Формулы Рунге для составной формуле левых прямоугольников, также как и для составной формуле правых прямоугольников Rh/2 ≈ Sh/2 − Sh , Rh ≈ 2(Sh/2 − Sh ). Для составной формулы трапеций и составной формулы средних прямоугольников, где p = 2 формулы Рунге имеют вид Rh/2 ≈ (Sh/2 − Sh ) , 3 97 Rh ≈ 4(Sh/2 − Sh ) . 3 Для составной формулы Симпсона, где p = 4 формулы Рунге имеют вид (Sh/2 − Sh ) , 15 16(Sh/2 − Sh ) Rh ≈ . 15 Rh/2 ≈ 8.8. Формулы наивысшей алгебраической степенью точности В дальнейшем будем рассматривать задачу о вычислении интеграла с весом Zb I[f ] = p(x)f (x)dx, a с помощью квадратурных формул Zb p(x)f (x)dx ≈ n X Ak f (xk ) k=0 a Весовая функция, которая входит в коэффициенты, используется, например, при вычислении интегралов с особенностями, см. ниже. Интерполяционные квадратурные формулы обладают алгебраической степенью точности N ≥ n, так как их коэффициенты Ak выбираются особым образом, раздел . Но можно выбирать и узлы интерполяции xk , что еще может повысить алгебраическую степень точности. Пример. Рассмотрим формулу Z1 √ √ 3 3 f (x)dx ≈ f (− ) + f ( ) 3 3 −1 98 Её алгебраическая степень точности (как можно проверить) N = 3, а у аналогичной (с тем же количеством узлов) формулы трапеций Z1 f (x)dx ≈ f (−1) + f (1) −1 алгебраическая степень точности N = 1. Перейдем к нахождению условий, обеспечивающих наивысшую алгебраическую точность. Сначала заметим, что общее число параметров (коэффициентов и узлов) 2n + 2, они обеспечивают условия, а многочлен степени 2n+1 содержит столько же коэффициентов. Следует ожидать, что N = 2n + 1. Это в самом деле так при определённых условиях, но доказать этот факт непросто. Справедливо утверждение. Теорема 14. Для того, чтобы квадратурная формула с весом была точна для любого многочлена степени 2n+1 и ниже, необходимо и достаточно, чтобы а) коэффициенты находились по формулам Zb Ak = p(x) a n Y j=0,j6=k x − xj dx xk − xj б) узлы такие, что многочлен ωn (x) = (x − x0 )(x − x1 ) · · · (x − xn ) ортогонален с весом любому многочлену g(x) степени n и ниже. Ортогональность с весом многочленов ωn (x) и g(x) означает, что Zb p(x)ωn (x)g(x)dx = 0. a Необходимость. Пусть квадратурная формула с весом точна для любого многочлена степени 2n + 1 и ниже, возьмем в качестве такого много99 члена n Y Qi (x) = j=0,j6=i x − xj xi − xj и получим формулу Zb p(x)Qi (x)dx = n X Ak Qi (xk ) k=0 a Но Qi (xk ) = 0, k 6= i, Qi (xi ) = 0, и последняя сумма вырождается в одно слагаемое Ai . Подставляя в эту формулу Qi (x), получаем Zb p(x) n Y j=0,j6=i a x − xj dx = Ai . xi − xj Возьмем любой многочлен g(x) степени не выше n, тогда произведение ωn (x)g(x) – многочлен, степени не выше 2n + 1, следовательно для этого произведения квадратура с весом точна и, учитывая что ωn (xk ) = 0, получаем Zb p(x)ωn (x)g(x)dx = n X Ak ωn (xk )g(xk ) = 0 k=0 a Достаточность. Сначала заметим, что т.к. коэффициенты Ak выбираются из условия замены интегрируемой функции её интерполяционным многочленом, то квадратура с весом точна для любого многочлена, степени n и ниже. Возьмем любой многочлен Q(x), степени не выше 2n + 1, поделим его на ωn (x), т.е. представим в виде Q(x) = g(x)ωn (x) + r(x), где многочлены g(x) и r(x) имеют степени не выше n. Тогда, в силу условий а) и б) Zb Zb p(x)Q(x)dx = a Zb p(x)g(x)ωn (x)dx + a p(x)r(x)dx = a 100 = n X Ak r(xk ) = n X Ak Q(xk ) k=0 k=0 Доказанная теорема дает условия для формул наивысшей алгебраической степени точности (квадратуры Гаусса), но не решает вопрос о её существовании и единственности. 8.9. Существование и единственность квадратуры Гаусса Представим многочлен ωn (x) в виде разложения по степеням x: ωn (x) = xn+1 + a1 xn + . . . + an x + an+1 = Ω(x) Обозначения Ω(x) будем применять, т.к. класс таких произвольных многочленов степени n + 1 шире, чем ранее рассмотренные многочлены ωn (x). Распишем условие ортогональности с весом многочлена Ω(x) степеням xl , l = 0, 1, . . . , n: Zb p(x)(xn+1 + a1 xn + . . . + an x + an+1 )xl dx = 0, l = 0, 1, . . . , n a Получили линейную неоднородную систему относительно коэффициентов ai . Эта система имеет единственное решение тогда и только тогда, когда соответствующая однородная система Zb p(x)(a1 xn + . . . + an x + an+1 )xl dx = 0, l = 0, 1, . . . , n a имеет только тривиальное решение. Умножим l-ое уравнение последней системы на al и сложим все уравнения, получим одно Zb p(x)(a1 xn + . . . + an x + an+1 )2 dx = 0. a 101 Установим условия, при котором это уравнение имеет лишь тривиальное решение. Для этого рассмотрим вспомогательное утверждение. Лемма. Пусть p(x) непрерывная неотрицательная функция, такая, что вия Rb a Rb p(x)dx > 0, а многочлен Q(x) неотрицателен на [a, b]. Тогда из услоp(x)Q(x)dx = 0 следует Q(x) ≡ 0. a Доказательство. Так как Rb p(x)dx > 0, а функция p(x) непрерывная a неотрицательная, то найдется отрезок [a1 , b1 ] ⊂ [a, b], такой, что p(x) > 0 для всех x ∈ [a1 , b1 ]. Так как Q(x) многочлен, то он имеет на отрезке [a1 , b1 ] конечное число корней xi , i = 1, . . . , m. Окружим каждый корень интервалом m S (x − ε, x + ε) (ε > 0 мало) и введем множество Xε = [a1 , b1 ]\ (x − ε, x + ε) i=1 На множестве ограниченном замкнутом множестве Xε положительная непрерывная функция Q(x) ограничена снизу: найдется δ > 0, такое, что Q(x) ≥ δ для всех x ∈ [a1 , b1 ]. Из равенства Zb Z p(x)Q(x)dx = 0= a Z p(x)Q(x)dx + Xε p(x)Q(x)dx [a,b]\Xε и неотрицательности подинтегральных функций следует, что Z p(x)Q(x)dx = 0 Xε Но тогда Z Z p(x)Q(x)dx ≥ δ 0= Xε p(x)dx > 0, Xε получили противоречие, следовательно Q(x) имеет бесконечное число корней, т.е. Q(x) – нулевой многочлен. Если применить результат этой леммы к рассматриваемой задаче, то получаем утверждение 102 Теорема 15. Пусть вес p(x) непрерывная неотрицательная функция, Rb такая, что p(x)dx > 0, тогда существует единственный многочлен a Ω(x) = xn+1 + a1 xn + . . . + an x + an+1 , ортогональный с весом степеням xl , l = 0, 1, . . . , n. Теорема 16. В условиях предыдущей теоремы многочлен Ω(x) имеет ровно n + 1 различных вещественных корней на отрезке [a, b]. Доказательство состоит из двух частей. Сначала докажем, что многочлен Ω(x) имеет n + 1 корней на отрезке [a, b]. Предположим противное, т.е. что корней меньше, чем n + 1, обозначим корни x0 , x1 , . . . , xs , s ≤ n − 1. Тогда Ω(x) = (x − x0 )(x − x1 ) · · · (x − xs )r(x), где r(x) не обращается в ноль на отрезке [a, b], для определённости будем считать, что r(x) > 0. Обозначим g(x) = (x−x0 )(x−x1 ) · · · (x−xs ), степень этого многочлена s + 1 ≤ n, воспользуемся ортогональностью с весом многочлену Ω(x), получим Zb p(x)g 2 (x)r(x)dx = 0, a применяя лемму, получим r(x) ≡ 0, противоречие, т.е. многочлен Ω(x) имеет n + 1 корней на отрезке [a, b]. Теперь докажем, что у многочлена Ω(x) кратных корней нет. Предположим противное, т.е. что имеется корень, скажем x0 кратности 2 (или выше). Тогда Ω(x) = (x − x0 )2 Q(x), где степень многочлена Q(x) равна n − 1, воспользуемся ортогональностью 103 Q(x) с весом многочлену Ω(x), получим Zb p(x)(x − x0 )2 Q2 (x)dx = 0, a применяя лемму, получим Q(x) ≡ 0, противоречие, т.е. многочлен Ω(x) не имеет кратных корней. Теоремы 15 и 16 дают условия существования и единственности квадратуры Гаусса. 8.10. Алгоритм построения квадратуры Гаусса Подведем итог изложенному выше в виде алгоритма для построения квадратуры Гаусса. Этап 1. Составляем систему Zb p(x)(xn+1 + a1 xn + . . . + an x + an+1 )xl dx = 0, l = 0, 1, . . . , n a и вычисляем её коэффициенты. Этап 2. Решаем эту систему, определяя ai , i = 0, 1, . . . , n. Этап 3. Находим корни xk , i = 0, 1, . . . , n многочлена ωn (x) = xn+1 + a1 xn + . . . + an x + an+1 . Этап 4. Находим коэффициенты Ak , i = 0, 1, . . . , n квадратуры Гаусса Zb p(x)f (x)dx ≈ n X Ak f (xk ), k=0 a например, методом неопределенных коэффициентов. Проиллюстрируем примером. Построим квадратуру Гаусса вида Z1 f (x)dx ≈ A0 f (x0 ) + A1 f (x1 ). −1 104 Так как n = 1, составим многочлен ω1 (x) = x2 +a1 x+a2 и запишем систему из условия ортогональности этого многочлена 1 и x:  R1 2     (x + a1 x + a2 )dx = 0 −1 R1  2    (x + a1 x + a2 )xdx = 0 −1 Считая интегралы, получаем   2a2 dx = 2 3   2 a1 = 0 3 2 Находим корни уравнения x + √ 1 3 = 0, получаем x0 = − 3 3 , √ x1 = 3 3 . Записываем формулу Z1 √ √ 3 3 f (x)dx ≈ А0 f (− ) + А1 f ( ) 3 3 −1 и методом неопределенных коэффициентов находим A0 = AA = 1. Алгебраическая степень точности N = 2n + 1 = 3. 8.11. Погрешность квадратуры Гаусса Будем, как и раньше, для краткости обозначать Zb I[f ] = p(x)f (x)dx, S[f ] = n X Ak f (xk ), k=0 a Пусть H2n+1 (x) – интерполяционный многочлен Эрмита, построенный по двукратным узлам x0 , x1 , . . . , xn , т.е. построенный по интерполяционным данным f (x0 ),f 0 (x0 ), . . . , f (xn ), f 0 (xn ). Так как в узлах функции f (x) и H2n+1 (x) совпадают, то для квадратуры Гаусса S[f ] = S[H2n+1 ]. Так как H2n+1 (x) – многочлен степени 2n + 1, а алгебраическая степень точности формулы N ≥ 2n + 1, то I[H2n+1 ] = S[H2n+1 ]. 105 Получаем цепочку соотношений для погрешности интегрирования Rинт [f ] = I[f ] − S[f ] = I[f ] − S[H2n+1 ] = I[f ] − I[H2n+1 ] = = I[f − H2n+1 ] = I[R2n+1 ]. Посчитаем интеграл от погрешности интерполяции многочленом H2n+1 (x), используя теорему о среднем Zb Rинт [f ] = p(x)f (x, x0 , x0 , · · · , xn , xn )(x − x0 )2 · · · (x − xn )2 dx = a Zb = f (c, x0 , x0 , · · · , xn , xn ) p(x)ωn2 (x)dx = a f (2n+2) (η) = (2n + 2)! Zb p(x)ωn2 (x)dx, η ∈ [a, b] a Таким образом, получаем формулу для погрешности f (2n+2) (η) Rинт [f ] = (2n + 2)! Zb p(x)ωn2 (x)dx, η ∈ [a, b] a и её оценку M2n+2 |Rинт [f ]| ≤ (2n + 2)! Zb p(x)ωn2 (x)dx. a Из этих оценок вытекают интересные следствия. Замечание 1. Алгебраическая степень точности квадратуры Гаусса N = 2n + 1. Нужно показать, что для многочленов степени 2n + 2 формула не точна. Возьмем такой многочлен Q(x) = x2n+2 . Тогда его производная поRb рядка 2n + 2 равна (2n + 2)!, а p(x)ωn2 (x)dx > 0, т.е. Rинт [Q] 6= 0. a Замечание 2. Коэффициенты квадратуры Гаусса Ak положительны. 106 Для доказательства возьмем многочлен Qi (x) = [ n Y j=0,j6=i x − xj 2 ] xi − xj степени 2n + 2. Тогда Zb 0< p(x)Qi (x)dx = n X Ak Qi (xk ) = Ai , k=0 a что и требовалось доказать. 8.12. Вычисление интегралов с особенностями Под интегралами с особенностями понимаются интегралы от неограниченных функций или интегралы на бесконечных промежутках интегрирования. Кроме того, особенностью называется недостаточная гладкость подинтегральной функции при применении той или иной квадратурной формулы, так, например, для применения составной формулы Симпсона требуется непрерывность четвертой производной. Рассмотрим на примерах некоторые приемы, позволяющие устранять особенности. Аналитические способы. 1. Требуется вычислить ZΠ/2 ln sin xdx, имеющий особенность в нуле. Добавим и вычтем sin x к подинтегральной функции и получим ZΠ/2 ZΠ/2 ZΠ/2 sin x ln sin xdx = ln xdx + ln dx. x 107 Первый интеграл может быть вычислен аналитически, а второй особенности не содержит и может быть вычислен численно. 2. Чтобы сдвинуть особенность в старшую производную, часто применяют формулу интегрирования по частям. Рассмотрим Za xα g(x)dx, где функция g(x) "хорошая" (достаточное число раз дифференцируемая). При −1 < α < 0 подинтегральная функция неограничена, при 0 < α < 1 подинтегральная функция ограничена, но её первая производная неограничена и так далее. Проинтегрируем по частям Za 1 1 xα g(x)dx = xα+1 g(x)|a0 − α+1 α+1 Za xα+1 g 0 (x)dx, таким образом сдвинули особенность на единицу в старшую производную. Метод усечений. Суть метода состоит в том, что мы отделяемся от особенности с нужной точностью. В качестве примера рассмотрим задачу о вычислении Z∞ 2 е −x dx, 2 с точностью ε = 0, 0001. Имеем Z∞ 2 е −x dx = 2 ZB 2 е −x dx + 2 Z∞ 2 е −x dx, B Последний интеграл дает погрешность усечения. Сделаем оценку Z∞ B 2 е −x dx ≤ Z∞ 2x −x2 1 −A2 е dx = е < ε/2. 2A 2A B 108 Вычисления показывают, что при B = 3 эта оценка выполняется. Оценка погрешности составной формулы Симпсона показывает, что если взять n = 8, то выполняется Z3 2 е −x dx < ε/2. 2 Введение весовой функции. Мультипликативный метод выделения особенностей. Если подинтегральная функция есть произведение двух функций, "хорошей" и содержащей особенность, то имеет смысл рассматривать формулы численного интегрирования с весом, вес включается в коэффициенты квадратуры. Например, для вычисления Z1 √ xf (x)dx можно использовать интерполяционные квадратурные формулы с весом Z1 √ xf (x)dx ≈ n X Z1 Ak f (xk ), Ak = k=0 √ x n Y j=0,j6=k x − xj dx xk − xj К этому же классу методов относятся специальные формулы для быстро осцилирующих функций. Требуется вычислить Zb g(x) sin ωxdx. a Попробуем применить какую-либо составную формулу, например, составную формулу средних прямоугольников Zb g(x) sin ωxdx ≈ h a n−1 X i=0 109 g(xi+1/2 ) sin ωxi+1/2 В оценку погрешности этой формулы входит производная, однако f 0 (x) = (g(x) sin ωx)0 = g 0 (x) sin ωx + ωg(x) cos ωx и эта функция при больших ω может принимать большие значения. Выведем аналог составной формулы средних прямоугольников заменяя константой в средней точке не всю подинтегральную функцию, а только g(x), включая sin ωx в весовую функцию. Zb x n−1 Zi+1 X g(x) sin ωxdx = g(x) sin ωxdx ≈ a i=0 x i ≈ n−1 X Zxi+1 g(xi+1/2 ) sin ωxdx = i=0 xi n−1 1X =− g(xi+1/2 )(cos ωxi+1 − cos ωxi ) = ω i=0 n−1 2 hX = sin ω g(xi+1/2 ) sin ωxi+1/2 ω 2 i=0 При этом погрешность оценивается величиной x n−1 Zi+1 X |Rинт [f ]| ≤ |g(x, xi+1/2 )||x − xi+1/2 || sin ωx|dx ≤ i=0 x i ≤ M1 [g](b − a) h, M1 [g] = max |g 0 (x)| 4 x∈[a,b] которая стремится к нулю при шаге h стремящемся к нулю равномерно по ω. Аддитивный метод устранения особенностей. Идея метода состоит в том, чтобы представить исходную подинтегральную функцию в виде суммы f (x) = ϕ(x) + ψ(x), где интеграл от ϕ(x), содержащей те же особенности, что и исходная функция, считается 110 аналитически, а функция ψ(x) не содержит особенности и интеграл от неё считается численно. Основным инструментом при этом является тейлоровское разложение в окрестности особенности. Рассмотрим пример. Требуется сдвинуть особенность в третью производную для подсчета Z1 sin x dx. x1/2 (1 + x2 /2) Разложим функции в нуле x3 sin x = x − + ··· 6 x2 x4 1 = 1 − + + ··· (1 + x2 /2) 2 4 sin x x5/2 x2 x4 1/2 = (x − + · · · )(1 − + + ···) = 6 2 4 x1/2 (1 + x2 /2) 2 = x1/2 − x5/2 + · · · 3 Обозначим 2 ϕ(x) = x1/2 − x5/2 3 sin x 2 5/2 1/2 ψ(x) = 1/2 − x + x 3 x (1 + x2 /2) Интеграл ϕ(x) считается аналитически, а функция ψ(x) по крайней мере дважды непрерывно дифференцируема и интеграл от неё может быть вычислен, например, с помощью составной формулы трапеций. 111 ЛИТЕРАТУРА 1. Бабенко К.И. Основы численного анализа. М.: Наука, 1986. 452 с. 2. Бахвалов Н.С. Численные методы. М.: Наука, 1973. 632 с. 3. Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. М.: Наука, 1987. 598 с. 4. Березин И.С., Жидков Н.П. Методы вычислений. Т. 1. М.: Физматгиз, 1959. 464 с. 5. Березин И.С., Жидков Н.П. Методы вычислений. Т. 2. М.: Наука, 1966. 430 с. 6. Вержбицкий В.М. Численные методы. Линейная алгебра и нелинейные уравнения. М.: ОНИКС 21 век, 2005. 432 с. 7. Вержбицкий В.М. Численные методы. Математический анализ и обыкновенные дифференциальные уравнения. М.: ОНИКС 21 век, 2005. 400 с. 8. Волков Е.А. Численные методы. М: Наука, 1982. 248 с. 9. Калиткин Н.Н. Численные методы. 2-е издание. СПб.: БХВПетербург, 2011. 586 с. 10. Крылов В.И., Бобков В.В., Монастырный П.И. Начала теории вычислительных методов (в 5 томах). Интерполирование и интегрирование. Минск, Наука и техника, 1983. 287 с. 11. Крылов В.И., Бобков В.В., Монастырный П.И. Начала теории вычислительных методов (в 5 томах). Линейная алгебра и нелинейные уравнения. Минск, Наука и техника, 1985. 279 с. 112 12. Марчук Г.И. Методы вычислительной математики. М.: Наука, 1989. 608 с. 13. Петров И.Б., Лобанов А.И. Лекции по вычислительной математике. М.: БИНОМ, 2006. 524 с. 14. Самарский А.А., Гулин А.В. Численные методы. М.: Наука, 1989. 430 с. 15. Фадеев А.К., Фадеева В.Н. Вычислительные методы линейной алгебры. СПб.: Лань, 2002. 736 с. 16. Формалев В.Ф., Ревизников Д.Л. Численные методы. М.: ФИЗМАТЛИТ, 2006. 400 с. Библиографический комментарий Учебную литературу условно можно разбить на несколько групп. К первой группе относятся классические университетские учебники [4, 5, 2, 3, 14, 9, 10, 11], данный курс читается ближе всего к учебникам [14, 9]. Ряд книг [12, 1, 15] предназначен, прежде всего, для специалистов и они содержат много дополнительных сведений. Имеется очень много учебников для технических университетов, особенно появившихся в последнее время, отметим некоторые [6, 7, 13, 16], которые содержат описание разнообразных методов, но зачастую в них отсутствуют доказательства, например, совсем простой учебник [8]. 113