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

Лекция 1 Введение. Типы и источники погрешностей при приближенных вычислениях Целью курса «Численные методы» является рассмотрение вопросов, связанных с приближенными методами решения различных прикладных задач. Численные методы разрабатывают и исследуют, как правило, высококвалифицированные специалисты по вычислительной математике. Что же касается прикладников, использующих численные методы как готовый инструмент в своей практической деятельности, то для них главной задачей является понимание основных идей методов, их особенностей и областей их применения. Обычно решение прикладной задачи состоит из нескольких этапов. На первом этапе формулируется содержательная постановка задачи, которую требуется решить. При этом решается вопрос, какие факторы следует учесть, а какими можно пренебречь. На втором этапе содержательной модели ставится в соответствие математическая модель, т.е. математическое описание естественного процесса с помощью различного рода уравнений и неравенств. Третий этап состоит в построении приближенного метода решения математической задачи, т.е. в выборе вычислительного алгоритма. Специалисту-прикладнику для решения задачи, как правило, необходимо выбрать из имеющегося арсенала методов тот, который наиболее пригоден в каждом конкретном случае. На четвертом этапе осуществляется программирование вычислительного алгоритма для ЭВМ. На пятом этапе готовятся исходные данные задачи и проводятся расчеты по составленной и отлаженной программе. В курсе численного анализа мы в основном будем иметь дело с вопросами, касающимися третьего этапа, т.е. с численными методами решения математических задач. При решении задачи на ЭВМ мы всегда получаем не точное, а приближенное решение. Мерой точности решения является допущенная при его получении погрешность. Можно выделить три основные причины возникновения погрешностей при численном решении задачи. Прежде всего, исходные данные задачи (начальные и граничные условия, коэффициенты и правые части уравнений и т.п.) всегда задаются с некоторой погрешностью. Эту погрешность принято называть неустранимой. Численный метод, используемый при решении математической задачи, также является источником погрешностей. Это связано, например, с заменой интегралов конечными суммами, усечением бесконечных рядов при вычислении значений функций и т.п. Как правило, погрешность метода регулируема, т.е. может быть уменьшена до некоторого разумного предела путем изменения параметров численного метода, например, уменьшения шага при численном интегрировании и т.д. Наконец, в силу ограниченности разрядной сетки ЭВМ, т.е. числа разрядов, предназначенных для хранения в ЭВМ чисел, неизбежными являются погрешности округления. В процессе выполнения алгоритма погрешности округления накапливаются. Результирующая погрешность, вносимая в решение за счет погрешностей округления, называется вычислительной. Величина вычислительной погрешности определяется двумя факторами: точностью представления чисел в ЭВМ и чувствительностью алгоритма к погрешностям округления. В связи со сказанным введем понятие устойчивости алгоритма. Алгоритм называется устойчивым, если в процессе его выполнения вычислительная погрешность возрастает незначительно, и неустойчивым – в противном случае. В качестве примера неустойчивого алгоритма рассмотрим следующий метод вычисления определенного интеграла . Интегрируя по частям, находим: , . Итак, для вычисления интеграла при всех натуральных получаем рекуррентную формулу Пользуясь этой формулой, находим: Очевидно, что значение интеграла не может быть отрицательным, поскольку подынтегральная функция неотрицательна на всем промежутке интегрирования . В чем же причина возникновения столь большой вычислительной погрешности? Погрешность при вычислении не превосходит 5 единиц последнего отброшенного разряда, т.е. . На каждом следующем шаге эта погрешность умножается на число, модуль которого больше . В итоге начальная погрешность, не превосходящая величины , умножается на , что и приводит к результату, не имеющему смысла. Таким образом, причиной возникновения столь большой вычислительной погрешности в данном примере является неустойчивость предложенного алгоритма к погрешностям округления. Отметим, что в данном примере на каждом шаге алгоритма вычисления выполнялись абсолютно точно, и только один операнд содержал наследственную погрешность, которая появилась в результате неточного вычисления числа . В общем случае все операнды содержат погрешности, т.е. являются приближенными числами, и все операции выполняются приближенно. Поэтому естественным представляется вопрос о том, как погрешности отдельных операндов сказываются на точности результата той или иной операции. Чтобы ответить на этот вопрос, введем понятия абсолютной и относительной погрешности приближенного числа. Пусть – точное, вообще говоря, неизвестное, значение некоторой величины, а – известное его приближение. Абсолютной погрешностью приближения называют такую величину , относительно которой известно, что . Таким образом, под абсолютной погрешностью приближенного значения величины понимается некоторая оценка модуля разности между приближенным и точным значениями этой величины. Относительной погрешностью приближения называется такая величина , относительно которой известно, что . Очевидно, что в качестве можно взять, например, . Сформулируем основную задачу теории погрешностей. Пусть известны погрешности некоторой совокупности величин; требуется определить погрешность заданной функции этих величин. Конкретнее, пусть – дифференцируемая функция своих аргументов; – известные приближения точных значений аргументов соответственно, причем абсолютные погрешности этих приближений заданы: . Требуется оценить абсолютную погрешность значения как приближения к точному значению . Эта задача имеет следующее решение: , (1) где , а , т.е. – некоторое промежуточное значение частной производной функции по аргументу в точке, расположенной на отрезке, соединяющем точную точку и приближенную точку . Так как точная точка и значение параметра нам неизвестны, то на практике пользуются приближенными оценками погрешности, следующими из формулы (1). Приведем конкретные оценки погрешностей для простейших функций. 1. . В этом случае . Следовательно, для абсолютной погрешности суммы и разности имеет место оценка . Таким образом, абсолютная погрешность алгебраической суммы приближенных слагаемых не превосходит суммы абсолютных погрешностей слагаемых. 2. . Легко показать, что в этом случае из (1) и определения относительной погрешности следует оценка . Следовательно, относительная погрешность произведения и частного приближенно равна сумме относительных погрешностей аргументов. 3. . Для этой функции имеем оценку , т.е. абсолютная погрешность синуса от приближенного значения аргумента не превосходит абсолютной погрешности аргумента. 4. . В этом случае имеет место приближенная оценка , т.е. относительная погрешность квадратного корня из приближенного значения не превосходит относительной погрешности аргумента. Лекции 2-4 Глава 1. Приближение функций § 1. Постановка задачи приближения функций Пусть некоторая величина является функцией аргумента . На практике явная связь между и зачастую неизвестна, т.е. невозможно записать эту связь в виде некоторой формульной зависимости . В некоторых случаях даже при известной зависимости она настолько сложна, например, содержит трудно вычисляемые выражения, специальные функции, сложные интегралы и т.п., что ее использование в практических расчетах затруднительно. Наиболее распространенным и практически важным случаем, когда зависимость между и неизвестна, является задание этой зависимости в виде таблицы , где – значения неизвестной функции в точках . Эти значения являются либо результатами предшествующих расчетов, либо экспериментальными данными. Пусть требуется найти значение величины в точках, отличных от точек , которые мы будем называть узлами. Этой цели и служит задача о приближении функций: неизвестную функцию требуется заменить некоторой функцией такой, что ее отклонение от функции является в некотором смысле наименьшим. Функция называется аппроксимирующей функцией. Одним из основных способов приближения функций является так называемое интерполирование. Оно состоит в следующем: аппроксимирующая функция подбирается из условия совпадения ее значений в узлах с заданными значениями : . На практике наибольшее распространение получил случай интерполирования алгебраическими многочленами, т.е. случай, когда функция выбирается в виде . Условия интерполяции в этом случае принимают вид: . (1.1) Условия (1.1) представляют собой систему из линейных алгебраических уравнений относительно неизвестных . Мы знаем, что в случае, когда число уравнений больше числа неизвестных, система, как правило, несовместна (она является переопределенной). В случае, когда число уравнений меньше числа неизвестных, система является неопределенной. Обычно такие системы имеют бесчисленное множество решений. Поэтому естественно рассмотреть случай, когда число уравнений в системе (1.1) совпадает с числом неизвестных, т.е. случай, когда , или . В этом случае аппроксимирующая функция является многочленом степени , т. е. а система (1) принимает вид . (1.2) Для выяснения вопроса о разрешимости системы (1.2) выпишем ее основной определитель: . Это известный из курса алгебры определитель Вандермонда, он равен . Ясно, что тогда и только тогда, когда все узлы различны. Таким образом, если все узлы интерполяции попарно различны, то система (1.2) имеет и притом единственное решение. Тем самым доказаны существование и единственность интерполяционного многочлена степени , совпадающего в заданных точках с заданными значениями . Поскольку многочлен используется для приближения функции на всем рассматриваемом интервале изменения переменной , то в этом случае говорят о глобальной интерполяции. Если количество узлов интерполяции велико, то глобальная интерполяция приводит к многочленам высокой степени, что не совсем удобно с вычислительной точки зрения, так как вычисление значений многочленов высокой степени сопровождается большой вычислительной погрешностью. Поэтому часто интерполяционные многочлены строят отдельно для различных небольших частей рассматриваемого интервала изменения аргумента . В этом случае говорят о кусочной, или локальной, интерполяции. В любом случае при интерполировании основным условием является прохождение графика аппроксимирующей функции через данные табличные точки. Однако в ряде случаев выполнение этого условия затруднительно или даже нецелесообразно. Допустим, например, что табличные значения получены путем измерений и, следовательно, содержат погрешности. Построение аппроксимирующей функции с условием обязательного прохождения ее графика через экспериментальные точки означало бы тщательное повторение допущенных при измерениях погрешностей. В подобных случаях часто используется другой подход к построению аппроксимирующей функции, называемый методом наименьших квадратов. Мы рассмотрим идею этого метода для случая, когда в качестве аппроксимирующей функции выбирается многочлен степени . Мерой отклонения многочлена от функции на множестве точек назовем величину , равную сумме квадратов разностей между значениями аппроксимирующего многочлена и приближаемой функции в узлах. Согласно методу наименьших квадратов коэффициенты многочлена подбираются так, чтобы функция принимала наименьшее возможное значение. На практике обычно ограничиваются многочленами небольшой степени, чаще всего 1, 2 или 3. Упражнение. Методом наименьших квадратов построить аппроксимирующие многочлены первой и второй степени для функции , заданной таблицей , где – значения функции в узлах . § 2. Интерполяционный многочлен Лагранжа Мы доказали существование и единственность интерполяционного многочлена степени , (1.3) совпадающего с заданными значениями в узлах интерполяции . Его коэффициенты являются решением системы . На практике нет необходимости решать эту систему, тем более что ее решение является очень чувствительным к погрешностям округления. Как мы сейчас увидим, для интерполяционного многочлена (1.3) можно выписать ряд явных представлений. Рассмотрим многочленов , каждый из которых имеет степень и удовлетворяет следующим условиям: (1.4) Условия (1.4) однозначно определяют многочлены . Действительно, так как при , то точки являются корнями многочлена . По теореме Безу это означает, что многочлен делится на двучлены , т.е. его можно представить в виде . Заметим, что степень многочлена по определению должна быть равна . Так как степень произведения , очевидно, уже равна , то . Значение этой постоянной определим из условия . Имеем: . Тогда . Покажем теперь, что многочлен (1.5) является искомым интерполяционным многочленом. Действительно, он имеет степень не выше чем , так как является линейной комбинацией многочленов степени . Далее для всех имеем , т.е. многочлен совпадает с заданными значениями функции в узлах интерполяции. Полученное выражение (1.5) интерполяционного многочлена называется интерполяционным многочленом Лагранжа. В дальнейшем нам потребуется оценка остаточного члена интерполяционного многочлена. Пусть функция раз непрерывно дифференцируема. Тогда имеет место следующее представление: . Здесь , а – некоторая промежуточная точка, расположенная между точками и . § 3. Разделенные разности. Интерполяционный многочлен Ньютона с разделенными разностями Пусть дана совокупность узлов . Значения функции в этих узлах будем называть разделенными разностями нулевого порядка. Разделенными разностями первого порядка назовем величины . Разделенными разностями второго порядка называются величины . Разделенными разностями -го порядка назовем величины . Таким образом, разделенные разности более высокого порядка определяются рекуррентно через разделенные разности порядка, меньшего на . Вместе с тем можно выписать выражение разделенной разности любого порядка непосредственно через значения функции в узлах, т.е. через разности нулевого порядка. А именно, имеет место формула . (1.6) Поскольку значение суммы не зависит от порядка слагаемых, то из формулы (1) следует, что разделенная разность не изменяется при перестановке ее аргументов, например, . Все разделенные разности, построенные по совокупности узлов , принято располагать в виде следующей треугольной таблицы: Перейдем теперь к другой форме записи интерполяционного многочлена, использующей разделенные разности. Преобразуем остаточный член многочлена Лагранжа: = . Согласно формуле (1.6) выражение в скобках есть разделенная разность , построенная по узлам . Таким образом, . (1.7) Обозначим через интерполяционный многочлен степени , построенный по узлам . Представим многочлен в виде (1.8) и рассмотрим разность . Она представляет собой многочлен степени , обращающийся в в точках . Действительно, по определению интерполяционного многочлена , . Следовательно, при имеем: . Отсюда по теореме Безу следует, что многочлен делится на двучлены . Значит, многочлен можно записать в виде . Заметим, что степень многочлена равна , т.е. совпадает со степенью многочлена . Поэтому . Для определения этой постоянной вычислим значение многочлена в узле : . Но , поэтому полученное равенство можно переписать в виде . (1.9) С другой стороны, воспользуемся равенством (1.7) при и : , или . (1.10) Сравнивая равенства (1.9) и (1.10), получаем: . Следовательно, . Замечая, что , и подставляя полученные выражения в (1.8), приходим к интерполяционному многочлену Ньютона с разделенными разностями: . В заключение этого параграфа сравним формулу (1.7) и представление остаточного члена многочлена Лагранжа, приведенное в конце предыдущего параграфа. В результате приходим к формуле , где . Таким образом, оказывается, что разделенная разность порядка лишь числовым множителем отличается от значения производной порядка в некоторой промежуточной точке. Отсюда, в частности, следует, что если функция является многочленом степени , то разделенные разности порядка постоянны, а все разделенные разности более высоких порядков равны нулю. § 4. Интерполирование с кратными узлами Иногда приходится решать более общую задачу интерполирования, когда известны не только значения некоторой функции в узлах , но и значения ее производных до некоторого порядка. Строго эта задача формулируется следующим образом. Пусть требуется построить многочлен степени , удовлетворяющий условиям (1.11) где все узлы попарно различны и . (1.12) Числа называются кратностями узлов соответственно. Можно показать, что поставленная задача имеет и притом единственное решение. Это решение можно записать в следующем виде: . Выписанное выражение называется интерполяционным многочленом Эрмита. При его записи мы использовали так называемые разделенные разности с кратными узлами, например, . Ясно, что обычное определение разделенных разностей для вычисления не подходит. Разделенные разности с кратными узлами определяются следующим образом. Если все аргументы разделенной разности одинаковы, то следует воспользоваться формулой . Если же среди аргументов есть различные, то разделенную разность можно свести к разделенным разностям более низкого порядка с помощью таких же формул, что и в случае, когда все узлы различны: = . В частности, , , . Можно показать, что для остаточного члена многочлена Эрмита имеет место представление, аналогичное представлению остаточного члена многочлена Ньютона: , где . Имеет место и представление, аналогичное представлению остаточного члена многочлена Лагранжа: , где . Пример. Построить многочлен Эрмита , удовлетворяющий условиям Решение. Здесь . Имеем: . Так как , то окончательно получаем: . § 5. Конечные разности. Интерполяционные формулы для случая равноотстоящих узлов Пусть узлы расположены друг от друга на равном расстоянии : . Число называется шагом таблицы узлов. Соответствующие этим узлам значения функции будем обозначать через и называть конечными разностями нулевого порядка. Конечными разностями первого порядка называются величины , которые в зависимости от обстановки обозначаются и именуются по-разному: – центральная разность; – разность вперед; – разность назад. В общем случае конечные разности -го порядка определяются через разности порядка следующим образом: ; ; . Конечные разности принято располагать в виде треугольной таблицы: Имеют место следующие свойства конечных разностей. 1. Конечные разности -го порядка выражаются через значения функции в узлах по формуле , где при четном число – целое, а при нечетном число – полуцелое. В частности, , , . 2. Если , то имеет место следующая связь между разделенными и конечными разностями: . (1.13) Воспользуемся указанной в конце § 3 связью между разделенной разностью и производной , (1.14) где . Сравнивая (1) и (2), получаем: . Следовательно, конечная разность порядка лишь числовым множителем отличается от значения производной порядка в некоторой промежуточной точке. Отсюда, в свою очередь, следует, что если – многочлен степени , то конечные разности порядка и выше будут равны нулю. Перейдем к построению интерполяционных формул, использующих конечные разности. Если данная таблица состоит из большого количества узлов, нам нет необходимости использовать их все. Это привело бы к увеличению степени интерполяционного многочлена, а значит, и к увеличению объема вычислений. На практике чаще всего достаточно использовать интерполяционный многочлен невысокой степени. Однако с целью достижения приемлемой точности интерполирования желательно выбирать те узлы таблицы, которые ближе всего расположены к точке, в которой выполняется интерполирование. Действительно, в выражении остаточного члена фигурирует величина , равная произведению расстояний от точки , в которой выполняется интерполирование, до узлов таблицы. Очевидно, что величина будет минимальной, если изо всех имеющихся узлов использовать ближайших к точке . Если точка находится внутри таблицы узлов, то выгодно взять половину узлов слева от и половину – справа от . Если же точка находится в начале или в конце таблицы, то такой совокупности узлов выбрать нельзя. Сейчас мы рассмотрим как раз такие случаи. Пусть точка расположена в начале таблицы узлов . Построим интерполяционный многочлен Ньютона с разделенными разностями по этой совокупности узлов: . Сделаем замену переменной и учтем, что . Тогда . Переходя в выражении для от разделенных разностей к конечным, получаем: . Полученная формула называется первой интерполяционной формулой Ньютона, или формулой Ньютона для интерполирования вперед. Рассмотрим теперь случай, когда точка находится в конце таблицы узлов . В этом случае удобнее привлекать узлы в следующей последовательности: . Построим по этой совокупности интерполяционный многочлен Ньютона с разделенными разностями: . Сделаем замену переменной и учтем, что . Тогда . Перейдем от разделенных разностей к конечным: , , , . Подставляя эти формулы в выражение для , получаем вторую интерполяционную формулу Ньютона, или формулу Ньютона для интерполирования назад: . Пусть теперь точка расположена внутри таблицы узлов. Обозначим через узел, ближайший к точке . Если , то оптимальным набором узлов будет совокупность . Построим по этой совокупности интерполяционный многочлен Ньютона с разделенными разностями: . Сделаем замену переменной и учтем, что . Тогда . Воспользуемся связью между разделенными и конечными разностями: , , , , , . Следовательно, . Полученная формула называется первой интерполяционной формулой Гаусса, или формулой Гаусса для интерполирования вперед. Если же , то оптимальной будет совокупность узлов . В этом случае получаем формулу , которая называется второй интерполяционной формулой Гаусса, или формулой Гаусса для интерполирования назад. § 6. Интерполирование с помощью кубических сплайнов Пусть функция определена на отрезке . На этом отрезке рассмотрим сетку , т.е. множество узлов . Функция называется интерполяционным кубическим сплайном для функции относительно сетки , если она удовлетворяет трем условиям: 1) , т.е. непрерывна на отрезке вместе со своими производными до второго порядка включительно; 2) на каждом частичном отрезке она является многочленом третьей степени, который мы обозначим через ; 3) . Замечание. Условия 1)-2) определяют кубический сплайн; условие 3) есть условие интерполирования сплайном функции в узлах сетки . Условие 2) означает, что при . Многочлен полностью определяется четырьмя коэффициентами . Поскольку общее число промежутков равно , то для полного определения сплайна нужно найти значения неизвестных . Пользуясь определением сплайна, выпишем уравнения, из которых определяются эти коэффициенты. Многочлен является функцией, бесконечное число раз непрерывно дифференцируемой, поэтому условие 1) сводится к требованию непрерывности самого сплайна и его первой и второй производных и во внутренних узлах , т.е. в точках, являющихся общими для соседних промежутков и : , (1.15) , (1.16) . (1.17) В силу условия 3) определения сплайна значения и не только совпадают между собой, но и равны значению . Поэтому, заменив соотношения (1.15) равенствами , (1.18) , (1.19) мы учтем как условие непрерывности сплайна во внутренних узлах , так и условие интерполяционности сплайна в этих узлах. К соотношениям (1.18) добавим условие интерполяционности сплайна в узле , а к соотношениям (1.19) – условие интерполяционности в узле . Соберем вместе все полученные уравнения: , (1.20) , (1.21) , (1.22) . (1.23) Соотношения (1.20)-(1.23) представляют собой систему из линейных уравнений относительно неизвестных. Для замыкания этой системы требуется еще два уравнения. Для их получения привлекают некоторые дополнительные условия, обычно индуцируемые физической постановкой задачи и называемые краевыми условиями. Рассмотрим наиболее употребительные краевые условия. а) , где и – известные числа; б) , где и – известные числа; в) , где и – кубические многочлены, графики которых проходят соответственно через первые четыре и последние четыре узла сетки; г) если известно, что приближаемая функция является периодической с периодом , то и сплайн также должен быть периодическим; в этом случае краевые условия задают в виде . В любом из перечисленных случаев мы получаем два дополнительных уравнения, необходимых для замыкания системы (1.20)-(1.23). На практике нет необходимости в явном виде составлять и затем решать полученную систему. Сейчас мы покажем, что сплайн полностью определяется значениями своей первой или второй производной в узлах сетки, что позволяет вместо системы из уравнений относительно неизвестных решать систему только из уравнений с таким же числом неизвестных. Введем обозначения: , . Числа называются наклонами сплайна , а числа – его моментами. Покажем, что сплайн целиком определяется, например, своими моментами . Для этого достаточно показать, что все коэффициенты можно выразить через моменты , а также через данные значения и . На отрезке сплайн имеет вид . Коэффициенты найдем из условий: , (1.24) . (1.25) Имеем: . Введем обозначение . Тогда из (1.24) находим , откуда . Из (1.25) находим , откуда . В результате получаем следующее выражение для через моменты: . Обычно выражение для преобразуют и используют в следующем виде: . В таком виде выражение для многочлена легче запоминается и может быть эффективно использовано для целей интерполирования функции , численного дифференцирования и интегрирования, если, конечно, предварительно определены моменты . Для определения моментов воспользуемся пока еще не использованным условием непрерывности первой производной сплайна во внутренних узлах и краевыми условиями. Имеем: . Аналогично . Запишем теперь условия (1.20): , или . После умножения полученного равенства на оно принимает вид . Так как условия (1.20) должны быть выполнены для всех , то мы получили систему из уравнений для определения неизвестных . Недостающие два уравнения дает использование краевых условий. Наиболее простые уравнения получаются при использовании краевых условий типа б). В этом случае они имеют вид , где числа и заданы. В случае использования краевых условий типа а) дополнительные уравнения получаются немного сложнее. Имеем: , . В результате элементарных преобразований получаем: , . После умножения первого из этих равенств на , а второго – на они принимают вид , . Аналогичным образом получают дополнительные уравнения и при использовании других краевых условий. Итак, в случае использования краевых условий типа б) для определения моментов приходим к системе (1.26) а в случае краевых условий типа а) получаем систему (1.27) И система (1.26) и система (1.27) являются системами линейных алгебраических уравнений с трехдиагональной матрицей вида , причем ее элементы удовлетворяют условиям: , . Последние условия означают, что матрица обладает свойством диагонального преобладания. Известно, что такая матрица невырождена. Следовательно, решения обеих систем (1.26) и (1.27) существуют и определяются единственным образом. § 7. Метод прогонки решения системы с трехдиагональной матрицей Опишем эффективный метод решения системы с трехдиагональной матрицей, имеющей диагональное преобладание, – метод прогонки. Рассмотрим систему (1.28) Из первого уравнения выразим : . Введем обозначения: . Тогда запишется в виде . Предположим, что , и подставим это выражение в -ое уравнение системы: . Отсюда выразим : . Обозначим: . Тогда . Тем самым мы доказали, что для всех справедливо представление , и получили рекуррентные формулы для так называемых прогоночных коэффициентов , : (1.29) При имеем: . Подставим это выражение в последнее уравнение системы: . Отсюда найдем : . (1.30) Затем для находим по формулам . (1.31) Вычисление прогоночных коэффициентов , по формулам (1.29) называется прямым ходом метода прогонки; вычисление неизвестных системы (1.28) по формулам (1.30)-(1.31) называется обратным ходом метода прогонки. Известно, что выполнение условий диагонального преобладания является не только достаточным для осуществления метода прогонки, но и обеспечивает вычислительную устойчивость этого метода, заключающуюся в том, что погрешности, допускаемые на отдельных шагах метода, не увеличиваются на следующих шагах. Лекция 5 Глава 2. Численное дифференцирование Простейшие формулы численного дифференцирования, т.е. формулы для приближенного вычисления значений производных таблично заданных функций, получаются дифференцированием интерполяционных формул. Пусть функция задана таблично, т.е. известны ее значения в узлах , и – ее интерполяционный многочлен. Положим . (2.1) Получим представление остаточного члена построенной формулы (2.1). Имеем: , где – остаточный член интерполяционной формулы . По формуле Ньютона-Лейбница получаем: . Обозначим: . По определению разделенной разности с кратными узлами имеем: . Используя формулу связи между разделенной разностью и производной, находим: , где – некоторая точка на промежутке , , . Окончательно получаем следующее представление остаточного члена формулы численного дифференцирования (2.1): . (2.2) Отсюда легко следует оценка остаточного члена формулы численного дифференцирования (2.1): . (2.3) На практике особый интерес представляют формулы численного дифференцирования для функций, заданных таблично на равномерной сетке, т.е. на сетке, узлы которой являются равноотстоящими друг от друга. В качестве примера получим приближенную формулу для вычисления значения первой производной с погрешностью на основе первой интерполяционной формулы Ньютона. Здесь – шаг сетки. Выпишем первую интерполяционную формулу Ньютона для узлов , где , : . Соответствующая ей формула численного дифференцирования есть . Остаточный член этой формулы согласно (2.2) имеет вид . Имеем: , . Следовательно, . Поскольку мы поставили задачу построения формулы численного дифференцирования, имеющей погрешность , то следует принять . В таком случае мы должны дифференцировать многочлен . Имеем: . Так как , то , и, следовательно, . Это и есть искомая формула для приближенного вычисления первой производной с погрешностью . Обычно нас интересуют приближенные значения производной таблично заданной функции в узлах таблицы. Например, пусть . Тогда , и мы получаем: . Итак, с погрешностью имеем формулу численного дифференцирования . Совершенно аналогичным образом на основе различных интерполяционных формул с конечными разностями можно получить другие формулы численного дифференцирования с различной величиной погрешности. Приведем без вывода несколько таких формул: с погрешностью ; с погрешностью ; с погрешностью ; с погрешностью . Лекции 6-7 Глава 3. Численное интегрирование § 1. Квадратуры Ньютона-Котеса Пусть требуется вычислить определенный интеграл , (3.1) где – некоторая интегрируемая функция. Мы рассмотрим ряд методов приближенного вычисления интеграла , в которых интеграл заменяется линейной комбинацией значений функции в конечном числе точек отрезка интегрирования . Другими словами, мы будем рассматривать методы, в которых . (3.2) Формулы типа (3.2) называются квадратурными, числа – коэффициентами квадратуры (3.2), а точки – узлами квадратуры. Простейшие квадратурные формулы получаются при помощи аппарата интерполирования. На отрезке выберем узлов . С помощью замены переменной (3.3) отрезок переходит в отрезок , а точки – в некоторые точки . Таким образом, выбор узлов на отрезке индуцирует некоторый выбор узлов на отрезке . По узлам построим интерполяционный многочлен для функции . В общем случае некоторые из узлов будем считать кратными. Напомним, что узел называется кратным, если в нем требуется совпадение значений не только самой функции, но и ее производных до некоторого порядка . При этом число называется кратностью узла . Как и ранее, через обозначим суммарную кратность всех узлов, т.е. . Итак, в общем случае будем строить интерполяционный многочлен Эрмита степени . Положим теперь . Получаемые при этом формулы называются квадратурными формулами Ньютона-Котеса. Рассмотрим остаточный член квадратуры Ньютона-Котеса: . Здесь – остаточный член интерполяционного многочлена Эрмита. Воспользуемся известным из теории интерполирования представлением , где . Так как все узлы принадлежат отрезку , то для всех справедлива оценка , где . Следовательно, . Воспользуемся заменой (3.3). Имеем: , . Здесь . Тогда . Обозначим: . Тогда окончательно получаем следующую оценку: . (3.4) Рассмотрим теперь частный случай, когда все узлы являются простыми. В этом случае ; . Снова воспользуемся заменой переменной (3.3). Имеем: , . Следовательно, . В результате получаем следующую квадратуру , (3.5) где . (3.6) Заметим, что поскольку при построении квадратуры (3.5) все узлы являлись простыми, то для нее оценка (3.4) остаточного члена принимает вид . Из этой оценки следует, что квадратура (3.5) точна в случае, когда функция является произвольным многочленом степени , ибо для такой функции производная порядка тождественно равна нулю. § 2. Симметричные квадратуры Ньютона-Котеса На практике узлы квадратурной формулы обычно выбирают так, чтобы они располагались симметрично относительно середины отрезка интегрирования : . Для соответствующих им узлов на отрезке это означает, что . В дальнейшем будем рассматривать квадратуры именно с таким расположением узлов. Такие квадратуры называются симметричными. По сравнению с общим случаем симметричные квадратуры Ньютона-Котеса (3.5), в которых коэффициенты определяются по формулам (3.6), обладают рядом дополнительных свойств. Рассмотрим эти свойства. 1. Коэффициенты квадратуры, соответствующие симметрично расположенным узлам, равны: . Доказательство. В интеграле , определяющем коэффициент , сделаем замену переменной . Тогда , , . 2. Симметричная квадратура точна для любой нечетной относительно середины отрезка интегрирования функции. Доказательство. Пусть функция нечетна относительно середины отрезка : . В этом случае . С другой стороны, . Если в этой сумме число членов четно, то все они взаимно уничтожаются. Если же нечетно, то центральное слагаемое остается, но оно равно нулю в силу нечетности функции относительно точки . Таким образом, в любом случае имеем: . Следовательно, , а это и означает, что квадратура точна. 3. Симметричная квадратура с нечетным числом узлов точна для всех многочленов степени . Доказательство. Рассмотрим произвольный многочлен степени . Очевидно, что его можно представить в виде , где – некоторый многочлен степени . Покажем, что первое слагаемое – функция, нечетная относительно точки : . В силу свойства 2 симметричная квадратура точна для такой функции, т.е. . Так как – многочлен степени , то для него будет точна любая квадратура Ньютона-Котеса, построенная по узлам. Следовательно, и . Остается показать, что остаточный член квадратуры Ньютона-Котеса является аддитивной функцией, т.е. для любых функций и имеет место равенство . Имеем: . Отсюда следует, что . Тем самым свойство 3 доказано. Естественно ожидать, что оценка остаточного члена симметричной квадратуры при нечетном должна быть более тонкой, чем в общем случае. Именно с целью получения более тонких оценок остаточного члена в данном и ему подобных случаях мы и допустили с самого начала, что некоторые узлы могут быть кратными. Конкретнее, в случае, когда квадратура симметрична, а число узлов нечетно, мы будем использовать центральный узел как двукратный, а все остальные узлы – как простые. При этом оказывается, что квадратура, получающаяся при замене функции многочленом Эрмита , совпадает с квадратурой, построенной только по простым узлам, т.е. получающейся при замене функции многочленом Лагранжа . Действительно, пусть узел – двукратный, а остальные узлы – простые. Построим многочлен Эрмита по совокупности узлов : . Покажем, что функция в силу нечетности числа и симметричного расположения узлов нечетна относительно средней точки : . Отсюда следует, что . Поэтому . Таким образом, мы показали, что если число узлов нечетно и квадратура симметрична, то она может быть получена двумя способами: 1) заменой функции многочленом Лагранжа и 2) заменой функции многочленом Эрмита . Соответствующие этим способам получения квадратуры оценки остаточного члена имеют вид: , (3.7) . (3.8) Естественно из двух оценок, справедливых для одной и той же формулы, пользоваться более тонкой оценкой (3.8). § 3. Употребительные квадратуры Ньютона-Котеса 1. Формула левосторонних прямоугольников. Положим , , . Следовательно, . Имеем: . Поскольку , то формула левосторонних прямоугольников имеет вид . Согласно (3.7) оценка остаточного члена этой формулы имеет вид . 2. Формула правосторонних прямоугольников. В этой формуле , , , . Имеем: . Поскольку , то формула правосторонних прямоугольников имеет вид . Согласно (3.7) оценка остаточного члена этой формулы имеет вид . 3. Формула центральных прямоугольников. В этой формуле , , , . Имеем: . Так как , то формула принимает вид . В силу (3.8) оценка остаточного члена этой формулы имеет вид . 4. Формула трапеций. В этой формуле , , , , , . Имеем: , . Так как , мы получаем формулу , оценка остаточного члена которой согласно (3.7) имеет вид: . 5. Формула Симпсона. В формуле Симпсона , , , , , , , . Тогда , , . Мы получаем формулу с оценкой остаточного члена . Заметим, что при малой длине отрезка интегрирования эта оценка на два порядка лучше оценок формул центральных прямоугольников и трапеций. § 4. Составные квадратуры Ньютона-Котеса До сих пор, сравнивая квадратурные формулы по порядку погрешности, мы предполагали, что длина отрезка интегрирования мала (меньше 1). Если же это не так, то рассмотренные формулы применяют к отдельным малым по длине частям этого отрезка. Получаемые при этом формулы называются составными квадратурными формулами. Разобьем отрезок на частей точками и для вычисления интеграла по каждой из частей применим одну из квадратурных формул Ньютона-Котеса: . Здесь . Пусть – кратность узла . Тогда оценка остаточного члена полученной составной формулы имеет вид: . Заметим, что, как и в простых формулах, здесь при четном все , а при нечетном и симметричном расположении узлов равны 1 все , кроме одного: . В простейшем и наиболее употребительном случае точки выбираются равноотстоящими: . Тогда полученная формула принимает вид , а оценка остаточного члена – вид: . Приведем конкретные примеры составных квадратурных формул. 1. Составная формула левосторонних прямоугольников: . Оценка остаточного члена составной формулы левосторонних прямоугольников имеет вид . 2. Составная формула правосторонних прямоугольников: . Оценка остаточного члена составной формулы правосторонних прямоугольников совпадает с оценкой погрешности составной формулы левосторонних прямоугольников: . 3. Составная формула центральных прямоугольников: . Здесь . Оценка остаточного члена составной формулы центральных прямоугольников имеет вид . 4. Составная формула трапеций: . Оценка ее остаточного члена: . 5. Составная формула Симпсона: . Оценка ее остаточного члена: . Замечание. Обычно формулу Симпсона записывают в терминах расстояний между соседними узлами квадратуры : . Оценка остаточного члена принимает в этом случае вид: . § 5. Правило Рунге Обозначим: . Зададим некоторое целое и будем последовательно вычислять при . После каждого вычисления найдем величину , где – суммарная кратность узлов используемой простой формулы, и проверим, выполнено ли неравенство , где – требуемая точность вычисления интеграла . Если при некотором это неравенство выполняется, то в качестве приближенного значения с требуемой степенью точности принимаем или, что, как правило, лучше, . Описанное правило достижения требуемой точности вычисления определенного интеграла называется правилом Рунге. § 5. Численное интегрирование с помощью кубических сплайнов Рассмотрим другой способ приближенного вычисления определенного интеграла . На отрезке рассмотрим сетку . Пусть сетка является равномерной, т.е. ее узлы находятся на равном расстоянии друг от друга: . Для функции построим интерполяционный кубический сплайн относительно сетки . Заменяя подынтегральную функцию ее сплайном , получим приближенную формулу . (3.9) Имеем: . Здесь – кубический многочлен, представляющий сплайн на частичном отрезке . Напомним, что интерполяционный кубический сплайн однозначно определяется своими наклонами или своими моментами . Подробно рассмотрим случай, когда сплайн задается моментами. В этом случае кубический многочлен , представляющий сплайн на частичном отрезке , имеет вид . Подставим это выражение в интеграл : . Тогда . Окончательно формула (3.9) принимает вид . Лекции 6-7 Глава 5. Прямые методы решения систем линейных алгебраических уравнений § 1. Прямые методы решения СЛАУ. Метод Гаусса Мы рассмотрим методы решения СЛАУ двух типов – прямые и итерационные. В прямых методах точное решение системы при отсутствии погрешностей округления находится за конечное число арифметических действий. В итерационных методах определяется некоторая последовательность приближений, которая в пределе стремится к точному решению системы. Одним из основных прямых методов решения СЛАУ является метод Гаусса. Под методом Гаусса понимается целая совокупность методов, в основе которых лежит использование элементарных преобразований матрицы системы с целью приведения ее к треугольному виду. Мы рассмотрим одну из возможных схем этого метода, в результате применения которой матрица системы принимает верхнюю треугольную форму. Рассмотрим систему линейных уравнений относительно неизвестных . Предположим, что все угловые миноры матрицы отличны от нуля: . На первом шаге матрица преобразуется к матрице , все элементы первого столбца которой, начиная со второго, равны нулю. Согласно предположению угловой минор первого отличен от нуля. Поэтому мы можем определить множители первого шага . Для всех вычтем из -ой строки матрицы первую, умноженную на . В результате получим матрицу требуемого вида: . При этом значения элементов находятся по формулам . Элемент , на который приходится делить при определении множителей первого шага, называется ведущим элементом первого шага. Отметим для дальнейшего, что описанные преобразования матрицы не изменяют значений ее угловых миноров. Поэтому все угловые миноры полученной матрицы отличны от нуля. В частности, отличен от нуля угловой минор второго порядка этой матрицы. Но он равен произведению . Следовательно, отличен от нуля элемент , который мы назовем ведущим элементом второго шага. Целью второго шага является получение нулей во втором столбце матрицы , начиная с третьего элемента этого столбца. Определим множители второго шага , и для всех из -ой строки матрицы вычтем вторую, умноженную на . Ясно, что описанные преобразования не изменяют первых двух строк и первого столбца матрицы . В результате приходим к матрице , где Пусть в результате выполнения шагов мы пришли к матрице . Так как преобразования, выполняемые на каждом из шагов, не меняют значений угловых миноров, то все угловые миноры матрицы совпадают с соответствующими угловыми минорами исходной матрицы . В частности, отличен от нуля угловой минор -ого порядка матрицы . Но он равен произведению . Следовательно, . Назовем элемент ведущим элементом -ого шага и перейдем к описанию этого шага. Его цель – получение нулей в -ом столбце матрицы , начиная с -ого элемента. Определим множители -ого шага , и для всех из -ой строки матрицы вычтем ее -ую строку, умноженную на . Так как при этом первые строк и первые столбцов матрицы остаются неизменными, то мы приходим к матрице вида , где . Ясно, что при мы придем к матрице , имеющей верхнюю треугольную форму. В силу невырожденности исходной матрицы матрица также невырождена. Описанный процесс приведения матрицы системы к верхней треугольной форме называется прямым ходом метода Гаусса. Отметим, что мы описали процесс преобразования только матрицы системы . Для того чтобы получающаяся в результате система оставалась эквивалентной исходной системе, следует, очевидно, с компонентами правой части выполнить те же самые преобразования, что и со строками матрицы . Так, на первом шаге для всех из следует вычесть первую компоненту , умноженную на соответствующий множитель . В результате мы придем к вектору , где . Точно таким же образом поступаем на всех следующих шагах. Именно, пусть в результате выполнения шагов мы пришли к вектору . На -ом шаге для всех из компоненты вычтем компоненту , умноженную на . В результате получим вектор , где . При придем к окончательно преобразованной правой части . При этом система преобразуется к эквивалентной ей системе . В дальнейшем матрицу будем обозначать через , а вектор – через . Так как матрица системы является верхней треугольной, то решение полученной системы находится очень просто. Последнее уравнение системы имеет вид , причем . Отсюда находим . Пусть значения неизвестных уже определены. Подставляя их в -ое уравнение системы и учитывая, что , находим . Описанный процесс решения системы с верхней треугольной матрицей называется обратным ходом метода Гаусса. § 2. -разложение матрицы При выполнении прямого хода метода Гаусса мы выполняли элементарные преобразования со строками матрицы . Выполнение каждой такой операции равносильно умножению матрицы слева на некоторую невырожденную матрицу. Действительно, выполнение первого шага, т.е. вычитание из -ой строки матрицы ее первой строки, умноженной на , проделанное для всех , равносильно умножению матрицы слева на матрицу . Таким образом, . Преобразования второго шага равносильны умножению матрицы слева на матрицу . Поэтому . По аналогии легко выписать матрицу , осуществляющую преобразования, выполняемые на -ом шаге: . При этом . В результате выполнения шагов мы приходим к невырожденной верхней треугольной матрице . Легко выписать выражение матрицы через исходную матрицу : . Обозначив произведение через , получим . В силу невырожденности множителей матрица невырождена. Обозначим через матрицу, обратную к . Оказывается, что матрица имеет вид . Доказательство этого утверждения состоит в непосредственной проверке равенства , где – единичная матрица. Представление матрицы в виде произведения нижней треугольной матрицы с единицами на главной диагонали и невырожденной верхней треугольной матрицы называется ее -разложением. -разложение играет важную роль в численных методах линейной алгебры. Так, если для матрицы системы получено -разложение , то решение этой системы сводится к последовательному решению двух систем с треугольными матрицами . При наличии -разложения матрицы особенно выгодным является решение большого количества систем с одной и той же матрицей и различными правыми частями. Дело в том, что по количеству необходимых арифметических действий вычисление -разложения является более трудоемкой операцией, чем решение систем с треугольными матрицами. Конкретнее, вычисление -разложения требует арифметических действий, а решение системы с треугольной матрицей – . К решению большого числа систем с одной и той же матрицей сводится, например, задача обращения матрицы. Пусть требуется найти матрицу , обратную к , т.е. решить матричное уравнение . Обозначив через -ый столбец матрицы , матричное уравнение запишем в виде совокупности систем уравнений . Таким образом, для вычисления обратной матрицы требуется решить систем уравнений с одной и той же матрицей . В заключение параграфа упомянем одно достоинство -разложения, связанное с возможностью экономичного использования машинной памяти в процессе его вычисления на ЭВМ. Вспомним, что на первом шаге мы получаем нули в первом столбце, начиная со второго элемента этого столбца. Этих нулей ровно столько, сколько имеется множителей первого шага, определяющих первый столбец матрицы . Поскольку хранить нули в памяти нет необходимости, то на их местах можно как раз разместить множители первого шага. Заметим также, что первая строка матрицы , которая вообще не изменяется, является в то же время и первой строкой результирующей матрицы . Таким образом, в результате выполнения первого шага первая строка и первый столбец двумерного массива , в котором первоначально размещалась матрица системы, содержат полную информацию о первой строке матрицы и первом столбце матрицы . Остальные элементы этого массива суть элементы матрицы . На втором шаге, оставляя первые две строки неизменными, мы преобразуем остальные строки, получая нули во втором столбце ниже главной диагонали. Этих нулей ровно столько, сколько множителей определяют второй столбец матрицы . Поэтому мы можем разместить множители на местах этих нулей. При этом вторая строка, начиная со второго элемента, полностью определяет вторую строку матрицы . Продолжая действовать подобным образом, мы сможем разместить полную информацию о сомножителях и в одном двумерном массиве , так что при вычислении -разложения не требуется никакой дополнительной памяти. § 3. Метод квадратного корня Предположим, что все угловые миноры матрицы отличны от нуля. Такая матрица имеет -разложение , где – нижняя треугольная матрица с единицами на главной диагонали, а – невырожденная верхняя треугольная матрица. Невырожденность матрицы означает, что ее диагональные элементы отличны от нуля. В таком случае матрицу можно представить в виде , где , , а , . Нетрудно проверить, что данное представление определяется единственным образом. -разложение матрицы приобретает теперь вид . (4.1) Представление матрицы в виде произведения нижней треугольной матрицы с единицами на главной диагонали, невырожденной диагональной матрицы и верхней треугольной матрицы с единицами на главной диагонали называется ее -разложением. Таким образом, равенство (4.1) является -разложением матрицы . Легко видеть, что -разложение определяется единственным образом. Рассмотрим -разложение симметричной матрицы . Из равенства следует , или . Так как – нижняя треугольная матрица с единицами на главной диагонали, а – верхняя треугольная матрица с единицами на главной диагонали, то представление также является -разложением матрицы . В силе единственности -разложения , так что . (4.2) Предположим теперь, что матрица не только симметрична, но и положительно определена. Согласно критерию Сильвестра все угловые миноры матрицы положительны. Поэтому для симметричной положительно определенной матрицы -разложение вида (4.2) заведомо существует. Покажем, что матрица в этом разложении также положительно определена. Воспользуемся невырожденностью матрицы и в выражении сделаем замену переменной : . Очевидно, что из условия вытекает и условие . Но матрица предполагается положительно определенной, следовательно, , а это и означает, что матрица также положительно определена. В силу критерия Сильвестра все угловые миноры матрицы положительны, откуда следует, что при всех . Поэтому мы можем рассмотреть матрицу , которая, очевидно, удовлетворяет равенству . Введем в рассмотрение еще одну матрицу , которая, очевидно, является невырожденной верхней треугольной матрицей. С использованием этого обозначения разложение (4.2) принимает вид . (4.3) Это разложение приводит к одному из наиболее употребительных методов решения систем линейных уравнений с симметричными положительно определенными матрицами – методу квадратного корня. Метод квадратного корня заключается в получении разложения (4.3) и последующем решении двух систем с треугольными матрицами . Лекция 8 Глава 5. Итерационные методы решения систем линейных алгебраических уравнений § 1. Двухслойные итерационные методы В отличие от прямых методов, в которых при отсутствии погрешностей округления точное решение системы находится за конечное число шагов, в итерационных методах всегда определяется некоторая последовательность векторов , которая только в пределе при дает точное решение системы: , где – точное решение рассматриваемой системы . Сходимость последовательности векторов к вектору будем понимать в смысле сходимости по норме: . Из курса функционального анализа известно, что в конечномерном пространстве из сходимости последовательности в одной норме вытекает ее сходимость и в любой другой норме. Поэтому для доказательства сходимости достаточно использовать какую-нибудь одну, наиболее удобную, норму. Из всего многообразия итерационных методов мы рассмотрим только так называемые двухслойные итерационные методы решения системы линейных уравнений , каждый из которых можно записать в следующем каноническом виде , где – некоторая невырожденная матрица, – итерационный параметр -ого шага, – произвольное начальное приближение. Если при всех , то метод называется стационарным. Далее мы будем рассматривать только стационарные двухслойные методы . (5.1) Предполагая матрицу невырожденной, покажем, что если метод (5.1) сходится, то он непременно сходится к точному решению системы . Действительно, пусть существует предел . переходя в равенстве (5.1) к пределу при , получаем . В силу единственности решения системы отсюда следует, что . Таким образом, выбирая тем или иным образом матрицу и параметр , нам следует позаботиться только о самом факте сходимости последовательности приближений . При этом, конечно, желательно обеспечить наиболее быструю сходимость. Это желание, как правило, вступает в противоречие с другим естественным желанием – простоты вычисления каждого очередного приближения по предыдущему приближению . В самом деле, для вычисления приближения на каждом шаге необходимо решить систему . (5.2) Трудоемкость решения этой системы определяется структурой матрицы . Ясно, что наиболее просто решаются системы с треугольными и диагональными матрицами. С другой стороны, можно обеспечить сходимость метода к точному решению за один шаг, выбрав в качестве матрицы матрицу . В этом случае метод фактически оказывается прямым. Но тогда теряет смысл сам итерационный метод – ведь мы отказываемся от применения прямого метода как раз из-за трудоемкости решения исходной системы. На практике в качестве матрицы , как правило, используют треугольную или диагональную матрицу. Вектор называется вектором погрешности -ого шага, а вектор – вектором невязки -ого приближения. В терминах векторов погрешностей метод (5.1) принимает вид , а условие сходимости: . Выразим вектор через вектор : , или , где . Матрица называется матрицей перехода итерационного метода (5.1). Как мы увидим, свойства сходимости метода полностью определяются его матрицей перехода. Важной характеристикой итерационного метода является так называемая асимптотическая скорость сходимости. Пусть – наименьшее число шагов, достаточное для уменьшения нормы начальной погрешности в раз, т.е. . (5.3) Так как , то . Отсюда следует, что неравенство (5.3) будет заведомо выполнено, если потребовать, чтобы имело место неравенство , или . (5.4) Естественно, что говорить о скорости сходимости имеет смысл, только если метод сходится. В дальнейшем будет установлено, что достаточным условием сходимости итерационного метода является неравенство . Имея в виду это условие, из неравенства (5.4) получаем: . Величина и называется асимптотической скоростью сходимости итерационного метода (5.1). Теорема 1 (критерий сходимости итерационного метода). Для того чтобы итерационный метод (5.1) сходился при любом начальном приближении, необходимо и достаточно, чтобы все собственные значения матрицы перехода этого метода были по модулю меньше 1. Теорема 2 (достаточное условие сходимости итерационного метода). Для сходимости итерационного метода (5.1) при любом начальном приближении достаточно, чтобы хотя бы одна из норм матрицы перехода, согласованных с какой-либо векторной нормой, была меньше 1. Доказательство. Пусть – произвольное собственное значение матрицы перехода и – соответствующий ему собственный вектор, т.е. . Если согласована с векторной нормой , то . Отсюда, учитывая, что , получаем: , т.е. модуль любого собственного значения матрицы перехода не превосходит любой ее согласованной нормы. Если , то в силу теоремы 1 итерационный метод сходится. § 2. Примеры двухслойных итерационных методов 1. Метод простой итерации. В этом методе , а – любое положительное число. Следовательно, он записывается в следующем каноническом виде: . В методе простой итерации приближение находится по явным формулам . Поставим следующую оптимизационную задачу: подобрать параметр так, чтобы асимптотическая скорость сходимости была наибольшей (т.е. была наименьшей). Мы приведем решение этой задачи в простейшем случае, когда матрица системы является симметричной и положительно определенной. Пусть . Тогда решение поставленной задачи есть . При этом . 2. Метод Якоби. Представим матрицу в виде , где – нижняя треугольная матрица, диагональные элементы которой равны 0, а поддиагональные элементы совпадают с соответствующими элементами матрицы ; – верхняя треугольная матрица, диагональные элементы которой равны 0, а наддиагональные элементы совпадают с соответствующими элементами матрицы ; – диагональная матрица, на главной диагонали которой стоят соответствующие элементы матрицы . Предположим, что матрица невырождена, т.е. . Полагая и , получаем метод Якоби: , или . Так как – диагональная матрица, то решение находится по явным формулам . Матрица перехода метода Якоби имеет вид . Приведем одно достаточное условие сходимости метода Якоби, которое формулируется в терминах самой матрицы системы, а не матрицы перехода . Для этого нам понадобится определение. Говорят, что матрица имеет строгое диагональное преобладание, если для всех выполняются неравенства . Теорема 3. Если матрица системы имеет строгое диагональное преобладание, то метод Якоби для такой системы сходится при любом начальном приближении. 3. Метод Зейделя. Положим и . В результате получим метод Зейделя , или . Поскольку матрица является треугольной, то решение находится по явным формулам . Матрица перехода метода Зейделя имеет вид . Теорема 4. Если матрица системы имеет строгое диагональное преобладание, то метод Зейделя для такой системы сходится при любом начальном приближении. Теорема 5. Если матрица системы является симметричной и положительно определенной, то метод Зейделя для такой системы сходится при любом начальном приближении. Лекция 9 Глава 6. Численное интегрирование обыкновенных дифференциальных уравнений § 1. Метод Эйлера Рассмотрим задачу Коши для обыкновенного дифференциального уравнения первого порядка (6.1) Общее решение уравнения имеет вид , (6.2) где – произвольная постоянная. Геометрически общее решение (6.2) представляет собой семейство так называемых интегральных кривых, т.е. совокупность линий, соответствующих различным значениям постоянной . Интегральные кривые обладают тем свойством, что в каждой их точке угловой коэффициент касательной равен значению функции в этой точке. Если задать точку , через которую должна проходить интегральная кривая, то с помощью условия из бесконечного семейства интегральных кривых выделяется некоторая определенная интегральная кривая , соответствующая частному решению уравнения . Это частное решение и есть решение задачи Коши (6.1). Задачу Коши не всегда удается решить аналитически, т.е. найти ее точное решение возможно не всегда. В таких случаях используют приближенные (численные) методы. Рассмотрим один из наиболее употребительных численных методов решения задачи Коши (6.1) – метод Эйлера. Разделим отрезок на равных частей точками , где , а . В точке к интегральной кривой уравнения проведем касательную . Учитывая, что , уравнение касательной перепишем в виде . Будем считать, что значение искомого решения в точке приближенно равно значению ординаты точки на касательной, соответствующей абсциссе , т.е. . Обозначая это приближенное значение через , приходим к формуле , называемой формулой метода Эйлера. Заметим, что точка лежит уже не на искомой интегральной кривой, а на некоторой другой, близкой к искомой, интегральной кривой уравнения. Исходя из этой точки и повторяя описанные рассуждения, получим приближенное значение решения задачи Коши в следующей точке : . В общем случае формула метода Эйлера имеет вид . Она позволяет вместо аналитической функции , являющейся точным решением задачи Коши (6.1), найти таблицу ее приближенных значений в точках сетки . Известно, что метод Эйлера относится к методам первого порядка точности, т.е. его погрешность аппроксимации имеет порядок . Рассмотрим несколько более точных методов, являющихся модификациями метода Эйлера. 1. Первый улучшенный метод Эйлера описывается формулами . 2. Второй улучшенный метод Эйлера описывается формулами . Первый и второй улучшенные методы Эйлера имеют погрешность . § 2. Методы Рунге-Кутты Методы Эйлера являются частными случаями методов Рунге-Кутты, к описанию которых мы приступаем. Явный -этапный метод Рунге-Кутты состоит в следующем. Пусть приближенное значение решения задачи (6.1) в точке уже известно. Задаются числовые коэффициенты , , , и последовательно вычисляются значения , , . Очередное приближенное значение решения задачи (6.1) в точке находится по формуле . Коэффициенты , , выбираются из соображений точности получаемой формулы. При и получаем метод Эйлера, рассмотренный в предыдущем параграфе. При получаем семейство методов (6.3) Исследуем погрешность аппроксимации методов (6.3) в зависимости от выбора параметров. Подставляя в последнее равенство выражения для и , получаем: . (6.4) По определению погрешностью аппроксимации называется выражение (6.5) получаемое заменой в (6.4) приближенных значений и точными значениями и соответственно. Найдем порядок погрешности в предположении, что функции и обладают достаточной степенью гладкости. Для этого разложим все величины, входящие в выражение (6.5) по формуле Тейлора в точке . Опуская достаточно громоздкие выкладки, окончательно получаем: Отсюда видно, что если , то методы (6.3) имеют не менее чем первый порядок аппроксимации. Если же дополнительно потребовать, чтобы выполнялись условия и , то эти методы будут иметь уже второй порядок аппроксимации. Таким образом, имеется однопараметрическое семейство двухэтапных методов Рунге-Кутты второго порядка аппроксимации. Это семейство можно записать в виде , (6.6) где , или . В частности, при , получаем первый улучшенный метод Эйлера, рассмотренный в предыдущем параграфе: При , получаем второй улучшенный метод Эйлера: Приведем еще несколько примеров методов Рунге-Кутты более высокого порядка аппроксимации. 1. Метод третьего порядка: , , . 2. Метод третьего порядка: , , . 3. Метод четвертого порядка: , , . 4. Метод четвертого порядка: , , .