Вычислительная математика

Г.В. ВАЩЕНКО ВЫЧИСЛИТЕЛЬНАЯ МАТЕМАТИКА  КУРС ЛЕКЦИЙ Красноярск, 2011 МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ ГОУ ВПО «Сибирский государственный технологический университет» Г. В. ВАЩЕНКО ВЫЧИСЛИТЕЛЬНАЯ МАТЕМАТИКА КУРС ЛЕКЦИЙ для студентов специальности 230105 (230100) “Программное обеспечение вычислительной техники и автоматизированных систем”, 230201 (230400)“Информационные системы и технологии” и направления 230100 “Информатика и вычислительная техника” очной, очной сокращенной, заочной форм обучения Красноярск, 2011 УДК 519.612.2 ББК 22.193 ВАЩЕНКО Г. В. ВЫЧИСЛИТЕЛЬНАЯ МАТЕМАТИКА: Курс лекций для студен-тов специальности 230105(230100), 230201(230400) и направления 230100 очной, очной сокращенной, заочной форм обучения. - Красноярск: СибГТУ, 2011. - 323 с. Излагаются основы численных методов для систем линейных и нелинейных уравнений, приближенного интегрирования и дифференцирования, интерполирования и другие. Описание каждого метода сопровождается представлением вычислительной схемы метода, форм записи, возможного алгоритма реализации, рассмотрением примеров. Рекомендуется студентам, специализирующимся в области информационных технологий. Рецензент: ст. преп. Е. В. Касьянова ( СибГТУ ) Утвержден на заседании кафедры системотехники от 16. 05.06г., протокол № 5 ©ГОУ ВПО “СИБИРСКИЙ ГОСУДАРСТВЕННЫЙ ТЕХНОЛОГИЧЕСКИЙ УНИВЕРСИТЕТ”, 2006 © ВАЩЕНКО Г. В., 2006 СОДЕРЖАНИЕ ВВЕДЕНИЕ.........................................................................................................………. 8 ЛЕКЦИЯ 1 Введение. Основные понятия вычислительной математики......….9 1.1 Предмет вычислительной математики....................................………………….9 1.2 Этапы решения задач на вычислительных машинах………………………….11 1.3 Устойчивость, корректность, сходимость……………………………………..13 1.4 Основные задачи курса………………………………………………………….14 1.5 Заключение……………………………………………………………………….14 1.6 Контрольные вопросы........................................................…………………… 15 ЛЕКЦИЯ 2 Погрешности результатов численного решения задач……………16 2.1 Машинная арифметика……………………..…………………………………...16 2.2 Характеристика точности числа………………………………………………...17 2.3 Классификация и источники погрешностей ………………………………….18 2.4 Заключение……………………………………………………………………….23 2.5 Контрольные вопросы …………………………………………………………..23 ЛЕКЦИЯ 3 Основные сведения алгебры матриц. Классификация численных методов решения систем линейных алгебраических уравнений …………………. 3.1 Системы линейных алгебраических уравнений ………….. 3.2 Нормы векторов и матриц…………………………………………………………. 3.3 Классификация численных методов решения систем линейных алгебраичес-ких уравнений……………………………………………………………………………... 3.4 Обусловленность и число обусловленности…………………………………….. 3.5 Заключение………………………………………………………………………… 3.6 Контрольные вопросы…………………………………………………………….. ЛЕКЦИЯ 4 Конечные методы решения систем линейных алгебраических уравнений. Метод исключения Гаусса. … 4.1 Конечные методы решения систем линейных алгебраических уравнений…… 4.2 Метод исключения Гаусса. ……………………………………………………… 4.4 Заключение……………………………………………………………………… 4.5 Контрольные вопросы………………………………………………………….. ЛЕКЦИЯ 5 Метод исключения Жордана – Гаусса… 5.1 Метод исключения Жордана – Гаусса…………………………………………… 5.2 Заключение……………………………………………………………………… 5.3 Контрольные вопросы………………………………………………………….. ЛЕКЦИЯ 6 Компактная схема исключения………. 6.1 Компактная схема исключения………………………………………………….. 6.2 Метод квадратных корней……………………………………………………….. 6.3 Заключение………………………………………………………………………. 6.4 Контрольные вопросы…………………………………………………………….. ЛЕКЦИЯ 7 Итерационные методы решения систем линейных алгебраических уравнений ………………………………………………………………………………… 7.1 Общий подход к построению итерационных методов………………………… 7.2. Метод простых итераций…………………………………………………………. 7.3 Заключение…………………………………………………………………………. 7.4 Контрольные вопросы……………………………………………………………... ЛЕКЦИЯ 8 Метод итераций Гаусса – Зейделя…………………………………… 8.1 Метод итераций Гаусса – Зейделя………………………………………………… 8.2 Заключение…………………………………………………………………………. 8.3 Контрольные вопросы……………………………………………………………... ЛЕКЦИЯ 9 Численные методы нахождения собственных значений и векто-ров…………………………………………………………………………………………. 9.1 Постановка общей и частичной задачи нахождения собственных значений и векторов……………………………………………………………………………….. 9.2 Классификация методов………………………………………………………….... 9.3 Метод Крылова нахождения собственных векторов……………………………. 9.4 Заключение………………………………………………………………………… 9.5 Контрольные вопросы.........................................................………………………. ЛЕКЦИЯ 10 Основные сведения о приближении функций. Постановка задачи об интерполировании функций алгебраическими многочленами……………….. 10.1 Основные сведения о приближении функций. Постановка задачи………… 10.2 Критерии близости…………………………………………………………….. 10.3 Постановка задачи об интерполировании функций алгебраическими много-членами. Алгебраический многочлен Лагранжа………………………………………. 10.4 Оценка погрешности многочлена Лагранжа. Вычислительная сложность….. 10.5 Заключение……………………………………………………………………….. 10.6 Контрольные вопросы....................... ………………………. ЛЕКЦИЯ 11 Конечные и разделенные разности, свойства. Алгебраические полиномы Ньютона…………………………………………………………………….... 11.1 Конечные и разделенные разности и их свойства……………………. 11.2 Алгебраические полиномы Ньютона……………………………………………. 11.3 Заключение………………………………………………… 11.4 Контрольные вопросы...............................……………………….. ЛЕКЦИЯ 12 Многочлены Чебышева и их свойства. Интерполяция сплай-нами………………………………………………………………………………………. 12.1 Многочлены Чебышева и их свойства…………………………………………. 12.2 Интерполяция сплайнами……………………………………………………….. 12.3 Заключение……………………………………………………………………….. 12.4 Контрольные вопросы.........................................................…………………….. ЛЕКЦИЯ 13: Тригонометрическая интерполяция. Полиномы Фурье ………… 13.1 Основные сведения о тригонометрической интерполяции……………………. 13.2 Интерполяция функций на отрезке……………………………………………… 13.3 Полиномы Фурье…………………………………………………………………. 13.4 Интерполирование таблично заданных функций полиномами Фурье………... 13.5 Заключение………………………………………………………………………. 13.6 Контрольные вопросы.........................................................……………………… ЛЕКЦИЯ 14: Среднеквадратичные приближения. Метод наименьших квадра-тов……………………………………………………………………………………… 14.1 Основные сведения об ортонормированных системах функций…………… 14.2 Ортогональные системы многочленов. Примеры…………………………… 14.3 Метод наименьших квадратов. Примеры……………………………………… 14.4 Заключение………………………………………………………………………. 14.5 Контрольные вопросы.........................................................……………………. ЛЕКЦИЯ 15: Постановка задачи приближенного интегрирования. Квадратур-ные формулы ……………………………………………………………………………. 15.1 Основные сведения теории интегрирования функций………………………… 15.2 Квадратурные формулы прямоугольников, трапеций и Симпсона………….. 15.3 Оценки погрешности квадратурных формул…………………………………. 15.4 Оценка погрешности методом Рунге…………………………………………. 15.5 Заключение…………………………………………………………………….. 15.6 Контрольные вопросы.........................................................……………………. ЛЕКЦИЯ 16: Квадратурные формулы Гаусса и Чебышева……………………… 16.1 Построение квадратурных формул Гаусса. Остаточный член формул Гаусса 16.2 Построение формул Чебышева. Остаточный член формул Чебышева……….. 16.3 Заключение……………………………………………………………………….. 16.4 Контрольные вопросы…………………………………………………………… ЛЕКЦИЯ 17: Постановка задачи численного дифференцирования. Формулы численного дифференцирования……………………………………………………. 17.1 Постановка задачи численного дифференцирования………………………… 17.2 Использование интерполяционных формул …………………………………. 17.3 Вычислительная погрешность формул численного дифференцирования…… 17.4 Заключение………………………………………………………………………. 17.5 Контрольные вопросы.........................................................……………………. ЛЕКЦИЯ 18: Численные методы решения нелинейных уравнений……………….. 18.1 Постановка и основные этапы решения задачи……………………………….. 18.2 Отделение корней. Вычислительная схема. ………………………………….. 18.3 Метод половинного деления. Вычислительная схема…………………………. 18.4 Метод простых итераций. Вычислительная схема. Условия сходимости. 18.5 Заключение……………………………………………………………………….. 18.6 Контрольные вопросы.........................................................……………………… ЛЕКЦИЯ 19: Методы Ньютона, секущих, комбинированный решения нелинейных уравнений…………………………………… 19.1 Метод Ньютона. Вычислительная схема. Выбор начального приближения…. 19.2 Метод секущих. Вычислительная схема. Выбор начального приближения. Условия сходимости……………………………………………………………………… 19.3 Комбинированный метод. Вычислительная схема. Выбор начального приближения. Условия сходимости……………………………………………………… 19.4 Заключение………………………………………………………………………. 19.5 Контрольные вопросы.........................................................…………………….. ЛЕКЦИЯ 20: Численные методы решения систем нелинейных уравнений….. 20.1 Постановка задачи……………………………………………………………….. 20.2 Метод последовательных приближений. Вычислительная схема …………… 20.3 Теорема о существовании и единственности решения. Оценки погрешности 20.4 Заключение……………………………………………………………………….. 20.5 Контрольные вопросы.........................................................……………………… ЛЕКЦИЯ 21: Метод Ньютона решения систем нелинейных уравнений ……….. 21.1 Метод Ньютона. Вычислительная схема…………………………………. 21.2 Выбор начального приближения. Оценки и условия сходимости……………. 21.3 Модификации вычислительных схем метода…………………………………. 21.4 Заключение………………………………………………………………………. 21.5 Контрольные вопросы........................................................………………………. ЛЕКЦИЯ 22: Приближенное интегрирование обыкновенных дифференциаль-ных уравнений…………………………………………………………………………. 22.1 Основные сведения теории обыкновенных дифференциальных уравнений 22.2 Постановка задачи численного интегрирования задачи Коши……………… 22.3 Классификация приближенных методов………………………………………. 22.4 Метод последовательных приближений…………………………………….. 22.5 Заключение……………………………………………………………………. 22.6 Контрольные вопросы.........................................................…………………… ЛЕКЦИЯ 23: Численные методы решения задачи Коши: метод Эйлера, методы Рунге-Кутты ……………………………………………………………………………. 23.1 Метод Эйлера и его модификации для решения задачи Коши. Вычислительная схема. Геометрическая интерпретация. Примеры…………………... 23.2 Методы Рунге-Кутты. Вычислительная схема. Явные и неявные схемы. Оценки погрешности. Контроль точности. Примеры………………………………….. 23.3 Заключение………………………………………………………………………. 23.4 Контрольные вопросы.........................................................……………………… ЛЕКЦИЯ 24: Численные методы решения систем обыкновенных дифферен-циальных уравнений …………………………………………………………………… 24.1 Постановка задачи Коши для систем обыкновенных дифференциальных уравнений………………………………………… ………………………………………. 24.2 Метод последовательных приближений. Методы Эйлера, Рунге –Кутты… 24.3 Многошаговые методы Адамса. Вычислительная схема методов …………. 24.4 Заключение……………………………………………………………………… 24.5 Контрольные вопросы.........................................................……………………. ЛЕКЦИЯ 25: Численные методы решения уравнений в частных производных.. 25.1 Основные сведения об уравнениях в частных производных…………………. 25.2 Классификация дифференциальных уравнений в частных производных. Примеры……………………………………………………………... 25.3 Начальные и краевые (граничные)условия, как условия однозначной характеристики описываемого физического процесса. Примеры………………….. 25.4 Конечно -разностные схемы. Простейшие разностные схемы. Примеры……. 25.5 Заключение…………………………………………………………… 25.6 Контрольные вопросы........................................................... ЛЕКЦИЯ 26: Конечно-разностные методы численного решения……. 26.1 Аппроксимация -сетка и шаблон…………………………………………………………... 26. 2 Методы составления схем. Устойчивость. Примеры………………………………… 26.3 Заключение………………………………………………………………………. 26.4 Контрольные вопросы.........................................................……………………. ЛЕКЦИЯ 27: Ортогональные системы функций…………………………………. 25.1 Ортогональные системы функций……………………………………………… 25.2 Обобщенный полином Фурье……………………………………………………. 25.3 Прямое и обратное преобразование Фурье. Свойства преобразований………. 25.4 Вычислительные схемы численного преобразования. Примеры……………. 25.5 Заключение…………………………………………………………. 25.6 Контрольные вопросы……………………………………………………… ЛЕКЦИЯ 28: Дискретное преобразование Фурье……………………………….. 26.1 Быстрое преобразование Фурье. Примеры…………………………….. 26.2 Вычисление преобразования Фурье. Вычислительные схемы………………. 26.3 Алгоритмы преобразований Фурье…………………………………………….. 26.4 Заключение………………………………………………………………. 26.5 Контрольные вопросы.........................................................…………………… ЛЕКЦИЯ 29: Параллельные алгоритмы матричных вычислений……….. 29.1 Основные сведения о параллельных моделях программирования…………. 29.2 Модель передачи сообщений. Этапы разработки параллельного алгоритма 29.3 Количественные характеристики быстродействия. Закон Амдала……… 29.4 Примеры топологий с использованием модели передачи сообщений (MPI). Кластеры…………………………………………………………………. ……………… 29.5 Основные сведения о программировании в технологии MPI……………….. 29.6 Заключение……………………………………………………………………… 29.7 Контрольные вопросы........................................................... ЛЕКЦИЯ 30: Параллельные алгоритмы матрично – векторных произведений . 30.1 Параллельные алгоритмы матрично – векторных произведений. Схемы размещения и хранения матриц и векторов: блочная по строкам, блочная по столбцам, циклические схемы по строкам и столбцам…………………….. 30.2 Параллельные алгоритмы матрично – векторных произведений. Произве-дение в топологии звезда, кольцо и двумерная решетка……………………………….. 30.3 Заключение………………………………………………………………………... 30.4 Контрольные вопросы.........................................................……………………. ЛЕКЦИЯ 31: Параллельные алгоритмы матричных произведений……….. 31.1 Параллельные алгоритмы матричных произведений. Произведение в топологии звезда, кольцо и двумерная решетка……………………………………. 31.2 Параллельный алгоритм метода исключения Гаусса и итерационных методов решения систем линейных алгебраических уравнений. Представление алгоритмов в топологии звезда, кольцо……………………………….……. 31.3 Заключение……………………………………………………………. 31.4 Контрольные вопросы .........................................................…………….. ЛЕКЦИЯ 32: Общая схема метода статистических испытаний ………………… 32.1 Метод статистических испытаний или метод Монте-Карло. Обобщенная схема. Построение статистических оценок для искомых величин…………….. 32.2 Методы моделирование случайных чисел с различными законами распре-деления……………………………………………………………………………………. 32.3 Случайные процессы и их моделирование. Примеры………………………. Заключение……………………………………………………………………………. Контрольные вопросы.........................................................…………………………. ЛЕКЦИЯ 33: Моделирование на основе метода статистических испытаний…. 33.1 Имитационное моделирование процессов в информационных системах. Примеры……………………………………… ………………………………………….. 33.2 Формирование входных и выходных данных моделирования………………. 33.3 Схема функционирования алгоритмов и их моделирование. Примеры…….. 33.4 Заключение………………………………………………………………………. 33.5 Контрольные вопросы.........................................................…………………… ЛЕКЦИЯ 34: Обзор современных задач вычислительной математики, математического моделирования и вычислительных технологий……………….. 34.1 О задачах математической физики……………………………………………… 34.2 О задачах моделирования информационных систем……………………….. 34.3 О вычислительных технологиях……………………………………………….. 34.4 Заключение………………………………………………………………………... 34.5 Контрольные вопросы.........................................................……………………… ЗАКЛЮЧЕНИЕ....................................................................................…………………. ПЕРЕЧЕНЬ ОБОЗНАЧЕНИЙ.......................................……………………………….. БИБЛИОГРАФИЧЕСКИЙ СПИСОК...........................................................…………. ВВЕДЕНИЕ Настоящий курс лекций является важным элементом в структуре учебного плана дисциплины “Вычислительная математика” и раскрывает содержание основных компонентов подготовки по вычислительной математики специалистов в области информационных технологий. Цель курса - организация самостоятельной работы студентов по овладению методов вычислительной математики, подходов в алгорит-мизации и разработке программного математического обеспечения для решения инженерных задач. Задачи курса лекций состоят в формировании представлений о возможностях вычислительных технологий и математического моделирования, а также способствовать приобретению навыков проведения численных экспериментов, моделирования и разработки программного математического обеспечения в системах визуального программирования. Приводимые в курсе основные сведения о численных методах решения различных прикладных задач могут быть востребованы студентами как при освоении дисциплин «Методы оптимизации», «Системы цифровой обработки сигналов» «Компьютерная графика» и других, так и при выполнении курсовых работ по данным дисциплинам и написании дипломной работы. Поскольку данный курс лекций ориентирован на студентов, специализирующихся в области информационных технологий, программного обеспечения и автоматизированных систем, то основной упор делается на алгоритмическую основу реализации вычислительной схемы метода и условиях применения. Изложение проводится на доступном для студентов этих специальностей уровне. При необходимости напоминаются основные сведения из курса высшей математики. Для лучшего понимания и активного усвоения каждая лекция завершается перечнем контрольных вопросов. Лекция 1 Введение. Основные понятия вычислительной математики План: 1.1 Предмет вычислительной математики 1.2 Этапы решения задач на вычислительных машинах 1.3 Устойчивость, корректность, сходимость 1.4 Основные задачи курса 1.5 Заключение 1.6 Контрольные вопросы 1.1 Предмет вычислительной математики Вычислительная математика - раздел математики, включающий в себя вопросы, связанные с применением вычислительных машин для решения задач науки, техники и практически всех областей человеческой деятельности. Современная вычислительная математика включает в себя три основные части: 1 теория вычислительных методов; 2 средства автоматизации вычислений; 3 вспомогательные средства для организации вычислений. Рисунок 1- Основные составные части вычислительной математики Охарактеризуем содержательно каждую этих частей вычисли-тельной математики. Теория вычислительных методов. Теория вычислительных методов является наукой, которая имеет применение там, где приходиться встречаться с числами или рассматривать явления и процессы, подчиняющиеся количественным законам. Такие законы изучаются во многих науках: математике, физике, механике, астрономии, экономике, биологии, медицине, теории управления и регулирования и других областях человеческой деятельности. В каждой из этих областей применения могут возникать свои специфические задачи, требующие разработки особых, приспособленных к ним вычислительных методов. Средства автоматизации вычислений. Средства автоматизации вычислений дают возможность хранить информацию, выполнять арифметические операции, выводить результата и др. Среди этих средств, на сегодняшний день, основным является вычислительная машина. Средства для управления вычислительной машиной. Средства для управления вычислительной машиной классифицируются исходя из состава вычислительной машины, разделяющегося на аппаратное и программное обеспечение. К программному обеспечению, требующемуся для решения задач вычислительной математики относятся алгоритмические языки программирования и системы программирования, системы предназначенные для расчета и решения математических задач (такие системы, как MathCad, Matlab, Maple и др.). Для управления самим процессом вычислений и работой вычислительной машины требуется операционная система, подсистемы диспетчеризации и т.п. Содержание предмета вычислительной математики, в общем случае, состоит в том, чтобы задачи, связанные с недоступными вычислительной машине понятиями, перевести на язык, понятный машине. Рассмотрим пример, иллюстрирующий содержание предмета вычислительной математики. Пусть требуется вычислить значение определенного интеграла , где f(x) – некоторая непрерывная на отрезке [a, b] функция. Этот интеграл непонятен машине, если функция не определяется однозначно конечным набором чисел. Пусть функция является многочленом степени n, в этом случае величина интеграла полностью определяется набором коэффициентов многочлена и, по сути, является отображением набора из (n +1) чисел на множество действительных чисел. При этом, коэффициенты многочлена и значение интеграла должны быть понятными вычислительной машине, другими словами, должны относиться к области понимания машины. Таким образом, вычислительная машина может иметь дело только с таким соответствием, которое отображает подмножество заданного конечного множества чисел в другое его подмножество. Этот пример показывает и содержание метода вычислительной математики, именно, всякая задача, решаемая на вычислительной машине, сводится к нахождению такого соответствия. Более строго содержание метода вычислительной математики можно определить следующим образом. Метод вычислительной математики. Если считать, что соответствие между условием задачи U и ее решением Y есть оператор A, причем Y = A(U), то речь идет о замене области определения, области значений и самого оператора дискретной областью определения Ũ, дискретной областью значений Ỹ и оператором Ã, осуществляющим соответствие между ними. 1.2 Этапы решения задач на вычислительных машинах Последовательность шагов при использовании вычислительных машин для решения задач. Постановка задачи. Данный этап состоит в содержательной (физической) постановке задачи, определение конечных целей. Построение математической модели или математическая формулировка задачи. На этом этапе формулируется математическая модель, адекватно описывающая законы явлений и процессов. Построение или выбор модели основан на знании соответствующих разделов математики. Выбор или разработка численного метода. На этом этапе выбирается или разрабатывается численный метод, позволяющий свести задачу к вычислительному алгоритму. Разработка алгоритма и построение схемы программы. На данном этапе процесс решения задачи (вычислительный процесс) записывается в виде последовательности элементарных арифметических и логических операций, т.е. формируется алгоритм решения задачи. На этом этапе должна быть разработана обобщенная конструктивная схема программного обеспечения, связывающая функциональные блоки разрабатываемого программного обеспечения. Программирование. Составление программы, реализующей алгоритм решения задачи, производится в системе программирования, в которой проводится компиляция для получения исполняемого кода программы. Отладка. Отладка составленной программы включает контроль программы по входным данным, диагностику, решение контрольных задач с целью получения достоверных результатов работы программы. Проведение расчетов. На этом этапе формируются исходные данные для расчетов и проводятся расчеты. Анализ результатов. Формирование отчетов и анализ, полученных результатов. Это общий подход, объем каждого этапа которого зависит от размерности решаемой задачи, целей и конечного результата. Сделаем несколько замечаний о численных методах. Численные методы. Прежде всего напомним, что для решения математических задач используются такие основные группы методов, как графические, аналитические и численные. Суть графических методов состоит в том, что решение находится с помощью геометрических построений. При использовании аналитических методов решение задачи удается выразить с помощью формул, но на практике такие задачи встречаются редко. Для решения сложных математических задач, являющихся моделью физических, технических процессов и явлений, основным способом решения на сегодняшний день являются численные методы, которые дают возможность свести решение задачи к выполнению конечного числа арифметических операций над числами, а результаты получаются в виде числовых значений или последовательности числовых значений. 1.3 Устойчивость, корректность, сходимость Рассмотрим на содержательном уровне такие понятия вычислительной математики, как устойчивость, корректность, сходимость. Строгие определения всех понятий будут формулироваться при рассмотрении соответствующих методов решения конкретных задач. Устойчивость. Задача называется устойчивой по исходному параметру x , если решение y непрерывно зависит от него, т.е. малое приращение исходной величины dx приводит к малому приращению dy. Другими словами, малые погрешности в исходной величине приводят к малым погрешностям в результате расчетов. Отсутствие устойчивости означает, что даже незначительные погрешности исходных данных приводят к большим погрешностям в решении или даже к неверному результату. Подобные неустойчивые задачи называют чувствительными к погрешностям исходных данных. Корректность. Задача называется поставленной корректно, если для любых значений исходных данных из некоторого набора ее решение существует, единственно и устойчиво по исходным данным. Все неустойчивые задачи являются некорректно поставленными. В последующих лекциях рассмотрим примеры и более строго определим это важное понятие. Сходимость. Понятие сходимости используется при анализе точности вычислительного процесса. Сходимость является критерием близости получаемого численного решения задачи к искомому. Рассмотрим понятие сходимости на примере метода последовательных приближений состоящего в том, что для решения задачи и нахождения искомого решения строится последовательность значений. Если эта последовательность имеет предел при неограниченном возрастании числа значений, то численный метод будет сходящимся. В общем случае, для получения решения задачи с требуемой точностью постановка задачи должна быть корректной, а численный метод, выбранный для решения должен быть устойчивым и сходящимся. 1.4 Основные задачи курса Задачами курса является формирование представлений о возможностях вычислительных технологий и математического моделирования, а также способствовать приобретению навыков проведения численных экспериментов, моделирования и разработки программного математического обеспечения в системах визуального программирования. В курсе будут рассмотрены следующие задачи: 1 Задача вычисления функций. 2 Численные методы решения задач линей ной алгебры (задача решения систем линейных алгебраических уравнений, обращение матриц, нахождение собственных векторов и собственных значений). 3 Интерполирование функций. 4 Численное интегрирование и дифференцирование. 5 Задача решения алгебраических и трансцендентных уравнений. 6 Задача решения систем нелинейных уравнений. 7 Решение обыкновенных дифференциальных уравнений и уравнений в частных производных. 8 Метод статистических испытаний и его применение в модели-ровании. 1.5 Заключение В данной вводной лекции были рассмотрены на содержательном уровне предмет вычислительной математики, определены основные понятия такие как устойчивость, корректность, сходимость. Рассмотрен основной метод вычислительной математики, определены задачи самого курса. В лекции не были приведены примеры математических моделей и математического моделирования, вычислительный эксперимент. Это будет сделано при рассмотрении соответствующих разделов вычислительной математики, некоторые примеры моделей процессов и явлений будут рассмотрены на практических и лабораторных занятиях. На следующей лекции поговорим о погрешности результатов решения задач на вычислительных машинах, об оценках этой погреш-ности. 1.6 Контрольные вопросы 1. Перечислить составные части вычислительной математики. 2. Сформулировать определение устойчивости. 3. Описать метод вычислительной математики. 4. Перечислить основные этапы решения задач на вычислительной машине. 5. Сформулировать определение сходимости и корректности. 6. Перечислить основные группы методов для решения математических задач. 7. Сформулировать определение вычислительного алгоритма. Лекция 2 Погрешности результатов численного решения задач План: 2.1 Характеристика точности числа 2.2 Машинная арифметика 2.3 Классификация и источники погрешностей 2.4 Заключение 2.5 Контрольные вопросы 2.1 Характеристика точности числа Для оценки и учета погрешностей действий используются разные подходы. В основе этих подходов лежат такие характеристики точности числа, как абсолютная и относительная погрешности и связанные с ними оценки. Рассмотрим эти два понятия. Характеристики точности. Для характеристики точности величи-ны используются понятия абсолютной и относительной погрешности. Пусть - приближенное значение числа, а x - точное его значение. Абсолютной погрешностью (ошибкой) приближенного числа называет- ся разность между точным его значением и приближенным, . Точное значение x, как правило, неизвестно, поэтому указываются границы или оценки, в которых изменяется абсолютная погрешность. Оценкой или предельной абсолютной погрешностью приближенного числа является такое число , что выполняется неравенство . Относительной погрешностью (ошибкой)  приближенного числа называется отношение абсолютной погрешности этого числа к модулю соответствующего приближенного значения, при  0. Оценкой или предельной относительной погрешностью приближенного числа является такое число , что выполняется неравенство . Для величин, близких по значению к единице абсолютная и относительная ошибки почти одинаковы. Для очень больших или очень малых величин относительная и абсолютная ошибки представляются разными числами. Пусть, например, точное значение некоторой величины x равно 0.00009, x = 0.00009, приближенное значение = 0.00008, тогда абсолют-ная погрешность  составит  = 0.00009 - 0.00008 = 0.00001, а относительная  = 0.00001/0.00008 = 0.125. Пусть точное значение некоторой величины x равно 100200, x = 100200, приближенное значение = 100000, тогда  = 100200 - 100000 = 200, а  = 200/100000 = 0.002. Замечание. Правила оценки абсолютной и относительной погрешностей для арифметических операций рассмотрим на практических занятиях. 2.2 Машинная арифметика Машинная арифметика. Известно, что любое неотрицательное число x может быть представлено в виде степенного ряда x = an10n + an-110n-1 + ... a0 +a-110-1 + a-210-2 + ..., где коэффициенты ai могут принимать значения 0,1, 2,..., 9. Если перечис- лить подряд все коэффициенты, указать положение запятой и приписать числу некоторый знак, то число x запишется в следующей системе записи x =  an an-1... a0 a-1 a-2.... . Запись в виде степенного ряда справедлива и для произвольной системы счисления. Если зафиксировать некоторое целое положительное число p и целые числа 0, 1, ...,p-1 , то для произвольного неотрицатель- ного числа справедливо представление x=bnpn +bn-1 pn-1 +...+b0+b-1 p-1 + b-2 p-2 +..., где каждый из коэффициентов bi может принимать одно из значений 0, 1, ...,p-1. Справедлива и запись x = bn bn-1... b0, b-1 b-2.... . Такие способы изображения чисел называется позиционными системами счисления. Каждое число здесь занимает определенную позицию, а отсчет определяется положением коэффициента b0. Значащими цифрами числа называются все цифры в его записи, начиная с первой ненулевой слева. Например, в числе 0.003070 первые три нуля не являются значащими, они только устанавливают десятичные разряды других цифр. Значащую цифру называют верной, если абсолютная погрешность не превосходит половины единицы разряда, соответствующего этой цифре. Например, для точного числа 46.97, число 47.00 является приближенным с тремя верными знаками, так как  =  46. 97- 47.00  = 0.03< 0.510-1. О целых машинных числах. Целые числа в вычислительных машинах представляются целыми числами, принадлежащими некоторому отрезку [N-, N+]. Например, для 16- битовой архитектуры вычислительной машины под целые числа без знака отводится отрезок [0, 65535] или [0, 216-1], а для целых со знаком [-32767, 32768] или [-215-1, 215]. Арифметические операции с целыми машинными числами выполняются точно, если результат не выходит за этот отрезок [N-, N+]. В противном случае фиксируется ошибка, а машинный результат операции не определен. О вещественных машинных числах. В вычислительных машинах используется, в основном, представление вещественных чисел с плаваю- щей запятой в экспоненциальной форме. Ненулевые вещественные числа имеют вид x = mb p, где b - основание системы счисления, p - порядок (или экспонента, или показатель числа), m - мантисса (дробная часть), b-1   m  b0 = 1. Если под мантиссу m отводится n элементов в b -ичной системе, а под порядок k элементов, то вещественное число x представляется конечным числом вида x  fl(x) = (1b-1 + 2 b-2 + ... + n b-n)be, где -bk  e  bk -1, 1 1  b-1, 0  i  b-1 для всех i от 2 до n, fl -сокращение от float - плавающий. На рисунке 1 схематично показана структура представления вещест- венного числа с плавающей запятой. Рисунок 1- Схема представления вещественного числа с плавающей запятой Значения bd, d = bk определяют диапазон изменения вещественных чисел с плавающей запятой. При этом числа из отрезка [-b-d , b-d] машина заменяет нулем, поэтому он называется машинным нулем, а числа, выходящие за границу отрезка [-bd-1 , bd-1], машина не воспринимает. Этот отрезок называется машинной бесконечностью. Еще одной характеристикой вещественного числа с плавающей запятой является характеристика, называемая машинным эпсилоном  или macheps, которая определяет расстояние между единицей и ближайшим следующим за ней числом с плавающей запятой. За машинный эпсилон  принимается величина  = 1 b1-n, здесь, как и выше, b - основание системы счисления, n-число разрядов, отводимых под мантиссу числа. Число  является оценкой относительной и абсолютной погрешности чисел с плавающей запятой. Пусть x -некоторое число, а fl(x) - запись этого числа с плавающей запятой, тогда для относительной погрешности справедлива оценка , а для абсолютной погрешности оценка . Как видно из этих неравенств, относительная точность одинакова и не зависит от количества разрядов, отводимых под мантиссу числа с плавающей запятой. Абсолютная погрешность, напротив, неодинакова в разных частях числового диапазона. Таким образом, множество машинных вещественных чисел характеризуется следующими четырьмя целыми параметрами: основанием системы счисления b  2, точностью n, нижней и верхней границами порядка. Машинные вещественные числа образуют конечное подмножество вещественных чисел, причем все машинные числа рациональны. Арифметические операции с вещественными машинными числами. Пусть x и z -два машинных вещественных числа, пусть знак  обозначает одну из четырех арифметических операций +, , , / и пусть результат выполнения операции с числами x и z на вычислительной машине обозначается через (x  z)м, тогда возможны следующие ситуации при выполнении операции на вычислительной машине: если x  z < b-d-1 , d = bk , x  zм= 0 есть машинный нуль; если x  z > bd-1(1 - b-n), то возникает переполнение (машинная бесконечность); если b-d-1  x  z  bd-1(1 - b-n) и x  z - машинное число, то x  zM=x  z; если b-d-1   x  z  bd-1(1 - b-n), но x  z - не машинное число, то x  zM будет равно одному из ближайших к x  z машинных чисел. Ошибкой округления при выполнении арифметических операций  называется величина  x  zM - x  z   b-d-1, если x  z < b-d-1 , d = bk и x  z, если b-d-1  x  z  bd-1(1 - b-n). Из определения ошибки округления следует оценка  x  zM - x  z   max{ b-d-1 , x  z} b-d-1 + x  z . Для вещественных машинных чисел арифметические операции не обладают свойствами дистрибутивности и ассоциативности, именно, ((1 + ( /b)M)M -1)M = 0, ((1-1)M + ( /b)M)M =  /b, здесь  -машинный эпсилон. В силу того, что арифметические операции не обладают свойствами ассоциативности и дистрибутивности, анализ ошибок округлений рассматривать не будем. 2.3 Классификация и источники погрешностей При численном решении математических задач, как правило, возникают погрешности, которые подразделяются на три основных типа: погрешность задачи, погрешность метода и погрешность округ ления. Охарактеризуем каждую из них на отдельных этапах решения задачи. Погрешность задачи. Эту погрешность порождают исходные данные задачи. Эти погрешности называются неустранимыми погреш ностями, так как они не могут быть уменьшены при вычислениях. Погрешность метода. Погрешность метода связана с заменой, например, при вычислении функции усеченным рядом, интерполиро-ванием табличных данных и т.п. Эта погрешность может быть умень шена, т.е. она регулируемая. Более подробно погрешности методов будем рассматривать при анализе конкретных численных методов. Погрешность округления. Погрешность округления связана с конечностью разрядной сетки вычислительной машины, т.е. для записи каждого числа количество разрядов ограничено. По сути, это означает, что подмножество действительных чисел, которым оперирует вычислительная машина, является конечным. Как было сказано в 2.2, наименьшее по модулю машинное число M0, отличное от нуля, называется машинным нулем. Далее, при появлении в ходе вычислений числа меньшего по модулю M0, ему присваивается нулевое значение. Наибольшее по модулю машинное число M называется машинной бесконечностью. При появлении числа, большего M, происходит переполнение. Таким образом, абсолютные значения всех ненулевых действительных чисел, представимых конкретной вычислительной машине, лежит между M0 и M. Если число a не представимо в машине точно, оно округляется, т.е. заменяется на близкое ему число a~, представимое точно. Одним из способов округления является отбрасы-вание разрядов мантиссы, выходящих за пределы разрядной сетки. 2.4 Заключение На сегодняшней лекции были рассмотрены характеристик точности результатов вычислений, классифицировали основные источ ники погрешностей при решении задач на вычислительных машинах и рассмотрели машинную арифметику, арифметические действия с числами, представимыми в вычислительной машине. Правила оценки абсолютной и относительной погрешностей рассмотрим на практичес ких занятиях. 2.5 Контрольные вопросы 1. Дать определение абсолютной погрешности. 2. Перечислить и охарактеризовать виды погрешностей. 3. Сформулировать понятие машинного нуля и машинной бесконечности. 4. Дать определение понятию машинный эпсилон . 5. Дать определение относительной погрешности. 6. Показать структура представления вещественного числа с плавающей запятой. 7. Дать определение значащей и верной цифре. Лекция 3 Основные сведения алгебры матриц. Классификация численных методов решения систем линейных алгебраических уравнений План: 3.1 Системы линейных алгебраических уравнений 3.2 Нормы векторов и матриц 3.3 Обусловленность и число обусловленности. 3.4 Классификация численных методов решения систем линейных алгебраических уравнений 3.5 Заключение 3.6 Контрольные вопросы 3.1 Системы линейных алгебраических уравнений Система вида , (3.1) где aij ( 1  i  m, 1  j  n) и b1, ..., bm- заданные числа, x1, ..., xn - неизвестные, называется системой линейных алгебраических уравнений. Линейная система (3.1) допускает следующие формы записи: , или Ax = b, где A = (aij) - m  n-матрица системы, x = (x1, ..., xn )T - вектор-столбец, b = (b1, ..., bm)T - вектор-столбец свободных коэффициентов; , или где ap, p = 1, 2, ..., m - вектор- строки, x = (x1, ..., xn)T - вектор-столбец неизвестных. Система (3.1) называется однородной, если все bi, i =1, 2, ..., m равны нулю. В противном случае, т.е. при b 0, система называется неоднородной. Система (3.1) называется совместной, если существует вектор x0 = (x10, x20, ..., xn0) такой, что Ax0 = b, тогда вектор x0 называется решением системы, в противном случае - несовместной. Множество решений системы (3.1) характеризуется следующим образом. Это множество может быть пустым, в этом случае система несовместна. Множество состоит из единственного решения, в этом случае система называется определенной. Множество имеет более одного решения, в этом случае система называется неопределенной. Если присоединить столбец b к матрице A, то полученная матрица B называется расширенной матрицей системы (3.1). B =. Условия совместности для линейной системы (3.1) даются следующей теоремой. Теорема Кронекера -Капелли. Система (3.1) совместна в том и только в том случае, когда ранги расширенной и исходной матриц совпадают, т.е. rangA = rangB. Системы с квадратной матрицей. Пусть матрица системы (3.1) квадратная, n = m и невырожденнная, detA  0 (это условие обеспечивает совместность и определенность), тогда единственное решение x1, x2, ..., xn может быть получено по формулам Крамера , где Aji - алгебраическое дополнение элемента aij. При больших n эти формулы не используются, так как число арифметических операций возрастает с ростом размерности и состав-ляет n2n!. Однородные системы линейных уравнений. Пусть у системы (3.1) все bi, i =1, 2, ..., m равны 0, тогда система однородная, т.е. , или Ax = 0, (3.2) где А = (аij), x - искомый вектор. Если detA  0, то в силу формул Крамера система (3.2) имеет единственное решение x = 0. Для того чтобы система (3.2) имела и ненулевые решения необходимо, чтобы detA = 0. В этом случае совокупность всех решений системы (3.2) образует векторное пространство, которое называется пространством решений. Размерность этого пространства определяется следующей теоремой. Теорема. Если n - число неизвестных однородной системы (3.2) и r - ранг матрицы А, то размерность пространства решений равна k = n- r. Любая линейно независимая система решений в числе k = n - r называется фундаментальной для системы (3.2). Для нахождения решения однородной системы (3.2) достаточно найти фундаментальную систему ее решений. Тогда все решения системы (3.2) будут линейными комбинациями фундаментальной системы. В общем случае неоднородной системы (3.1), общее ее решение есть сумма произвольного частного решения и общего решения соответствую щей ей однородной системы. Справедлива теорема. Теорема. Пусть система (3.1) совместна, то есть rangA = rangB, тогда общее ее решение представляется в виде суммы произвольного частного решения и общего решения соответствующей ей однородной системы. Любая линейная система Ax = b может быть симметризована умножением слева на матрицу АТ. 3.2 Нормы векторов и матриц Норма вектора. В пространстве Rn скалярной характеристикой вектора является норма, которая дает возможность сравнивать векторы по величине и, которая, в общем случае, обобщает понятие длины вектора. Понятие нормы определяется аксиоматически. Нормой вектора x называется действительное число (функционал) , удовлетворяющее следующим условиям: 1)  0, = 0  x = 0 (неотрицательность ); 2) а -произвольное число (скаляр); 3) для произвольных x и y. Норма вектора определяется различными способами в зависимости от целей исследования и условий задачи, но при любом определении норма должна удовлетворять вышеуказанным условиям 1), 2), 3). Наиболее распространенными в практике вычислений являются следующие три нормы. 1. Равномерная норма 2. Абсолютная норма 3. Евклидова (сферическая) норма Две векторные нормы и называются эквивалентными, если существуют такие положительные вещественные числа с1 и с2 , что для любого вектора x выполняются неравенства с1  с2, причем с1 и с2 не зависят от выбора x Неравенство Коши -Буняковского. Для любых векторов x, y справедли- во неравенство (x, y)  . Норма матриц. Рассмотрим теперь множество всех n x n -матриц. Для матриц понятие нормы также вводится аксиоматически. Нормой квадратной матрицы А называется действительное число (функционал) , удовлетворяющее следующим условиям: 1)  0, = 0  А = 0 (нулевая матрица ); 2) а - произвольное число (скаляр); 3) для произвольных A и B; 4) для произвольных A и B. Как и норму вектора, норму матрицы можно вводить различными способами, но поскольку матрицы и векторы, как правило, рассматривают- ся совместно, то матричная норма должна соотноситься с векторной, т.е. матричная норма должна удовлетворять условию согласованности и условию подчиненности. Условие согласованности. Норма матрицы А называется согласованной с нормой вектора , если для любых x и А справедливо неравенство . Условие подчиненности. Норма матрицы А называется подчиненной норме вектора x, если выполняется равенство . Условие подчиненности означает, что за подчиненную норму матрицы А при заданной норме вектора x принимается максимальная норма векторов Ax, когда x пробегает множество всех векторов, норма которых равна единице. Матричные нормы соответствующие векторным нормам , имеют вид , . Эти нормы легко вычисляются. Евклидовая норма матрицы вычисляется сложнее и определяется как Геометрическая интерпретация. Матричную норму можно интер-претировать как максимальную длину вектора, полученного в результате действия преобразования A к векторам единичной сферы. Энергетические нормы. Важный класс норм образуют так называ-емые энергетические или А -нормы, определяемые выражением где А - некоторая заданная симметрическая положительно определен ная матрица. Евклидовая норма, в этом случае, получается как частный случай при А= Е. 3.3 Обусловленность и число обусловленности Рассмотрим понятие обусловленности. В вычислительной алгебре важным понятием является число обусловленности матрицы, которое может рассматриваться как мера чувствительности решения к возму щениям системы. Под возмущениями системы понимаются погреш ности округления в результате арифметических действий и погрешнос ти в исходных данных. Рассмотрим линейную систему Ax = b (3.3) с квадратной невырожденной матрицей А. В результате решения системы получается, в большинстве случаев, приближенное решение = x +x. Для оценки того насколько приближенное решение отличается от точного решения x применяется обратный анализ погрешностей. Пусть в системе (3.3)вместо истинного вектора b используется приближенный вектор b+b, т.е. “возмущается” только вектор b. В результате решения системы (3.3)получается вектор = x +x, который можно рассматривать как точное решение возмущенной системы А= или А= b+ b. В этом случае справедлива оценка [4], связывающая относительные погрешности вектора-решения и вектора b . (3.4) Пусть теперь в исходной системе (3.3)получила возмущение матрица А, тогда вектор = x +x рассматривается как точное решение возмущенной системы = b или (А+А)( x +x) = b. В этом случае справедлива следующая оценка . (3.4) Если в системе (3.3)возмущены и матрица А и вектор b, т.е. по сути решается система = или (А+А)( x +x) = b+ b, то справедлива оценка , (3.5) при условии, что < 1 . В неравенствах (3.3), (3.4) и (3.5) положительное число называется числом обусловленности и обозначается как cond(A), т.е. cond(A) = (cond от англ. conditioned -обусловленный). Число обусловленности является коэффициентом связи между относи-тельными погрешностями вектора -решения, вектора b и матрицей А. Из неравенств (3.3), (3.4) и (3.5) следует, что чем больше число обусловленности, тем сильнее влияют на решение исходной системы погрешности в исходных данных. Таким образом, число обусловлен-ности можно рассматривать как меру чувствительности решения к погрешностям системы. В качестве примера рассмотим неравенство (3.4). Из этого неравен-ства видно, что даже если вектор -невязки r = b - A мал, то относительные возмущения в решении могут быть большими, если число cond(A) значительно больше единицы. Системы уравнений, у которых число обусловленности cond(A) велико, называются плохо обусловленными. Для числа обусловленности справедлива оценка снизу, именно, cond(A) =  , другими словами, число обусловленности не может быть меньше единицы. В общем случае, вопрос о том является ли система плохо обусловлен- ной или нет, решается для каждой конкретной системы. Геометрически понятия обусловленности можно рассмотреть на примере системы двух уравнений с двумя неизвестными. Плохая обусловленность системы в этом случае означает, что прямые, являющиеся геометрическим образом уравнений, пересекаются на координатной плоскости под очень острым углом. Таким образом, если возмущения правой части уравнения рассматривать как параллельный перенос, а возмущения матрицы как поворот прямых, то небольшие погрешности (возмущения) приводят к значительному перемещению точки пересечения прямых, которая является геометрическим образом решения системы уравнений. В многомерном случае геометрическая картина может быть сложней. Например, для трех переменных возможен случай плохой обусловленности, когда соответствующие трем уравнениям плоскости пересекаются под большими углами (т.е. далеки от параллельности), но линии их попарного пересечения почти параллельны. Используя ту или иную конкретную норму можно определять различные числа обусловленности. Для вырожденных матриц cond(A) = . 3.4 Классификация численных методов решения систем линейных алгебраических уравнений Решение систем линейных алгебраических уравнений является одной из основных задач как в линейной алгебре, так и во многих приложениях. Численные методы линейных систем подразделяются на конечные, итерационные и вероятностные. Мы будем рассматривать в нескольких последующих лекциях только первые две группы, конечные и итерационные. Конечным или прямым методом решения линейных систем называется такой метод, который дает возможность найти решение за конечное число арифметических операций при любых входных данных. Итерационные методы основаны на построение последователь-ности, сходящейся к искомому решению. 3.6 Контрольные вопросы 1. Дать определение понятию нормы вектора и нормы матрицы. 2. Сформулировать понятие числа обусловленности. 3. Что означает система плохо обусловленная? 4. Сформулировать определение фундаментальной системы решений. 5. Дать геометрическую интерпретацию плохой и хорошей обусловленности. 6. Выписать условие невырожденности заданной системы уравнений. 7. Сформулировать определение конечного метода. 8. Выписать правило Крамера для решения линейной системы. Лекция 4 Конечные методы решения систем линейных алгебраических уравнений. Метод исключения Гаусса План: 4.1 Вводные замечания 4.2 Метод исключения Гаусса 4.3 Вычисление определителя матрицы 4.4 Обращение матриц 4.5 Вычислительная сложность схемы Гаусса 4.5 Заключение 4.6 Контрольные вопросы 4.1 Вводные замечания Конечные методы состоят в преобразовании исходной системы уравнений в такую эквивалентную систему, для которой матрица легко обращается и, тем самым, легко находится решение системы. Рассмотрим метод исключения Гаусса, который исторически является одним из первых конечных методов. 4.2 Метод исключения Гаусса Метод исключения Гаусса для решения линейных систем уравнений является основой многих вычислительных схем и сводится к преобразованию исходной системы к равносильной системе с верхнетреугольной матрицей. Преобразования исходной системы осуществляются с помощью линейных операций. Пусть дана система линейных алгебраических уравнений Ax = b, (4.1) где А = (аij) - неособенная размерности nn матрица, detA  0, x = (x1, ..., xn)T и b =(b1, ..., bn)T - векторы-столбцы. Вычислительная схема метода состоит из двух этапов. Первый этап заключается в приведении системы к верхнетреугольному виду. Этот этап называется прямым ходом. Второй этап - определение неизвестных x1, x2, ..., xn называется обратным ходом. Прямой ход состоит в последовательном исключении (обнулении) коэффициентов при неизвестных x1, x2, ..., xn-1, начиная с первого столбца. Исключение коэффициентов при xk - ой переменной называется k -ым циклом метода исключения Гаусса. Прямой ход реализуется по следующим формулам ( индекс k в круглых скобках означает номер цикла (номер столбца) ). Умножение k - ой строки на число (4.2) Вычитание k - ой строки из m - ой строки (4.3) (4.4) В результате n-1-го цикла получается система с верхнетреугольной матрицей (4.5) или Элементы , стоящие на главной диагонали получен- ной системы с верхнетреугольной матрицей (4.5), называются ведущими элементами. Обратный ход - вычисление неизвестных xn, xn-1, ..., x1 реализуется по следующим формулам, начиная с последнего уравнения системы (4.5) (4.6) Формулами (4.2) - (4.6) определяется поэлементная форма записи основной вычислительной схемы метода исключения Гаусса. Необходи- мым и достаточным условием применимости метода является неравенство нулю всех ведущих элементов, 0. 4.3 Вычисление определителя матрицы Преобразования выполняемые по формулам (4.2)- (4.4) не изменяют определителя матрицы A исходной системы (4.1). Как известно определитель треугольной матрицы равен произведению ее элементов, стоящих на главной диагонали, поэтому из представления (4.5) следует, что определитель матрицы A равен произведению всех ведущих элементов в методе Гаусса, именно, (4.7) 4.4 Обращение матриц Применение метода исключения Гаусса для нахождения матрицы, обратной к исходной матрице A, основывается на том, что обратная матрица A-1 является единственным решением матричного уравнения AX=E, (4.8) где E - единичная матрица, X - n  n - матрица. Если обозначить A-1 = X, то матрица X с элементами xij становится искомой. Перемножение матриц A и X приводит к n системам уравнений относительно n2 неизвестных xij (4.9) где Выражение (4.9) определяет систему уравнений для вычисления обратной матрицы в поэлементной форме. Если обозначить столбцы матрицы X, как векторы xi = (x1i, x2i, ..., xni)T, i =1, 2, ..., n, а столбцы единичной матрицы E вектор -столбцами e1, e2, ..., en , то выражение (4.9) представляется в виде набора из n не связанных между собой систем матричных уравнений Ax1 = e1 ; Ax2 = e2 ; ... ; Axn = en, (4.10) с одной и той же матрицей A и правыми частями, отличающимися только расположением 1 в элементах векторов ei. Для решения каждого из уравнений (4.10) применяется вычислительная схема (4.2)-(4.6), а получаемые значения xi будут составлять столбцы обратной матрицы A-1 . 4.5 Вычислительная сложность схемы Гаусса Оценка числа арифметических операций в методе Гаусса проводится непосредственным подсчетом в каждом цикле исключения. Как следует из формул (4.2) - (4.4) в каждом k - ом цикле исключения выполняются последовательно операции: 1. Деления - вычисление множителя dmk = amk/akk , m = k+1, ..., n . 2. Исключения (вычитание и умножение)- amp = amp - dmk akp , p= k+1, ..., n. 3. Вычитания и умножения ( для правой части ) - bm = bm - dmk bk . Число делений в первом цикле составит (n - 1), на втором - (n - 2), и так далее, на последнем цикле одно деление. Таким образом, общее число операций деления для прямого хода составит Операция исключение состоит из умножений на число dmk и поэлементных вычитаний. Число умножений на первом цикле (k = 1) составляет n(n -1), на втором - (n -1)(n -2), и 1.2 умножений при k = n-1, а общее число умножений составит или Аналогичное число составляют и операции вычитания. Таким образом, число операций для исключения равно . Для вычисления правых частей потребуется - операций умножения и столько же операций вычитаний, а общее число операций для вычисления правой части составит Общая трудоемкость прямого хода будет равна . (или порядка арифметических операций). Трудоемкость обратного хода. Обратный ход состоит из операций деления, умножения и вычитания. Для обратного хода согласно формулы(2.6) потребуется - операций. (или порядка О(n2) арифметических операций). Общая трудоемкость метода исключения Гаусса составит . Таким образом, вычислительная трудоемкость метода Гаусса описывается полиномами от n степени не выше третьей с коэффициентом 2/3 при старшем члене. Рассмотрим примеры на применение схемы Гаусса для решения линейной системы и обращения матрицы. Пример 1. Дана система третьего порядка Требуется: а) найти решение данной системы методом исключения Гаусса и вычислить определитель матрицы исходной системы; б) выписать метод исключения в матрично - векторной форме. Решение. а) В данном случае количество циклов равно двум, k =2. Реализация прямого хода осуществляется по формулам (2.2) - (2.4). Для первого столбца, т.е. исключение коэффициентов первого столбца, ниже диагонального элемента, k = 1, имеем , , После первого цикла система уравнений запишется в виде . Далее, цикл по второму столбцу, k = 2, т.е. исключение коэффициен- тов второго столбца, ниже диагонального элемента, имеем После второго цикла получаем систему с верхнетреугольной матрицей . Определитель вычислим с помощью формулы (4.7) и полученной верхнетреугольной матрицы detA= (-5)( -3. 2) (-6. 5) = -104. Вычислим абсолютную и относительную погрешность вычисления определителя, имеем:  =  -104 + 104  = 0;  = 0/  -104  = 0. Для нахождения решения, т.е. вычисление значений x1, x2, x3 исполь- зуем формулы (2.6) обратного хода, получим x3 = (-6. 5)/ (-6. 5) = 1, x2 = (0.2 - 3.4·1)/(-3. 2) = 1, x1 = (3 - 1·1 -7·1)/(-5) = 1. Вычислим вектор - невязки r , имеем: r = Ax - b  . Таким образом, решение 1, 1, 1 является точным. б) Для записи метода исключения в матрично-векторной форме введем матрицы L(1), L(2) , . Процесс прямого хода запишется в виде , ; , . Система для получения вектора - решения принимает вид A(3) x = b(3)  . Пример 2. Для матрицы А найти обратную А-1 А =. Решение. Составим систему из трех матричных уравнений (4.10) , . Решая каждое из этих уравнений методом Гаусса получим столбцы обратной матрицы А-1. Рассмотрим подробно получение первого столбца матрицы А-1. Применяя формулы (4.2) -(4.5) и учитывая, что k =1, 2, имеем, цикл по первому столбцу, k =1, , , , цикл по второму столбцу, k = 2, , В результате получаем систему для вычисления значений x11 , x21 , x31 . Применяя формулы обратного хода (4.6), получим x31 = 1 / 1 =1, x21 = -2 - 2·(-2) ·1 / (-1) =0, x11 = (1- 2·0 - 3·1)/(-1) = -2. Таким образом, первый столбец обратной матрицы будет иметь вид . Для второго и третьего столбцов, проводим аналогичные вычисления, получим , . В итоге, обратная матрица имеет вид А-1=. 4.6 Контрольные вопросы 1. Дать словесное описание вычислительной схемы Гаусса. 2. Сформулировать необходимые условия применения схемы. 3. Выписать формулы прямого хода схемы. 4. Сформулировать общий подход к определению вычислительной сложности. 5. Показать использование обратного хода для обращения верхнетреугольной матрицы. 6. Как измениться вычислительная схема, если матрица и вектор будут размещаться в одном двумерном массиве? 7. Показать использование обратного хода для обращения нижнетреугольной матрицы. 8. Выписать формулы обратного хода схемы. Лекция 5 Метод исключения Жордана – Гаусса План: 5.1 Вводные замечания 5.2 Метод исключения Жордана -Гаусса 5.3 Вычисление определителя матрицы 5.4 Обращение матриц 5.5 Заключение 5.6 Контрольные вопросы 5.1 Вводные замечания Метод исключения Жордана -Гаусса или полного исключения для решения систем линейных алгебраических уравнений неизвестных сводится к преобразованию исходной системы к системе с единичной или диагональной матрицей. 5.2 Метод исключения Жордана –Гаусса Пусть дана система линейных алгебраических уравнений Ax = b, (5.1) где А = (аij) – неособенная размерности nn матрица, detA  0, x = (x1, …, xn)T и b =(b1, …, bn)T – векторы-столбцы. Вычислительная схема метода состоит из n циклов, в каждом из которых последовательно с помощью k-ой строки исключаются элементы при xk –ой неизвестной в каждой из n-1 –ой строки, кроме k –ой. Схема реализуется по следующим формулам: деление k – ой строки на число akk (5.2) вычитание k – ой строки из m – ой строки (5.3) (5.4) В результате n циклов получается система с единичной матрицей (5.5) или . Формулами (5.2) – (5.5) определяется поэлементная форма записи основной вычислительной схемы метода исключения Жордана –Гаусса. Необходимым и достаточным условием применимости метода является неравенство нулю всех ведущих элементов, 0. Другая форма метода исключения Жордана –Гаусса, называемая матрично-векторной, получается из (5.2) – (5.5), если ввести набор матриц L( k) вида , с элементами = - / , i = 1, …, n; k = 1, 2, …, n в k-ом столбце. Эти матрицы отличаются от единичной лишь k-м столбцом. Процесс же преобразований сводится к вычислению последовательных произведений A( k) = L( k)A( k –1); b( k) = L( k)b( k –1), k = 1, …, n, (5.6) и начинается с матрицы A(0) равной исходной матрице A, т.е. A(0) = A, вектора b(0), b(0) = b и матрицы L(1) вида , где = - /, i = 2, …, n. В результате n циклов получается система с единичной матрицей Е и вектором b(n), именно, , или Е x = b(n). (5.7) 5.3 Вычисление определителя матрицы Для метода Жордана-Гаусса вычисление определителя сводится к вычислению произведения ведущих элементов , т.е. . (5.8) 5.4 Обращение матриц Для вычисления обратной матрицы совершаются одинаковые преобразования Жордана -Гаусса в равенстве AА-1 = E одновременно и над исходной матрицей А и над матрицей Е до тех пор, пока А не превратится в единичную Е, а единичная в обратную А-1 При практических вычислениях к матрице А дописывается единичная матрица, затем проводятся преобразования исключения Жордана - Гаусса. Рассмотрим примеры на нахождении решения системы линейных уравнений и обращение матрицы методом исключения Жордана-Гаусса. Пример. Дана система третьего порядка . Требуется: а) найти решение данной систем методом исключения Жордана-Гаусса, получить обратную матрицу и вычислить определитель; б) выписать метод исключения Жордана-Гаусса в матрично -вектор- ной форме. Решение. а) В данном случае количество циклов равно трем. Для получения обратной матрицы одновременно реализуем преобразования и над единичной матрицей. Введем для удобства верхние индексы. Исключим коэффициенты при x1 во второй и третьей строках, для этого вычислим коэффициенты и вычтем первую строку, умноженную на , из второй, получим Аналогично для единичной матрицы: вычитание первой строки, умноженной на , из второй: , вычитание первой строки, умноженной на , из третьей: Аналогично для третьей строки единичной матрицы: . После первого цикла система и единичная матрица имеют вид и . Далее, исключим коэффициенты при x2 в первой и третьей строках. Для этого вычислим коэффициенты , и вычтем вторую строку, умноженную на , из первой, получим , и вторую строку, умноженную на , из третьей . Для единичной матрицы, имеем . Вычитание второй строки, умноженной на , из первой , и второй строки, умноженной на , из третьей . После второго цикла система и единичная матрица имеют вид и . Исключение коэффициентов при x3 в первой и во второй строках: , , Для единичной матрицы, имеем: . Вычитание третьей строки, умноженной на , из первой: , и третьей строки, умноженной на , из второй: . В результате трех циклов получаем решение и обратную матрицу и А -1 = . Для вычисления определителя матрицы применяем формулу (5.8), получим detA =. 5.6 Контрольные вопросы 1. Дать словесное описание вычислительной схемы Жордана -Гаусса. 2. Сформулировать необходимые условия применения схемы. 3. Выписать формулы схемы Жордана -Гаусса. 4. Сформулировать общий подход к определению вычислительной сложности. 5. Определить вычислительную сложность схемы. 6. Как измениться вычислительная схема, если матрица и вектор будут размещаться в одном двумерном массиве? Лекция 6 Компактная схема исключения План: 6.1 Вводные замечания 6.2 Компактная схема исключения 6.3 Вычисление определителя матрицы 6.4 Заключение 6.5 Контрольные вопросы 6.1 Вводные замечания Компактная схема исключения (LU –разложение или схема Холецкого) для решения систем линейных алгебраических уравнений сводится к преобразованию исходной системы к равносильной системе с верхнетреугольной матрицей. При этом матрица исходной системы представляется в виде произведения двух треугольных матриц. 6.2 Компактная схема исключения Компактная схема основана на следующей теореме, которую сформулируем без доказательства Теорема. Пусть квадратная матрица A порядка n имеет ненулевые главные миноры Ak, detAk 0, k =1, 2, ..., n -1. Тогда существуют единственные нижнетреугольная матрица с единичными диагональными элементами L = ( lij), l11= l22= ...= lnn= 1 и верхнетреугольная матрица U, такие, что A = LU. Кроме того, detA = u11 u22 ...  unn, где uii - диагональные элементы матрицы U. Пусть дана система линейных алгебраических уравнений Ax = b, (6.1) где А = (аij) - неособенная размерности nn симметрическая матрица, т. е. A = AT= (aji ), detA  0, x = (x1, ..., xn)T и b =(b1, ..., bn)T - вектор -столбцы. Согласно приведенной теореме матрица системы (6.1) может быть представлена в виде произведения двух треугольных матриц A = LU, (6.2) где L - нижнетреугольная матрица, U - верхнетреугольная матрица. Представление (6.2) дает возможность записать систему (6.1) в виде двух систем уравнений, именно, Ax = b  LUx=b  (6.3) тогда решение системы (6.1) сводится к последовательному решению двух систем с треугольными матрицами L и U. Элементы lij матрицы L и элементы uij матрицы U определяются из следующих уравнений , которые получаются из перемножения матриц L и U в равенстве (6.2). В поэлементной форме (6.3) записывается в виде . Вычислительная схема метода состоит из двух этапов. На первом этапе, называемом прямым ходом, определяются элементы lij матрицы L и элементы uij матрицы U, затем неизвестные z1, z2, ..., zn. Второй этап, называемый обратным ходом, состоит в нахождении неизвестных xn, xn - 1, ..., x1. Прямой ход реализуется по следующим формулам (6.4) (6.5) (6.6) (6.7) Для вычисления z1, z2, ..., zn справедливы формулы (6.8) Обратный ход аналогичен обратному ходу метода Гаусса, именно, (6.9) Формулами (6.4) - (6.9) определяется поэлементная форма записи компактной схемы исключения. Условием применимости метода является отличие от нуля элементов uii, uii 0, i = 1, 2, ..., n. 6.3 Вычисление определителя матрицы Определитель произведения квадратных матриц равен произведению определителей матриц, поэтому из равенства (6.2) следует, что определитель матрицы A равен произведению определителей треугольных матриц L и U, именно, (6.10) На рисунке 1 представлен один из возможных алгоритмов, реализующий нахождение решения системы линейных алгебраических уравнений с помощью компактной схемы исключения. Алгоритм компактной схемы исключения Прямой ход метода -формулы 6.4 - 6.9 вычисление элементов первой строки матриц U и L - формула 6.4 Для j =1 до n выполнять u1, j = a1, j; l ij =a1, j / u11 конец цикла по j вычисление остальных элементов матриц U и L - формулы 6.5 - 6.7 Для i = 2 до n выполнять цикл для вычисления внутренних слагаемых формулы 6.5 Для k = 1 до i - 1 выполнять s = s + li,k·uk, i конец цикла по k ui i = ( ai, i - s ) диагональный элемент матрицы U вычисление внедиагональных элементов матрицы U Для j = i +1 до n выполнять Для k = 1 до i - 1 выполнять вычисления внутр. слагаемых формулы 6.6 s = s + li,k ·uk, j конец цикла по k ui, j = ( ai, j - s) конец цикла по j вычисление внедиагональных элементов матрицы L Для j = i +1 до n выполнять Для k = 1 до j - 1 выполнять вычисление внутр. слагаемых формулы 6.7 s = s + li,k ·uk, j конец цикла по k lj, i = ( aj, i - s)/uii конец цикла по j конец цикла по i z1 = b1 / u11 начало вычисления по формулам 6.8 Для i = 2 до n выполнять Для k = 1 до i -1 выполнять цикл для вычисления формулы 6.8 s = s + li,k ·zk конец цикла по k zi = ( bi - s) / uii конец цикла по i xn = zn / unn начало обратного хода метода - формулы 6.9 Для i = n - 1 до 1 выполнять Для k = i + 1 до n выполнять цикл для вычисления формулы 6.9 s = s + uik·xk конец цикла по k xi = ( zi - s) / uii конец цикла по i Рисунок 1- Алгоритм компактной схемы исключения Замечания к алгоритму. Алгоритм реализован непосредственно по формулам (6.4) - (6.9). Матрицы U и L не сохраняются. Рассмотрим примеры на нахождении решения системы линейных уравнений и вычисление определителя. Пример. Дана система третьего порядка Требуется найти решение данной систем компактной схемой исключения и вычислить определитель матрицы исходной системы. Решение. Первый этап - нахождение элементов матриц L и U, имеем: . Далее, вычислим элементы получим: . Разложение (6.2) A=LU запишется в виде . Выпишем систему (первая система из (4.3)) для определения значений z1, z2, z3 и найдем эти значения, , тогда Для вычисления значений x1, x2, x3 запишем систему (вторая система из (6.3)) и используем формулы (6.9) обратного хода, получим и . Для вычисления определителя используем формулу (6.10), получим detA = detLdetU = 6.5 Контрольные вопросы 1. Дать словесное описание компактной схемы. 2. Сформулировать необходимые условия применения схемы. 3. Выписать формулы схемы. 4. Сформулировать общий подход к определению вычислительной сложности. 5. Определить вычислительную сложность схемы. 6. Как измениться вычислительная схема, если матрица и вектор будут размещаться в одном двумерном массиве? Лекция 7 Итерационные методы решения систем линейных алгебраических уравнений. Метод простых итераций План: 7.1 Вводные замечания 7.2 Общий подход к построению итерационных методов 7.3. Метод простых итераций 7.4 Условия сходимости метода 7.5 Оценка погрешности приближений 7.6 Заключение 7.7 Контрольные вопросы 7.1 Вводные замечания Итерационные методы дают решение системы уравнений в виде предела последовательности некоторых векторов. Построение этих векторов реализуется с помощью процесса, называемого процессом итераций. Введем основные понятия связанные с итерационными методами. В общем случае под итерационным алгоритмом понимается алгоритм, реализующий на каждом множестве W последовательность отображений Fk: W W, при помощи которых по некоторой начальной точке w (0)  W вычисляется последовательность точек по правилу w( k+1) = Fk w( k), k = 0, 1, ... . (7.1) При этом последовательность w(k) называется итерационной последовательностью, а правило (1) - итерацией. Наиболее разработаны и исследованы итерационные алгоритмы для решения систем линейных алгебраических уравнений Ax = b. Эти алгорит- мы и соответствующие им методы, подразделяются на линейные и нелинейные. К линейным относятся метод простых итераций (метод Якоби), метод итераций Гаусса-Зейделя, метод релаксаций, к нелинейным -метод градиента, метод минимальных невязок и другие. 7.2 Общий подход к построению итерационных методов Общий подход к построению отображения Fk для линейной системы Ax=b состоит в том, что к обеим частям системы добавляется вектор x, разность b - Ax умножается на некоторую неособенную матрицу H. В результате получается уравнение, называемое приведенным уравнением x= x + H(b - Ax), или x = (E - HA)x + Hb, которое может быть записано, в общем случае, как x = Cx + Hb. Матрица H может быть как постоянной, в этом случае, итерационный процесс x(k + 1) = C x(k) + Hb называется стационарным, так и изменяющейся от итерации к итерации, при этом итерационный процесс x(k + 1) = Ckx(k) + Hk b называется нестационарным. Отличие итерационных методов определяется выбором матрицы H и построением матрицы С. Итерационный процесс называется сходящимся при некотором начальном приближении x(0) к решению x*, если x(k)  x* при k  . В общем случае, для сходимости итерационного процесса вида x(k + 1) = Ckx(k) + Hk b необходимо и достаточно выполнение условия при k  . Для стационарного итерационного процесса x(k + 1) = C x(k) + Hb необходимо и достаточно , чтобы спектральный радиус матрицы С был меньше единицы, т.е.(С) < 1. Достаточным условия сходимости стационарного процесса является выполнение неравенства Для характеристики сходимости итерационного процесса вводится понятие средней скорости сходимости и асимптотической скорости сходимости. Так для стационарного итерационного процесса средняя скорость сходимости за k итераций при выполнении условия определяется как Rc = -ln/ k, а асимптотическая скорость как Rac = -ln(С). Для оценки качества приближений итерационного процесса и его точности используются понятия вектора ошибки и вектора невязки. Вектор ошибки определяется как (k) = x* - x(k) и характеризует величину отклонения очередного приближения x(k) от искомого точного вектора-решения x*. Вектор невязки rk = Ax(k) - b -характеризует качество приближенного решения x(k) системы уравнений Ax = b. Между этими характеристиками имеются следующие соотношения A = r,  = A-1r, (k+1) = Ck(k) = (E -HkA) (k), r(k+1) = ACk A-1r(k) = A(E -HkA)A-1r(k). С помощью вектора ошибок определяется средняя скорость сходи мости нестационарного итерационного процесса. Средней скоростью сходимости для нестационароного итерацион-ного процесса является величина Rk = -ln/ k, если При этом асимптотической скоростью сходимости называется Rac = limRk при k   . Особенностью итерационных методов является их самоисправля- емость и простота реализации на вычислительной машине. Если в конеч- ных методах ошибка в вычислениях, когда она не устраняется случайно, ведет к ошибке в результате, то в случае сходящегося итерационного процесса ошибка в каком-то приближении исправляется в последующих вычислениях. Рассмотрение итерационных методов начнем с метода простых итераций. 7.3 Метод простых итераций Метод простых итераций является одним из общих и простых методов приближенного решения системы линейных алгебраических уравнений Ax = b матрица которой имеет отличные от нуля диагональные элементы. Пусть дана система линейных алгебраических уравнений Ax = b, (7.1) где А = (аij) - неособенная размерности nn матрица, detA  0, x = (x1, ..., xn)T и b =(b1, ..., bn)T - векторы-столбцы. Метод простых итераций требует приведения системы (7.1) к виду x = Cx + D , (7.2) где C = - G-1(A - G) = ( E - G-1A), D = G-1b, G = diag(aii), i = 1, 2, ..., n, - диагональная неособенная матрица, detG 0, G-1 = diag(1/aii) - обратная к G матрица, E -единичная матрица, с единицами на главной диагонали, а остальные элементы нулевые. Для применения метода итераций к системе (2.2) задается вектор начального приближения x(0). Затем, строятся последователь-ные приближения x(1) = Cx(0) + D, x(2) = Cx(1) + D, x(3) = Cx(2) + D, ..., x(k+1) = Cx(k) + D, ... . В общем случае, любое (k +1)-е приближение вычисляется по формуле x (k +1) = Cx(k) + D, k = 0,1, 2,...- номера итераций. (7.3) Используя обозначение для матрицы C из (7.2) формула (7.3) записывается в виде x( k+1) = - G-1(A - G) x( k) + G-1b, (7.4) или x( k+1) = ( E - G-1A) x( k) + G-1b, (7.5) или x( k+1) = x( k) + G-1(b - A x( k) ), k = 0, 1, 2, ... . (7.6) x(1) = Cx(0) + D, x(2) = Cx(1) + D, x(3) = Cx(2) + D, ..., x(k+1) = Cx(k) + D, ... . Формулами (7.3) - (7.6) определяется метод простых итераций. 7.4 Условия сходимости метода Итерационный процесс (7.3) - (7.6) называется сходящимся, если последовательность x(k) сходится к x* при любом начальном значении x(0) . Итерационный процесс называется стационарным, если матрица С не зависит от номера k шага итерационного процесса, в противном случае процесс называется нестационарным. Необходимые и достаточные условия сходимости последователь-ности приближений x(k) к решению x* системы (7.2) определяется следующей общей для итерационных методов теоремой. Сформулируем ее без доказательства. Теорема. Последовательность приближений x(k) сходится при любом начальном приближении x(0) к решению x* системы (2.2) тогда и только тогда, когда все собственные значения матрицы C по модулю меньше единицы, т.е. i < 1 ( i =1, 2, ..., n ) Комментарий к теореме Теорема гарантирует получения единственного решения, но при этом требуется знание всех собственных чисел матрицы C. Нахождение же собственных чисел является непростой и самостоятельной задачей. Поэтому при практических вычислениях применяются более простые достаточные условия сходимости. Для метода простых итераций необходимые и достаточные условия устанавливает следующая теорема. Теорема 1. Последовательность приближений x(k) в методе простых сходится при любом начальном приближении x(0) к решению x* системы (2.2) тогда и только тогда, когда все корни уравнения по модулю меньше единицы. Комментарий к теореме Теорема является развернутой записью условий общей теоремы через элементы матрицы А системы (7.1). Теорема 2. Для сходимости последовательности приближений x(k) к единственному решению x* системы (7.2) при любом начальном приближении x(0), достаточно, чтобы какая-либо норма матрицы E -G-1A была по модулю меньше единицы, т.е. Комментарий к теореме. Теорема дает достаточный признак сходимости. Для его применения требуется вычислить какую-либо норму матрицы E -G-1A и проверить выполнение неравенства, указанного в теореме. Теорема 3. Если в матрице А системы (2.1) выполняются неравенст-ва или , то метод простых итераций сходится. Комментарий к теореме. Теорема утверждает, что диагональное преобладание в элементах исходной матрицы А системы (7.1) достаточно для сходимости процесса итераций. 7.5 Оценка погрешности приближений Для метода простых итераций оценка погрешности приближений определяется следующей теоремой. Теорема 4. Если какая-либо норма матрицы C, согласованная с рассматриваемой нормой вектора x, меньше единицы, то справедлива следующая оценка погрешности приближений . (7.7) Комментарий к теореме. Теорема дает количественную оценку нормы вектора ошибки. Здесь под матрицей С подразумевается матрица E -G-1A . Следствие 1. В условиях теоремы 4 и формулы (7.7) имеет место следующая оценка . (7.8) Комментарий к следствию. С помощью неравенства (2.8) можно определить число (априорную оценку) итераций необходимых для получения решения с заданной точностью. Следствие 2. Если в формуле (7.7) или (7.8) взять за начальный вектор x(0), вектор свободных коэффициентов D, то справедлива оценка (7.9) Комментарий к следствию. Неравенство (7.9) аналогично (7.8), но получаемая с помощью него априорная оценка более грубая. Один из возможных алгоритмов вычисления по методу простых итераций при выполнении достаточных условий сходимости показан на рисунке. Рассмотрим примеры на использование метода простых итераций для решения системы алгебраических уравнений, применение условий теоремы 2, неравенства (7.9) и условий теоремы 1. Пример 1. Дана система третьего порядка 2x1 + 0.25 x2 + 0.4x3 = 4 x1 + 3 x2 + 0.3x3 = 4.5 0.08x1 + 4.02 x2 + 5x3 = 9.1 . Требуется: а) привести систему к виду (7.2); б) показать, что процесс итераций сходится; в) получить число итераций, необходимых для получения решения с точностью 10-5; г) выписать формулу итерационного процесса. Решение. Используя формулу (7.5) приведем систему к виду x=Cx+D, получим систему x1 = 2 - 0.125x2 - 0.2x3 x2 = 1.5 - 0.3x1 - 0.1x3 x3 = 1.82 - 0.016x1 - 0.804x2 , которая в матричной форме имеет вид Для проверки сходимости используем достаточный признак (теорема 2). В качестве нормы возьмем -норму, получим = max{0.325; 0.4; 0.82 }= 0.82 < 1, условия теоремы 2 выполняются. Оценку числа итераций для получения решения с точностью 10-5, проведем используя неравенство (7.9). За начальное приближение x(0) возьмем вектор D, вычислим норму = max{2; 1.5; 1.82 } = 2, имеем < 10-5 или (0.82)k+1 < 0.0910-5 , отсюда следует оценка для числа k ( k + 1 ) > ln(910-7)/ln(0.82) 70.45. Для записи итерационного процесса используем формулу (7.6), получим представление итерационного процесса в координатной форме В результате расчета по программе, алгоритм которой показан на рисунке, получается сходимость к решению за 23 итерации, при началь-ном векторе x1(0) =2, x2(0) =1.5, x3(0) =1.82 и решение x1 1.671076, x2  0.830413, x3  1.125609 . Пример 2. Дана система третьего порядка 2x1 + x2 = 0.9 x1 + 3 x2 + x3 = 2.1 x2 + 4x3 = 2.1 . Требуется: а) проверить выполнение условий теоремы 1; б) выписать формулу итерационного процесса. Решение. Согласно теоремы 1 выпишем уравнение для нахождения корней, имеем , отсюда получаем 1, 2 = 0, 3 = 2/3, т.е.1,2,3 < 1. Таким образом условия теоремы 1 выполняются. Используя формулу (2.6), итерационный процесс записывается в виде: Алгоритм вычисления по методу простых итераций ввод исходных данных:матрица A,вектор b, вектор x0 и точность e повторять для i=1 до N выполнять начало для j=1 до N выполнять s= s +( A[i,j]*X0[j]) -вычисление внутренней суммы формулы (2.6) конец по j X1[i] =X0[i] - (s - B[i])/A[i,i] - вычисление очередной итерации (2.6) конец по i проверка условия окончания вычислений здесь используется m - норма max:=abs(X1[1] - X0[1]) для i=2 до N выполнять если abs(X1[i] - X0[i]) > max то max=abs(X1[i] - X0[i]) конец по i сохранение вычисленных значений для i=1 до N выполнять X0[i]= X1[i]; пока max < e конец цикла Рисунок - Алгоритм метода простых итераций Замечания к алгоритму. Как видно из алгоритма для реализации одного шага требуется сохранение массива X0 до тех пор, пока не сформируется полностью другой массив, массив X1- результат текущего итерационного шага. По этой причине метод простых итераций называется еще и как метод одновременных сдвигов (смещений). 7.7 Контрольные вопросы 1. Описать общую схему метода простых итераций. 2. Сформулировать достаточные условия сходимости метода простых итераций. 3. Сформулировать необходимые и достаточные условия сходимости мето да простых итераций. 4. Выписать неравенство для оценки сходимости последовательных приближений. 5. Дать основной ход вычислений априорного числа итераций для любой выбранной нормы. 6. Выписать матричную форму записи метода итераций. 7. Выписать формулу записи метода итераций удобную для программирования. Лекция 8 Метод итераций Гаусса - Зейделя План: 8.1 Вводные замечания 8.2 Метод итераций Гаусса - Зейделя 8.3 Условия сходимости метода 8.4 Оценка погрешности приближений 8.5 Заключение 8.6 Контрольные вопросы 8.1 Вводные замечания Метод Гаусса-Зейделя - итерационный метод решения системы линейных алгебраических уравнений Ax = b, матрица А которой с ненулевыми диагональными элементами. Метод Гаусса-Зейделя является некоторым видоизменением метода простых итераций. Рассмотрим этот итерационный метод, используя такой же подход как и в предыдущей лекции. 8.2 Метод итераций Гаусса – Зейделя Пусть дана система линейных алгебраических уравнений Ax = b, (8.1) где А = (аij) - неособенная размерности nn матрица, detA  0, x = (x1, ..., xn)T и b =(b1, ..., bn)T - векторы-столбцы и пусть система (8.1) приведена к виду x = Cx+ D. (8.2) Метод итераций Гаусса-Зейделя представляет собой такое преобразование метода простых итераций x (k +1) = Cx(k) + D для решения системы (8.1), приведенной к виду (8.2), при котором для вычисления i-ого элемента (k + 1) - го приближения используются уже найденные на этом (k + 1)-м шаге значения первых i -1 элементов. При этом итерационный процесс строится так же как и в методе простых итераций. В поэлементной форме метод Гаусса-Зейделя записывается в виде следующей системы равенств или x1(k+1) = (b1 - a1jxj(k) )/a11 , xi(k+1) = (bi -aijxj(k+1) -aijxj(k) )/aii, i =2, 3, ..., n; k = 0,1, ..., где первый элемент x1(k+1) вычисляется по методу простых итераций, aij и bi элементы исходной матрицы A и вектора b системы (8.1), x(0) - вектор начального приближения. Система равенств (8.3) определяет итерационный метод Гаусса-Зейделя, записанный по элементам матрицы A и вектора b исходной системы (8.1). Другая форма записи метода Гаусса-Зейделя называется матрично-векторной и получается из представления матрицы C и вектора D системы (8.2) в виде C = - ( L + G ) -1 R, D = (L+ G ) -1b, где L = (aij) -нижняя треугольная матрица, R = (aij), - верхняя треугольная матрица, G = (aii) - ненулевая диагональная матрица, ( L + G ) - нижняя треугольная часть матрицы A, обратная к ней существует в силу предположения о том, что G ненулевая диагональ, det(L+G)  0, тогда итерационный процесс x(k+1) = Cx(k) + D записывается следующим образом: x( k+1) = - ( L + G ) -1 R x( k) + (L+ G )-1b, k = 0, 1, 2, .... (8. 4) Если ввести единичную матрицу E, с единицами на главной диагонали, а остальными элементами равными нулю, то (8.4) примет вид x(k+1)=(E - ( L + G ) -1A )x(k ) + ( L + G ) -1b, (8.5) или x( k+1) = x(k ) + ( L + G ) -1 (b - A x(k ) ), k = 0,1, 2,.... (8.6) Формулами (8.4) - (8.6) определяется итерационный метод Гаусса-Зейделя в матрично-векторной форме. 8.3 Условия сходимости метода Сходимость метода итераций Гаусса-Зейделя определяется следующими теоремами. Необходимые и достаточные условия сходимости последователь-ности приближений x(k) к решению x* системы (8.2) устанавливает теорема. Теорема 1. Последовательность приближений x(k) в методе Гаусса-Зейделя сходится при любом начальном приближении x(0) к решению x* системы (8.1) тогда и только тогда, когда все корни уравнения по модулю меньше единицы.. Комментарий к теореме. Теорема является прямым следствием общей теоремы о сходимости последовательности приближений, где в качестве матрицы С принимается матрица С = - ( L + G ) -1 R. С увеличением размерности исходной систе- мы уравнений проверка выполнения условий теоремы затруднительна, поэтому при практических вычислениях пользуются условиями на основе оценок норм матрицы С. Теорема 2. Для сходимости последовательности приближений x(k) к единственному решению x* системы (8.2) при любом начальном прибли-жении x(0), достаточно, чтобы какая-либо норма матрицы ( L + G ) -1 R или матрцы ( E - ( L + G ) -1A ) была по модулю меньше единицы, т.е. Комментарий к теореме. Теорема устанавливает достаточный признак сходимости. Для практической проверки выполнения условий требуется представление матрицы (L + G ) -1 R в поэлементной форме исходя из системы равенств (8.3). При этом равномерная (- или m-) норма и l- норма могут вычисляться по формулам: (8.7) . (8.8) Теорема 3. Если в матрице А системы (8.1) выполняются неравенст-ва или , то метод итераций Гаусса-Зейделя сходится. Комментарий к теореме. Теорема устанавливает достаточный признак сходимости, аналогичный методу простых итераций и основанный на диагональном преобладании в элементах матрицы исходной системы. Теорема 4. Если система (8.1)- нормальная, то процесс итераций Гаусса-Зейделя для системы (8.3) сходится при любом начальном приближении x(0) . Комментарий к теореме. Теорема устанавливает еще один достаточный признак сходимости метода итераций Гаусса -Зейделя. 8.4 Оценка погрешности приближений Для метода итераций Гаусса-Зейделя оценка погрешности приближений определяется следующей теоремой. Теорема 5. Если какая-либо норма матрицы C, согласованная с рассматриваемой нормой вектора x, меньше единицы, то справедлива следующая оценка погрешности приближений . (8.9) Комментарий к теореме. Теорема дает количественную оценку величины нормы вектора ошибки метода Гаусса-Зейделя. Под матрицей C понимается матрица ( L + G ) -1 R . Следствие. В условиях теоремы 5 и формулы (8.9) имеют место следующие оценки , (8.10) где q = , вычисляется по формуле (3.7) и , (8.11) где t = , вычисляется по формуле (8.7), . Комментарий к следствию. Неравенства (8.10) и (8.11) дают возможность вычисления числа итераций (априорная оценка) необходимых для получения решения с заданной точностью с применением равномерной нормы и l- нормы. Следствие 2. Если в формуле (8.10) или (8.11) взять за начальный вектор x(0), вектор свободных коэффициентов D, то справедливы оценки (8.12) и . (8.13) Комментарий к следствию. Неравенства (8.12) и (8.13) аналогичны оценкам (8.11), (8.10) и дают возможность определения числа итераций, но при этом, получаемая оцен- ка более грубая. Один из возможных алгоритмов вычисления по методу итераций Гаусса-Зейделя при выполнении достаточных условий сходимости показан на рисунке. Рассмотрим примеры на использование метода итераций Гаусса-Зейделя для решения системы алгебраических уравнений, применение условий теоремы 2, неравенства (8.12) и условий теоремы 1. Пример 1. Дана система третьего порядка 9x1 + 3 x2 + x3 = 6.5 x1 - 5 x2 + 2x3 = -1 2x1 + x2 + 7x3 = 5 . Требуется: а) привести систему к виду (8.3); б) показать, что процесс итераций сходится; в) получить число итераций, необходимых для получения решения с точностью 10-5; г) выписать формулу итерационного процесса. Решение. Используя формулу (8.3) приведем систему к виду или в матрично - векторной записи Для проверки сходимости используем достаточный признак (теорема 2). В качестве нормы возьмем -норму, вычисляя по (8.7), получим = max{} 0.83 < 1, условия теоремы 2 выполняются. Оценку числа итераций для получения решения с точностью 10-5, проведем используя неравенство (8.12). За начальное приближение x(0) возьмем вектор D, вычислим норму = max{6.5; -1; 5 } = 6.5, имеем < 10-5 или (0.83) k+1 < 0.02610-5 , отсюда следует оценка для числа k ( k + 1 ) > ln(0.02610-5) / ln(0.83)  81. Для записи итерационного процесса используем формулу (8.6), получим В результате расчета по программе, алгоритм которой показан на рисунке, получается сходимость к решению за 6 итерации, при начальном векторе x1(0) = 0.1, x2(0) = 0.1, x3(0) = 0.1 и решение x1  0.19979, x2  0.19887, x3  0.18897 . Пример 2. Дана система третьего порядка Требуется: а) проверить выполнение условий теоремы 1; б) выписать формулу итерационного процесса. Решение. Согласно теоремы 1 выпишем уравнение для нахождения корней, имеем , отсюда получаем 1, 2 = 0, 3 = 1/4, т.е.1,2,3 < 1. Таким образом условия теоремы 1 выполняются. Используя формулу (8.6), итерационный процесс записывается в виде Алгоритм вычисления по методу итераций Гаусса-Зейделя ввод исходных данных: матрица A,вектор b, вектор x( 0) и точность e повторять для i =1 до N выполнять начало для j =1 до N выполнять s = s + ( A[i,j]*X0[j]) -вычисление внутренней суммы формулы (3.6) конец по j NX[i] = X0[i] X0[i] =X0[i] + ( B[i]-s) / A[i,i] - вычисление очередной итерации (3.6) конец по i проверка условия окончания вычислений здесь используется m - норма max = abs(X0[1] - NX[1]) для i =2 до N выполнять если abs(NX [i] - X0[i]) > max то max = abs(X0[i] - NX[i]) конец по i сохранение вычисленных значений пока max < e конец цикла Рисунок - Алгоритм метода итераций Гаусса-Зейделя Замечания к алгоритму. Как видно из приведенного алгоритма для реализации метода требуется один и тот же массив X0, который в начале алгоритма заполняется начальными значениями (компонентами начального вектора), а входе вычисления элементы массива X0 замещаются новыми элемента ми. По этой причине метод Гаусса-Зейделя называют еще и как метод последовательных сдвигов (смещений). 8.6 Контрольные вопросы 1. Описать общую схему метода итераций Гаусса - Зейделя. 2. Выписать матричную форму записи метода Гаусса -Зейделя. 3. Сформулировать необходимые и достаточные условия сходимости метода итераций Гаусса - Зейделя. 4. Сформулировать достаточный признак сходимости по одной из приведенных норм. 5. Дать основной ход вычислений теоретического ( априорного ) числа итераций и проверки достаточных условий сходимости. 6. Выписать формулу записи метода итераций удобную для программирования. 7. Выписать оценку между вектором ошибки и двумя последовательными приближениями. 8. Показать связь между вектором ошибки и двумя начальными приближениями. Лекция 10 Основные сведения о приближении функций. Постановка задачи об интерполировании функций алгебраическими многочленами План: 10.1 Основные сведения о приближении функций. Постановка задачи 10.2 Постановка задачи об интерполировании функций алгебраическими много-членами. Алгебраический многочлен Лагранжа 10.3 Оценка погрешности многочлена Лагранжа. Вычислительная сложность 10.4 Заключение 10.5 Контрольные вопросы 10.1 Основные сведения о приближении функций. Постановка задачи Задача о приближении функции состоит в приближенной замене (аппроксимации) данной функции обобщенным полиномом таким образом, чтобы отклонение данной функции от этого полинома было наименьшим. Предположим, что задано некоторое упо­ря­до­чен­­­­­ное множество вещественных абcцисс х1, х2, ..., хn и свя­зан­ное с ним множество вещественных ор­ди­­нат у1, у2, ..., уn. Пусть х1< х2< ...< хn и каждое уi есть некоторое ве­щес­т­венное число, отвечающее хi, которое определяется ма­тематически или в результате каких-либо на­блю­де­ний. Точ­ки (xi,yi) называются узлами ин­тер­поляции. Кривая, которая точно проходит через эти уз­лы, называется интерполяционной кривой. Требуется построить такую функцию F(x) значения которой в заданных узлах интерполяции совпадало бы с значениями у1, у2, ..., уn, т.е. F(xi) = yi, i = 0, 1, ..., n. Для однозначной разрешимости задачи качестве функций используют полиномиальные функции. 10.2 Постановка задачи об интерполировании функций алгебраическими многочленами. Алгебраический многочлен Лагранжа Задача интерполирования данной функции алгебраическими многочленами состоит в построении такого полинома, значения которого совпадают со значениями данной функции в узлах интерполяции. Сформулируем задачу более строго. Постановка задачи. Пусть для функции y = f(x) заданы значения yi = f(xi) в неравноотсоящих (n + 1)-ом узлах интерполяции, т.е. f(x0) = y0, f(x1) = y1, . . ., f(xn) = yn. Требуется поcтроить полином Ln(x) степени не выше n и принимающий в заданных узлах (точках) xi, i = 0, 1, ..., n значения, совпадающие со значениями функции f(x), т.е. Ln(xi ) = yi , i = 0, 1, 2, ..., n . Алгебраический (интерполяционный) полином Лагранжа для произвольно заданных узлов. Теорема 1. Пусть заданы абциссы xi, i= 0, 1,…, n, среди которых нет совпадающих, xi  xj, ij, и заданы значения y0= f(x0), y1= f(x1), . . ., yn= f(xn) функции f(x), yi = f(xi), i= 0, 1,…, n. Тогда существует и притом единственный полином (10.1) степени не выше n, принимающий в заданных xi, i = 0, 1, ..., n заданные значения yi, Ln(xi) = yi, i= 0, 1,…, n. Полином (10.1) называется интерполяционным полиномом Лагранжа для неравноотстоящих узлов, а коэффициенты при yi называются лагранжевыми коэффициентами. Доказательство.(без док-ва) Другая форма записи полинома Лагранжа. Если ввести обозначение wn+1(x) = (x- x0) (x- x1) … (x- xi-1) (x- xi) (x- xi+1) … (x- xn) , тогда коэффициент Лагранжа запишется в виде , где есть производная от wn+1(x) при x = xi, именно . (сформулировать как задачу) Полином Лагранжа (10.1) будет иметь вид (10.2) Интерполяционный многочлен Лагранжа для равноотстоящих узлов. Теорема 2. Пусть заданы равноотстоящие узлы xi= x0 +ih , i= 0, 1,…, n, и заданы значения y0=f(x0), y1= f(x1), . . ., yn= f(xn) функции f(x) в этих узлах. Тогда существует и притом единственный полином , (10.3) степени не выше выше n, принимающий в заданных узлах xk, i = 0, 1,…, n заданные значения yk, Ln(xk)= yk, k= 0, 1,…, n. , Доказательство. Введем переменную , тогда x = x0 + mh Тогда для произвольного xi имеем а разность запишется в виде: wn+1(x) = (x- x0) (x- x1) … (x- xi-1) (x- xi) (x- xi+1) … (x- xn) = hn+1 m(m- 1) (m- 2) … (m- n) и (-1)n –ii!(n- i)! Коэффициент Лагранжа в этом случае имеет вид: , i= 0, 1, …, n, где . Отсюда, . Например, для n = 1 многочлен Лагранжа (10.3) будет иметь вид 10.3 Оценка погрешности многочлена Лагранжа. Вычис-лительная сложность Теорема 3. Для интерполяционного полинома Лагранжа остаточный член (погрешность)) определяется формулой , [a, b] (10.4) где ξ зависит от x и находится внутри [a, b]. Погрешность в случае равноотстоящих узлов. Теорема 4. Для интерполяционного полинома Лагранжа от равноотстоящих узлов остаточный член (погрешность)) определяется формулой , [a, b]. (10.5) где ξ зависит от x и находится внутри [a, b]. Теорема 5. Оценка остаточного члена интерполяционного полинома Лагранжа определяется формулой . (10.6) Доказательство. Если обозначить через Cn+1 = max , то и получим оценку, именно, x [a, b], тогда подставляя Cn+1 в (10.5) получим искомую оценку для остаточного члена. Вычислительная сложность. На вычисление значений функции в одной точке требуется 2n2 +3n операций сложения, 2n2 +n-1 операций умножения и n+1 операций деления. На рисунке представлен один из возможных алгоритмов, реализующий интерполирование полиномом Лагранжа для произвольных узлов (формула (10.2.) ). Алгоритм интерполяции полиномом Лагранжа ввод значения z и число узлов n s = 0; Для k = 0 до n выполнять c =1; Для j = 0 до n выполнять h =xk -xj; //вычисление шага if ( k == j) { h = z -xk; } конец если if (h != 0 ) {c*=(float)(z-x[j])/h;} //вычисление k - го коэффициента else { printf("s=%f\n",f[k]); } //введенное число совпадает с узлом конец цикла по j s+=c*f[k]; //формирование результата конец цикла по k вывод результата Рисунок - Описание алгоритма интерполяции полиномом Лагранжа Комментарий к алгоритму. Алгоритм реализует вычисления непосредственно по формуле (10.3.) для таблично заданной функции. Рассмотрим пример на применение оценки погрешности (10.6). Пример. С какой точностью можно вычислить с помощью интерполяционной формулы Лагранжа для функции , выбрав узлы интерполирования x0=100, x1=121, x2 =144 ? Решение. Вычислим производные до третьего порядка, имеем , , . Отсюда С3= max  = при 100  x  144. На основании формулы (10.6) получим R2 . 10.5 Контрольные вопросы 1. Сформулировать задачу интерполирования функции полиномом Лагранжа. 2. Выписать полином Лагранжа в компактной форме для произвольных узлов. 3. Выписать оценку погрешности формулы Лагранжа для равноотстоящих узлов. 4. Выписать полином Лагранжа для трех узлов. 5. Выписать остаточный член формулы Лагранжа для равноотстоящих узлов и произвольных узлов. 6. Сформулировать теорему о полиноме Лагранжа для произвольной системы узлов. Лекция 11 Конечные и разделенные разности, свойства. Алгебраические полиномы Ньютона План: 11.1 Алгебраические полиномы Ньютона 11.2 Погрешности полиномов Ньютона 11.3 Заключение 11.4 Контрольные вопросы 11.1 Алгебраические полиномы Ньютона В силу единственности многочлена степени n, по­­­­стро­ен­ного по n + 1 значениям функции f(х), из предыдущей лекции, мно­­­гочлен Нью­то­на Рn(х) является раз­но­вид­нос­тью записи ин­тер­по­ля­ци­он­но­го мно­гочлена и по­л­нос­тью сов­па­да­ет с многочленом, по­стро­ен­ным по фор­­муле Ла­гран­жа. Однако с помощью мно­го­чле­­­на Нью­­тона удобнее про­во­дить интерполяцию в на­­ча­­ле и конце таблицы зна­че­ний функции. В на­ча­ле таб­­лицы используют пе­рвую ин­­тер­по­ля­ци­он­ную фор­­му­лу Ньютона, а в кон­це - вто­­рую. Рас­­смотрим один из спо­­со­бов по­стро­ения пер­вой ин­­тер­по­ля­ци­он­ной фор­му­­лы. Пер­вая ин­­тер­по­ля­ци­он­ная фор­му­­ла. Пусть некоторая функция f(х) задана таб­­лич­ны­ми зна­­чениями у0 = f(х0); у1= f(х1); ...; уn = f(хn) в рав­но­от­сто­­ящих узлах интерполяции {х0 , х1 = х0 + h;...; хn = x0+ nh}. Требуется по­стро­ить ин­тер­по­ля­ци­он­ный поли­ном Нью­тона Рn(х) сте­пе­ни n, при котором Рn(х0)  у0; Рn(х1)  у1; ...; Рn(хn)  уn . (11.1) Будем искать полином в виде Pn(x0) = a0 + a1(x - x0) +a2(x - x0) (x - x1) + ... ...+ an(x - x0) (x - x1)...(x - xn), (11.2) где неизвестные коэффициенты аi находятся из простых зависимостей. Для того чтобы най­­­ти а0, по­ло­жим х=х0. Очевидно, что при этом Рn(х0)  у0 а0 [это сле­ду­ет из формулы (11.1) ]. Но так как все члены урав­нения (11.2) , кроме пер­вого, со­де­р­­жат сомножитель (x - x0), сле­до­ва­­тель­но, они все станут равными нулю, и из формулы (11.2) с учетом формулы (11.1) имеем Рn(х0) = а0  у0 . Для нахождения а1, положим х = х1. Повторив все рас­суждения и учитывая, что значение по­ли­нома в ука­зан­ной точке будет тождественно рав­­но у1 (формула (11.1)), пос­­ле подстановки в формулу (11.2) х1 име­ем Pn(x1) = a0 + a1(x1 - x0) = у0 + а1 h = у1 . Все остальные со­мно­жи­те­ли при неизвестных ко­эф­фи­циентах аi будут рав­ны ну­лю. Преобразуя по­след­нее вы­ражение, на­хо­дим а1 как а1 = 1/h, здесь 1 = Рn(х0 + h) - Рn(х0) = у1 - у0 . Для определения а2, положим х = х2 и, рас­­суж­дая аналогично, определим третий ко­эф­фи­циент как a2  = 2 / (2! h2). Подставляя в выражение (11.2) последовательно все хn, при­хо­дим к общей формуле для получения коэф­фи­ци­ен­тов аi: ai = i / ( i! hi), которые подставим в фор­­мулу (11.2) и получим пер­вую ин­­тер­по­ля­ци­он­ную формулу Ньютона: (11.3) В формуле (11.3) обычно выполняют замену пе­ре­мен­ных: q = (х - х0)/h, где h - шаг ин­тер­­по­ли­ро­ва­ния. Тогда пер­вая интер­по­ля­ци­он­ная фор­мула Нью­то­на за­пи­сы­ва­ет­ся в виде: (11.4) При n = 1 получаем формулу линейного ин­тер­по­ли­рования; при n = 2 - параболического ин­тер­по­ли­ро­вания и т.д. Вторая интерполяционная формула Ньюто­на. Вторую формулу Ньюто­на по­­лу­ча­ют, если узлы интерполяции в Рn(х) бе­рут в другом порядке: Pn(x0) = a0 + a1(x - xn) + a2(x - xn) (x - xn-1) + ... ...+ an(x - xn) (x - xn-1) ... (x - x0) . Далее, рассуждая как и в случае пер­вой ин­тер­по­ля­ци­онной формулы, получаем ис­ко­мую форму записи по­ли­нома Ньютона, кото­рая из­вест­на как вторая ин­тер­по­ляционная фор­му­ла: Выполнив подстановку q = (х - хn)/h, получим иную за­пись интерполяционного многочлена Ньютона: (11.5) Первая интерполяционная формула Ньютона ис­поль­зу­ет­­ся для интерполирования в начале отрезка [хi, хi+1] и экстра­по­лирования до первой точки х0, т.е. для ин­тер­по­ли­ро­ва­ния впе­ред и экс­траполирования назад. При таком ин­тер­по­­ли­ро­ва­­нии q = (х - хi)/h > 0. При экс­т­ра­по­лировании на­зад по первой интерпо­ля­ци­­онной фор­муле q = (х - хi)/h < 0. При ин­тер­по­ли­ро­­вании в кон­це таблицы, т.е. при ин­­терпо­ли­ро­ва­нии назад, ког­да шаг ин­тер­по­ляции по­­сто­я­нен, при­меняют вто­рую ин­тер­по­ля­ци­он­ную фор­­му­лу Нью­тона, где q = (х - хi+1) / h < 0. Эту же фор­­­мулу ис­поль­зу­ют при экст­ра­­по­ли­ро­ва­нии впе­ред (в кон­це таблицы), но тогда q = (х -хi+1)/ h>0. Замечание. Полином Нью­то­­на является иной записью полинома Лагранжа для функ­ции, за­дан­ной в рав­но­отстоящих узлах ин­тер­по­ляции. При n = 2 полином Ньютона дает формулу квадратичного интерполирования , а при n = 3 формулу кубического интерполирования . 11.2 Погрешности полиномов Ньютона Оценка погрешности первой и второй ин­тер­по­ля­­ци­он­ных формул Ньютона получается из оцен­ки по­греш­ности интерполяционного поли­но­ма Ла­­гран­жа. Получим эти оценки. Вве­дем обоз­на­че­ние h = хi+1 - хi и полагая для пер­вой ин­тер­по­ля­ци­он­ной фор­мулы q = (х - х0)/h, а для вто­рой q = (х - хn)/h, по­лу­чим для первой ин­­тер­поляционной формулы и для второй интерполяционной формулы . Здесь x - некоторая точка из отрезка ин­­тер­по­ли­ро­ва­ния [х0, хn] в случае задачи ин­тер­по­ли­ро­ва­­ния или не при­над­ле­­жа­щая указанному от­рез­ку в случае за­­дачи экс­тра­по­ли­ро­ва­­ния. Выполняя пре­­­­дель­ный пе­­ре­ход для и под­­­ставляя по­­след­нее вы­ра­же­ние в формулу для по­­греш­ности, окон­­­ча­тель­но бу­дем иметь для пер­вой ин­тер­по­ля­ци­он­­­ной формулы и для второй интерполяционной формулы . Из последних формул получается оценка по­­­­грешности интерполирования, смотри предыдущую лекцию. 11.4 Контрольные вопросы 1. Получить оценку погрешности полиномов Ньютона. 2. Выписать первую формулу Ньютона и формулу погрешности. 3. Выписать вторую формулу Ньютона и формулу погрешности. Лекция 12 Многочлены Чебышева и их свойства План: 12.1 Многочлены Чебышева и их свойства 12.2 Построение полиномов Чебышева 12.3 Заключение 12.4 Контрольные вопросы 12.1 Многочлены Чебышева и их свойства На этой лекции поговорим о полиномах Чебышева. Рассмотрим свойства. Начнем с того, что одним из наиболее удобных приемов вычисления функций состоит в том, чтобы разложить функцию в ряд того или иного типа и заменить ее приближенно с нужной точностью частичной суммой. Есть и другие способы, но самым употребляемым способом приближенных вычислений функций является разложение функции в степенной ряд. Однако с точки зрения вычислительной математики это разложение обладает рядом существенных недостатков. Например, точность аппроксимации (приближение функции частной суммой ряда) велика вблизи центра разложения и мала у границ области сходимости. Поэтому приходится брать много членов ряда, чтобы обеспечить на всем интервале достаточную точность. Конечно, предпочтительнее было бы разложение, дающее одинаковую точность на сем интервале. Если учитывать эти соображения, то наилучшим является разложение по полиномам Чебышева. Русский математик Пафнутий Львович Чебышева поставил и решил задачу об отыскании полиномов «наименее уклоняющихся от нуля ». Более точная математическая постановка этой задачи состоит в следующем. Постановка задачи. Среди всех полиномов заданной степени n со старшим коэффициентом, равным 1, найти тот полином , для которого достигает наименьшего значения. Оказывается, для каждого n существует единственный такой полином. Он и называется полиномом Чебышева степени n. Обозначают его, как правило, Заметим, что ограничение области изменения переменной x отрезком [-1, +1] не уменьшает общности. Любой отрезок [a, b] путем введения новой переменной может быть переведен в отрезок [-1, +1] изменения нового переменного . При этом . у 1 -1 0 1 х рисунок 1 Очевидно, что (рисунок 1) 1) 2) 3) Перейдем к построению полиномов Чебышева в общем виде. 12.2 Построение полиномов Чебышева Полиномы Чебышева можно представить формулами Функция, действительно, является полиномом степени со старшим коэффициентом, равным 1. Покажем это. Для этого рассмотрим выражение Выразим действительную часть через . По формуле Муавра Раскрывая правую часть равенства, получаем Отделяя действительные части и заменяя через , получаем Полагая , получим Отсюда видно, что при нечетных n полином содержит только нечетные степени x, а при четных n - только четные степени x. Экстремальное свойство полиномов Чебышева. Докажем, что полиномы решают задачу, поставленную Чебышевым. Теорема. Среди всех полиномов степени n со старшими коэффициентами, равными 1, рассматриваемых на отрезке [-1, +1] , полином с наименьшей верхней гранью определяется формулой , т.е Доказательство. От противного, допустим, что существует полином степени n со старшим коэффициентом, равным 1, для которого – точная верхняя грань на отрезке [-1, +1] меньше, чем : Рассмотрим где - полином степени (n - 1) относительно n , так как коэффициенты при старших степенях обоих многочленов и равны 1. Рассмотрим на отрезке [-1, +1] точки , в которых достигает наибольшего значения Решая уравнения получим Придавая k значения 0, 1, 2,…, получим (n + 1) различных значений ; таким образом, достигает наибольших по модулю значений в (n + 1) точках По предположению отклонение полинома от нуля на отрезке больше, чем отклонение полинома . Тогда из этого предположения следовало бы, что Таким образом, полином принимал бы попеременно положительные значения в (n + 1) последовательных точках . Следовательно, имел бы, по крайней мере, n нулей. Но – полином степени (n - 1), а поэтому не может иметь n корней, если, конечно, он тождественно не равен нулю. Таким образом, . Кроме того, наибольшее и наименьшее значения равны . Доказанное свойство полиномов Чебышева используется при экономизации степенных рядов. Используя полиномы Чебышева, можно уменьшить число членов в заданных степенных рядах (точнее, в их частичных суммах) без уменьшения точности последующего числового расчета. Замечание. При вычислении частичных сумм ряда на ЦВМ, используя полиномы Чебышева, приходится хранить в памяти меньшее число коэффициентов. Уменьшается также число арифметических операций, нужных для вычисления значения частичной суммы в данной точке. Такая экономизация чрезвычайно желательна в случае подпрограмм, которые должны занимать наименьшее число ячеек памяти и должны быть по возможности более быстрыми. Как же достигается экономизация степенного ряда? Пусть . Вычисляется , с точностью . Полагаем где Оценив остаточный член R, получаем, что для достижения заданной точности расчетов необходимо взять (n + 1) членов ряда При этом . Таким образом, вместо можно брать значения с заданной точностью (точнее < ). Рассмотрим , имеем то есть . Обозначим , тогда . Заменим x в частичной сумме многочленом , т.е. рассмотрим При этом Если при этом , то вместо можно брать без ущерба для заданной погрешности , уменьшив, таким образом, число частичной суммы. Рассмотрим пример. Пример. На отрезке вычисляется функция с точностью . Разлагая sinx в степенной ряд и удерживая 5 членов ряда, будем иметь при . Удерживая лишь 4 члена ряда Тейлора, получили бы , т.е. для достижения нужной точности следует брать не меньше 5 членов ряда Тейлора: . Заменим на , именно: то есть . В заменим многочленом . Тогда . При этом Таким образом, с помощью полинома Чебышева удалось, не понижая точности, уменьшить на один число членов частичной суммы, приближающей sinx. Чтобы вычислять функции с помощью разложения их в ряды по полиномам Чебышева, надо уметь находить коэффициенты этих рядов и их частные суммы. Для вычисления частичных сумм нам будут полезны некоторые рекуррентные соотношения. Рекуррентные соотношения для полиномов Чебышева. Для вывода первого рекуррентного соотношения рассмотрим произведение . Запишем это произведение через сумму и перейдем к полиномам Чебышева с помощью подстановки : . Замечание: Так как полиномы отличаются от лишь множителем , для удобства вычислений впредь будем изучать свойства полиномов Тогда . При m = 1 формулы упрощаются . Подставляя вместо его выражение через x, получим . Зная и , можно получить полиномы Чебышева всех степеней, а именно: Часто бывают нужны разложения степеней x по полиномам Чебышева . Общая формула имеет вид , где – целая часть . Так для первых степеней имеем: Для получения 2-го рекуррентного соотношения между полиномами Чебышева и их производными, продифференцируем основное соотношение: Полагая n = k + 1 и n = k - 1 и вычитая из первого результата, деленного на (k + 1), второй, деленный на (k - 1), получим то есть . Метод Кленшоу для вычисления частичных сумм ряда по полиномам Чебышева Вычисляется с помощью частичной суммы ряда по полиномам Чебышева. Пусть . Для вычисления полагаем . Рассмотрим вспомогательную последовательность , где . По заданным значениям и найти все остальные члены этой последовательности. Задаем . Тогда Тогда . Перегруппируем слагаемые, объединив члены, содержащие множитель , В силу первого рекуррентного соотношения, все коэффициенты при равны нулю. Итак, . Пример. Вычислить с помощью ряда по полиномам Чебышева. При вычислении ограничимся тремя членами ряда Тейлора для функции : . Подставляя получим полагаем . Тогда Описанный метод вычисления частичной суммы ряда по полиномам Чебышева принадлежит Кленшоу. Метод может быть использован не только для вычисления функций, разложенных по полиномам Чебышева, но и при вычислении функций, выраженных через любые другие функции, определяемые рекуррентными соотношениями. Разложение функции в ряд по полиномам Чебышева Для непосредственного получения разложения заданной функции в ряд по полиномам Чебышева может быть использовано свойство ортогональности этих полиномов. Докажем, что полиномы Чебышева ортогональны на отрезке [-1, +1]. Для этого рассмотрим Полагая x = cost, получаем Таким образом, полиномы Чебышева ортогональны на отрезке [-1, +1] с весовой функцией . Коэффициенты разложения функции по любой ортогональной системе функций находятся легко. Пусть , где . Умножив левую и правую части выражения на , проинтегрировав обе части на отрезке [-1, +1], получаем . Все слагаемые в правой части, кроме обращаются в нуль, в силу ортогональности полиномов Чебышева, а потому . Для получения остальных коэффициентов an (n = 1, 2,…) достаточно проделать ту же операцию, домножив левую и правую части равенства на . В результате получим: . Интегралы, входящие в выражения для коэффициентов, могут быть сведены к интегралам вида , а подстановка x = cost преобразует формулы для вычисления коэффициентов к виду: В заключение рассмотрим два примера разложения функций в ряды по полиномам Чебышева. Пример. Разложить в ряд по полиномам Чебышева функцию y = arccos x и вычислить . . Замена переменного дает , , а поэтому для разложения функции y в ряд по полиномам Чебышева достаточно разложить в ряд Фурье. Ограничиваясь 8 членами ряда, получаем Вычисляя по методу Кленшоу, получим Пример. Разложить в ряд по полиномам Чебышева функцию y = |x| на отрезке [-1,+1). В силу четности функции Полагая , получаем ЗАКЛЮЧЕНИЕ Изложенные в курсе лекций методы составляют основу для реше-ния сложных математических задач, являющихся моделью различных отраслей науки и техники. В курсе изложены вычислительные основы этих методов, а для более глубокого освоения методов и их алгоритмического анализа, конечно, данного курса недостаточно. Необходимо изучение другой литературы, чтение статей, знакомство с алгоритмами и возможное участие в разработке и реализации систем, предназначенных для решения задач и их приложений. Приведенный список литературы вполне достаточен для доброт-ного освоения вычислительных методов и возможного углубленного их изучения. Кроме этого, в этих книгах содержится большой выбор интересных задач, которые могут использоваться при выполнении курсовых и дипломных работ, а также в научно – исследовательской работе студентов. полный объем курса находится на кафедре системотехники. БИБЛИОГРАФИЧЕСКИЙ СПИСОК 1. Бахвалов, Н.С. Численные методы/ Н.С. Бахвалов, Н.П. Жидков, Г.М. Кобельков. – М.: Лаборатория Базовых Знаний, 2002. – 632 с. 2. Вержбицкий, В.М. Основы численных методов/ В. М. Вержбицкий. – М.: Высшая школа, 2002. – 840 с. 3. Ващенко, Г.В. Вычислительная математика. Вычислительный лабораторный практикум по численному решению обыкновенных дифференциальных уравнений и систем/ Г.В. Ващенко. - Красноярск: СибГТУ, 2010. - 70 с. 4. Ващенко, Г.В. Вычислительная математика. Основы алгебраической и тригонометрической интерполяции/ Г.В. Ващенко. - Красноярск: СибГТУ, 2008. - 66 с. 5. Ващенко, Г.В. Вычислительная математика. Основы конечных методов решения систем линейных алгебраических уравнений/ Г.В. Ващенко. – Красноярск, СибГТУ, 2005. – 80 с. 6. Безруких, Н.С. Вычислительная математика. Итерационные методы решения систем линейных алгебраических уравнений/ Н.С. Безруких, Г.В. Ващенко. – Красноярск: СибГТУ, 2003. - 76 с. 7. Ващенко, Г.В. Численные методы и программирование. Численные методы решения уравнений. Численное интегрирование/ Г.В. Ващенко, И.С. Матвеева, Н.М. Зубарева.- Красноярск: КГТА, 1995. – 105 с. 8. Калиткин, Н.Н. Численные методы/ Н.Н. Калиткин. – М.: Наука, 1978. – 512 с. 9. Киреев, В.И. Численные методы в примерах и задачах/ В.И. Киреев, А.В. Пантелеев. -М.: Высш. шк., 2006.- 480 с. 10. Лебедев, В.И. Функциональный анализ и вычислительная математика/ В.И. Лебедев. –М.: Физматлит, 2000. – 296 с. 11. Марчук, Г.И. Методы вычислительной математики/ Г.И. Марчук. – М.: Наука, 1980. – 343 с. 12. Каханер, Д. Численные методы и программное обеспечение/ Д. Каханер, К. Моулер, С.Неш - М.: Мир, 1998. - 575 с. 13. Соболь, И.М. Численные методы Монте-Карло/ И.М. Соболь. –М.: Наука, 973.-387 с. 14. Тыртышников, Е.Е. Методы численного анализа/ Е.Е. Тыртышников.-М.: Academia, 2007. -320 с. 15. Фаддеев, Д.К. Вычислительные методы линейной алгебры/ Д.К. Фаддеев, В.Н. Фаддеева. - М.: Физматгиз, 1963.- 734 с. 16. Формалев, В.Ф.Численные методы/ В.Ф. Формалев, Д.Л. Ревизников. -М. :Физматлит, 2006. -400 с. 17. Численный анализ: http://www.szcc.msu.su/num_anal/url_tem/page_4u.htm ПЕРЕЧЕНЬ ОБОЗНАЧЕНИЙ Векторы и матрицы Rn - множество всех n-мерных вещественных векторов x - вектор из Rn (вектор-столбец) xT -вектор, транспонированный по отношению к вектору x из Rn (вектор-строка) - норма вектора x (x, y) = xT y - скалярное произведение векторов x и y из Rn ,- равномерная норма или m -норма вектора x ,- абсолютная норма или l- норма вектора x ,- евклидова или k- норма вектора x r = Ax - b - вектор - невязки A = (aij) - матрица из n строк и n столбцов; aij - элемент матрицы, стоящий в i-й строке и j-ом столбце A = (aij) - матрица верхняя треугольная; элементы aij равны нулю для i > j A = (aij) - матрица нижняя треугольная; элементы aij равны нулю для i < j А = АТ, т.е. если (аij ) = (aji ) - матрица A симметричная AT - матрица, транспонированная к матрице A AАТ= АТА -матрица А нормальная G = diag(aii) -матрица диагональная; диагональные элементы aii ненулевые, остальные равны нулю A-1 - матрица, обратная к матрице А E - единичная матрица; диагональные элементы равны 1, остальные - нулю - норма матрицы А - равномерная норма или m -норма матрицы А - абсолютная норма или l- норма матрицы А - евклидова или k -норма матрицы А AC - произведение матрицы A на матрицу C detA , A - детерминант (квадратной) матрицы А - относительная погрешность (ошибка) вектора cond(A) = - число обусловленности