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

Часть 1. ВВЕДЕНИЕ Лекция 1 ВВЕДЕНИЕ В ТЕОРИЮ ЧИСЛЕННЫХ МЕТОДОВ ЦЕЛЬ ЛЕКЦИИ: определить предмет курса; дать общую характерис­тику таких свойств численных методов как эффективность, точность, устойчивость, экономичность; пояснить эти свойства на простейших вычислительных правилах. Основные требования, предъявляемые к вычислительным алгоритмам Вычислительный алгоритм должен быть устойчивым к ошибкам, обеспечивать требуемую точность, быть эффективным по времени выполнения, быть экономичным по требуемым объемам памяти, не допускать аварийных остановов. Поясним эти требования. Устойчивость. Пусть последовательность величин вычисляется по рекурсивному правилу при заданных и d. Предположим, что при вычислении внесена ошибка i (например, за счет операции округления), т. е. вместо значения использу-ется приближенное значение . В соответствии с рекурсивным правилом Следовательно, и ошибка, допущенная на i-м шаге процесса, на следующем шаге не увеличивается, если операция сложения выполнена без новых округлений. Это означает, что алгоритм устойчив. Рассмотрим другое рекурсивное правило вычисления этой последовательности: Опять y0 и q считаем заданными. Пусть, как и в предыдущем примере, . Тогда , или . В этом случае , и при q>1 ошибка будет возрастать. Такой алгоритм является неустойчивым. Точность. Обозначим через y точное решение задачи, yk – ее приближенное значение, полученное на k-м шаге численного метода. Тогда составляет погрешность решения, а является абсолютной погрешностью. Вычислительный алгоритм должен давать решение с заданной точностью и, следовательно, критерием завершения процесса уточнения решения является выполнение неравенства . На практике в силу трудностей вычисления абсолютной погрешности вместо  используют ее оценку сверху – предельную абсолютную погрешность Δy. В качестве Δy выбирают как можно меньшее значение, удовлетворяющее неравенству , а критерием завершения процесса в этом случае является неравенство . Часто используют относительную погрешность т. е. абсолютную погрешность на единицу измерения, и предельную относительную погрешность . Критерий завершения процесса в этом случае имеет вид , где – заданная допустимая относительная ошибка. Рассмотрим в качестве примера, как может быть построена оценка предельной абсолютной погрешности вычисления значений функции , если известны абсолютные погрешности ее аргументов. Ошибка имеет вид . Тогда , и в качестве предельной относительной погрешности можно использовать величину . Пусть решением задачи является вектор , который будем рассматривать как элемент векторного пространства . Приближенное решение и его погрешность также являются элементами этого пространства. Рассмотрим, как оцениваются ошибки решения в этом случае. Для этого используем количественную характеристику вектора в виде нормы. Говорят, что в задана норма, если каждому вектору из сопоставляется вещественное число, называемое нормой вектора , для которого справедливы следующие свойства: 1. , причем тогда и только тогда, когда , 2. для любого вектора и любого числа , 3. для любых векторов и . В вычислительных методах наиболее употребительны следующие нормы: . Абсолютную и относительную погрешность вектора в любой из перечисленных выше норм можно определить следующим образом: . Имеет смысл говорить об абсолютной и относительной погрешности (mn)- матрицы решений A. В этом случае используется такая характеристика, как норма матрицы, согласованная с нормой вектора. Для нормы матрицы A, обозначаемой , справедливы следующие свойства: 1. , причем тогда и только тогда, когда, 2. для любой матрицы A и любого числа  , 3. для любых (mn)-матриц и. Каждой из векторных норм соответствует своя согласованная норма матриц. В частности, нормам соответствуют нормы , вычисляемые по правилам где j(ATA) – собственные числа матрицы ATA. Абсолютная и относительная погрешности матрицы решений вычисляются через соответствующие нормы: , причем – искомая матрица решений, – ее некоторое приближение. Эффективность. Эффективность вычислительного алгоритма оценивают по времени реализации либо по общему числу операций, требуемых для получения результата. Обозначим через nсл число операций сложения, через nу – число операций умножения, через nдел – число операций деления, через tу, tдел, tсл – время выполнения операций умножения, деления и сложения на ЭВМ. Объективный критерий сравнения вычислительных алгоритмов – требуемое время вычислений: T=nуtу+nслtсл+nделtдел. Если tу≈tдел≈tсл, то T≈(nу+nдел+nсл)tоп, где tоп – время выполнения на ЭВМ арифметической операции, и критерием оценки эффективности алгоритма может быть nΣ=nу+nдел+nсл . Если tу≈tдел>>tсл, то критерий эффективности – число так называемых длинных операций nдл=nдел+nу. Пусть Q – критерий эффективности метода. Значение Q зависит от требуемой точности . Из двух сравниваемых вычислительных алгоритмов решения задачи из них будет эффективнее тот, который при заданном значении ε характеризуется меньшим значением Q. В качестве примера построим два алгоритма вычисления полинома n-й степени . Непосредственная реализация: Расчет полинома по рекурсивному правилу , вытекающему из приведенной реализации, соответствует алгоритму 1. Алгоритм 1. 1. Начальная установка 2. Вычислить в цикле (k=1,2,…,n): Алгоритм 1 характеризуется таким числом арифметических операций: nу=2n; nсл=n. Преобразуем исходную запись полинома: Теперь расчет полинома можно выполнить следующим образом: Записанное правило реализуется алгоритмом 2. Алгоритм 2. 1. Начальная установка: 2. Вычислить в цикле (k=n,n-1,…,1): Значение . Затраты на выполнение алгоритма: nу=n , nсл=n. Очевидно, что алгоритм 2 эффективнее по времени выполнения алгоритма 1. Экономичность. Экономичность алгоритма оценивается по требуемым для его реализации объемам памяти ЭВМ. Аварийные остановы. Любое число ЭВМ принадлежит не всей числовой оси, а некоторому интервалу . Если какое-либо число x в процессе работы алгоритма выходит за указанный интервал, то происходит так называемый аварийный останов. Необходимо построить вычислительный алгоритм таким образом, чтобы вычисленные в процессе его работы числа . Часть 2. системЫ линейных алгебраичЕских уравнений Лекция 2 ПРЯМЫЕ МЕТОДЫ ЦЕЛЬ ЛЕКЦИИ: Определить два класса численных методов (прямые и итерационные); показать, как строятся прямые методы Гаусса, LU-факторизации, Холесского; выполнить оценку их эффективности. Постановка задачи. Основная задача вычислительной алгебры – решение систем линейных алгебраических уравнений (СЛАУ) В дальнейшем будем использовать запись этой системы в компактной форме: ( запись означает, что индекс i изменяется от 1 до n с шагом 1), или в векторном виде , где Предполагается, что матрица неособенная, т. е. , и решение единственно. Прямые и итерационные методы. Численные методы решения СЛАУ делятся на две большие группы: прямые и итерационные. Прямые методы при отсутствии ошибок округления за конечное число арифметических операций позволяют получить точное решение . В итерационных методах задается начальное приближение и строится последовательность , где k – номер итерации. В действительности итерационный процесс прекращается, как только становится достаточно близким к . Имеется промежуточный класс методов, в которых решение ищется итерационно, однако для них заранее известно, какое число итераций необходимо выполнить, чтобы в отсутствии ошибок округления получить точное решение. На практике при вычислении приближенного решения число итераций в наиболее эффективных методах оказывается значительно меньшим, чем этого требует теория точного решения. Какой класс методов лучше? Однозначно на этот вопрос ответить нельзя. Итерационные методы привлекательнее с точки зрения объема вычислений и требуемой памяти, когда решаются системы с матрицами высокой размерности. При небольших порядках системы используют прямые методы либо прямые методы в сочетании с итерационными методами. Метод Гаусса. В методе Гаусса линейная система решается в два этапа. На первом этапе система преобразуется к виду (см. рис. 2.1) , Рис. 2.1. Структура системы и портрет ее ненулевых элементов до (а) и после (б) прямого хода Гаусса где – верхняя треугольная матрица с единичной диагональю (это так называемый прямой ход Гаусса). На втором этапе (обратный ход Гаусса) решается система . Рассмотрим эти этапы подробнее. Прямой ход. Прямой ход Гаусса состоит из n шагов. Первый шаг. Полагаем, что и разделим на него первое уравнение. Перепишем систему с учетом этого преобразования: Умножим первое уравнение на и вычтем его из i-го уравнения преобразованной системы: Обозначим . Получим Второй шаг. На втором шаге из системы исключается аналогичным образом: K-й шаг. Запишем общий вид преобразованной системы после k-го шага прямого хода Гаусса: Здесь Проиллюстрируем, как меняется матрица системы в процессе прямого хода Гаусса на примере системы четвертого порядка (рис. 2.2; ненулевые элементы матрицы обозначены крестиками). Рис. 2.2. Преобразование матрицы системы 4-го порядка на прямом ходе Гаусса Оценим количество длинных операций (умножений и делений) на первом шаге прямого хода Гаусса. Преобразование первого уравнения требует n таких операций. Преобразование остальных n-1 уравнений – n(n-1) операций умножения и деления. Таким образом, первый шаг выполняется за длинных операций. Рассуждая по аналогии, нетрудно найти затраты на остальных n-1 шагах. Суммарные затраты прямого хода Гаусса определяются в итоге рядом . Последняя оценка имеет место для n>>1. Обратный ход. Запишем систему, решаемую на обратном ходе, в координатном виде Ее решение: Запись означает, что индекс k изменяется от значения n-1 до 1 с шагом 1. Требуемое число длинных операций на обратном ходе Приближенная оценка справедлива для n>>1. Общие затраты метода Гаусса: Таким образом, при больших n основные затраты в методе Гаусса приходятся на прямой ход. Метод LU-факторизации. В методе LU-факторизации (эту схему называют компактной схемой Гаусса) при решении системы выполняется следующая последовательность действий. Матрица представляется в виде произведения , Рис. 2.3. Структура матриц L и U в разложениях Дулиттла (а) и Краута (б) где L  нижняя треугольная матрица, U  верхняя треугольная матрица. Такое разложение единственно при условии предварительного выбора диагональных элементов одной из матриц. В этом случае число элементов в матрице A совпадает с суммарным числом неизвестных элементов матриц L и U. Если диагональ L принимается единичной, то такое разложение называют разложением Дулиттла (рис. 2.3,а), если единична диагональ U – разложением Краута (рис. 2.3,б). В дальнейшем при построении метода LU-факторизации будем привлекать разложение Краута. Система заменяется системой , легко решаемой за два шага: Шаг 1. . Принимая во внимание треугольный вид матрицы L, нетрудно получить, что в алгоритме Краута Шаг 2. . Решение этой системы в алгоритме Краута: . Суммарные затраты реализации обоих шагов при n>>1 составляют длинных операций. Получим соотношения для расчета элементов матриц L и U в алгоритме Краута. Для этого перемножим матрицы L и U и приравняем результат к A. По правилу перемножения матриц Учтем, что Рассмотрим элемент (рис. 2.4), расположенный на центральной диагонали либо в нижней треугольной части матрицы A. Для такого элемента i ≥ j. Из рисунка следует, что Рис. 2.4. Иллюстрация вычисления элемента матрицы, расположенного ниже главной диагонали , так как i ≥ j и . Отсюда Рассмотрим элемент (рис. 2.5), находящийся выше главной Рис. 2.5. Иллюстрация вычисления элемента матрицы, расположенного выше главной диагонали диагонали матрицы A (для него j>i). В этом случае Следовательно, Получили в итоге соотношения, которые позволяют вычислять элементы матриц L и U. Последовательность вычислений: сначала вычисляется столбец матрицы L, далее строка матрицы U, затем опять столбец матрицы L, далее строка матрицы U и т. д. (см. рис. 2.6, который иллюстрирует последовательность вычислений и схему хранения матриц L и U). Вычисление столбца матрицы L и строки матрицы U назовем шагом LU-разложения. Приведем в качестве примера схему хранения элементов матриц A,L,U после второго шага LU-разложения (рис. 2.7). Число длинных арифметических операций на этапе LU-разложения при n>>1 составляет величину , на шаге решения ли- нейных систем с треугольными матрицами – . Суммарное число длинных операций приближенно равно (как и в методе Гаусса), Рис. 2.6. Исходная матрица A (а), схема хранения L и U матриц (б), последовательность вычисления элементов в принятой схеме хранения (в) Рис. 2.7 Схема хранения элементов 44 матриц A, L, U после второго шага LU-факторизации т. е. основные затраты приходятся на LU-факторизацию матрицы A. Эта особенность делает особо привлекательным метод LU-факторизации при решении СЛАУ с одной и той же матрицей A, но разными правыми частями: В этом случае факторизация матрицы выполняется однократно, требуя длинных операций, а решение каждой системы с соответствующей правой частью реализуется за таких операций. Метод Холесского. Исключительно эффективную реализацию метода LU-факторизации можно получить, если ограничиться классом линейных систем с симметрической положительно определенной матрицей A, т. е. Такую реализацию называют методом Холесского, либо методом квадратного корня. Будем полагать, что решаемая система имеет симметрическую положительно определенную матрицу A. В этом случае матрица A представляется в виде Здесь – нижняя треугольная матрица. Такое разложение существует и единственно для положительно определенных симметрических матриц. Система преобразуется к виду . Вектор ищется путем последовательного решения двух систем с треугольными матрицами: ; . Для получения расчетных соотношений элементов матрицы рассмотрим произвольный элемент матрицы A: Суммирование здесь выполняется только до j, т. к. j≤i. Выделим член при значении k=j: . Теперь Эти соотношения позволяют вычислить по столбцам элементы матрицы . Эффективность такого метода достигается на этапе разложения матрицы, т. к. необходимо вычислить в этом случае только матрицу . Арифметические затраты в методе Холесского составляют длинных операций и n операций извлечения квадратного корня. Существует другой вариант разложения симметрической положительно определенной матрицы, в котором удается избежать операций извлечения квадратного корня. В этом варианте вводится новая матрица по правилу , причем – матрица, вычисленная ранее по схеме Холесского, – диагональная матрица с элементами матрицы . Матрица существует и является нижней треугольной с единичной диагональю. В этом случае , где . Расчетные соотношения для элементов матриц и можно получить, как и прежде, привлекая правило перемножения матриц , из которого следует, что т. к. матрица имеет единичную диагональ. Такой алгоритм потребует вдвое большего числа перемножений, чем схема Холесского. Однако, если ввести замену переменных , то расчетные соотношения примут вид Здесь сначала вычисляют вспомогательные величины , а затем их используют для определения искомых величин и . Количество умножений при такой организации алгоритма составляет приблизительно и не содержит операции извлечения квадратного корня. Лекция 3 УСТОЙЧИВОСТЬ И ТОЧНОСТЬ ПРЯМЫХ МЕТОДОВ ЦЕЛЬ ЛЕКЦИИ: Выполнить оценку устойчивости и точности прямых методов; показать, как перестановкой строк и столбцов обеспечивается устойчивость и точность прямых методов, каким образом осуществляется выбор матриц перестановки строк и столбцов в случае систем с разреженной матрицей. Устойчивость алгебраических методов. Прямые методы в отсутствии ошибок округления позволяют получить точное решение системы . Современные вычислительные машины оперируют с конечными десятичными дробями, представленными в форме с плавающей точкой. В этом случае уже на этапе запоминания элементов матрицы A и вектора в памяти ЭВМ вносится ошибка округления и реально решается возмущенная система . Здесь и – возмущения матрицы системы и вектора правой части. Для элементов матрицы и компонент вектора справедливы оценки , где– элемент матрицы A, – компонент вектора , ε – число, характеризующее относительную погрешность машинной ариф-метики. Например, для двоичных ЭВМ, использующих арифметику с плавающей точкой и t – разрядную мантиссу, . Перейдем к более общей числовой оценке возмущений – норме. Из записанных выше неравенств следует, что , где знаком обозначена какая-либо норма вектора и согласованная с ней норма матрицы. Поясним теперь суть одного из наиболее разработанных подходов к анализу устойчивости алгебраических методов. Пусть – приближенное решение СЛАУ, полученное применением некоторого прямого метода. Очевидно, что вследствие ошибок округления при реализации на ЭВМ прямого метода это решение не будет точно удовлетворять системе . Пусть, однако, удается показать, что является точным решением системы . Если для матрицы F и вектора , называемых эквивалентными возмущениями метода, можно получить оценки вида , где f(n), h(n)  некоторые степенные функции типа с небольшим показателем k, то метод считается устойчивым по Уилкинсону. Такой вид функций f(n), h(n) объясняется тем, что в процессе реализации прямых методов на ЭВМ неизбежно некоторое накопление ошибок округления, пропорциональное числу арифметических операций, за которое прямой метод приводит к решению. Точность алгебраических методов. Рассмотрим влияние возмущений матрицы системы и вектора правой части на точность решения. Относительная погрешность решения , где – точное решение, в любой согласованной норме подчиняется неравенству . Обозначим . Величины , являются оценками относительных возмущений матрицы A и вектора b. Число χ есть не что иное, как число обусловленности матрицы A. Естественно допустить, что . Такое неравенство означает, что норма матрицы возмущений существенно меньше нормы исходной матрицы, т. е. предполагается, что прямой метод решения линейной системы устойчив. В этом случае . Относительная ошибка решения системы прямым методом может достигать величины, определяемой произведением числа обусловленности и суммы относительных возмущений матрицы системы и вектора правой части. Если число обусловленности велико, то даже при малых эквивалентных возмущениях матрицы и вектора правой части погрешность решения линейной системы оказывается значительной. Задачи, для которых χ велико, называют плохо обусловленными. Число обусловленности зависит от выбора матричной нормы. Например, для симметричной положительно определенной матрицы A в качестве можно взять ее максимальное собственное число. Тогда число обусловленности , где , – максимальное и минимальное собственные числа матрицы. Устойчивость метода Гаусса. Опуская обременительные преобразования в методе обратного анализа ошибок округления, отметим, что возмущенная система метода Гаусса имеет вид . Запишем оценку нормы матрицы возмущения: . Вид этой оценки удовлетворял бы критерию устойчивости Уилкинсона, если бы множитель g(A) имел небольшое значение. Поясним смысл множителя g(A). Пусть обозначает матрицу, полученную из A после k шагов исключения. Обозначим . Тогда . Следовательно, g(A) показывает, во сколько раз могут возрасти элементы матрицы A в ходе исключения переменных по сравнению с их исходным уровнем. По этой причине g(A) называют коэффициентом роста матрицы A . Элементы активной части матрицы Ak в методе Гаусса вычисляются по формуле . Для ограничения роста элементов матрицы в процессе гауссова исключения желательно, чтобы поправочные члены в этой формуле были не слишком большими. Это достигается процедурой выбора элемента , который называют главным. Выбор главного элемента по столбцу. В этом случае ограничение роста элементов матрицы Ak на k–м шаге гауссова исключения достигается перестановкой строк таким образом, чтобы гарантировать неравенство . С этой целью при исключении переменной в качестве главного элемента выбирается элемент матрицы Ak-1 по правилу , т. е. наибольший по модулю элемент в k–м столбце матрицы Ak-1 (рис. 3.1). Строки r и k переставляются и только после этого выполняется k–й шаг исключения прямого хода Гаусса. При столбцовой стратегии выбора главных элементов справедлива такая оценка для значения параметра ak, определяющего коэффициент роста: . Она допускает, что и, следовательно, коэффициент роста . По этой причине метод Гаусса с выбором главного элемента по столбцам является условно устойчивым. Несмотря на это, он широко используется на практике, так как g(A) редко достигает своего верхнего предела. Выбор главного элемента по всей матрице. В этой стратегии в качестве главного элемента при исключении неизвестной xk выбирается элемент по правилу , т. е. наибольший по модулю элемент в квадратной подматрице матрицы Ak-1 (рис. 3.2). Строки k и r, а также столбцы k и l переставляются и далее выполняется k–й шаг исключения. Такая стратегия гарантирует выполнение неравенства и, следовательно, ограничивает рост элементов в процессе исключения Гаусса. Оценка коэффициента роста элементов матрицы A в этом случае имеет более благоприятный вид: . Точность метода Гаусса. Привлекая оценку нормы матрицы возмущения, можно записать, что . Анализ неравенства позволяет определить пути повышения точности метода Гаусса: выбор главных элементов, работа с числами удвоенной длины, переобусловливание системы линейных алгебраических уравнений. Устойчивость и точность метода LU-факторизации. Метод LU-факторизации, если не делать никаких специальных усилий, характеризуется такими же оценками нормы матрицы возмущения и относительной ошибки, как и метод Гаусса. В то же время этот метод наряду с упоминавшимися при рассмотрении метода Гаусса направлениями повышения точности решения имеет еще одну возможность: использование так называемой операции скалярного накопления. Эта операция во многих ЭВМ реализована аппаратно либо моделируется специальными подпрограммами. Операция скалярного накопления привлекается при вычислении скалярного произведения . Пусть основным режимом работы ЭВМ является арифметика с плавающей точкой и t-разрядной мантиссой машинного слова. Рассмотрим вычисление скалярного произведения. Произведение чисел βk и γk, выполненное арифметическим процессором, имеет 2t–разрядную мантиссу. Вместо того, чтобы округлить это произведение до t разрядов и прибавить результат к хранимой t–разрядной промежуточной сумме Sk-1, режим накопления предусматривает удержание 2t разрядов как для Sk-1, так и для βkk. Очередное суммирование выполняется с 2t-разрядными числами. Оно также сопровождается ошибкой округления, но эта ошибка затрагивает лишь последний разряд 2t–разрядного результата. Возвращение к t–разрядному представлению происходит лишь после вычисления скалярного произведения. Анализ расчетных соотношений метода LU-факторизации свидетельствует о том, что они реализуются с использованием операции скалярного накопления. Оценка нормы матрицы возмущений метода LU-факторизации при использовании операции скалярного накопления имеет вид . Устойчивость вычислительной схемы метода определяется величиной коэффициента роста элементов матрицы A. Чтобы уменьшить величину g(A), в методе LU-факторизации используется столбцовая процедура выбора главных элементов. Погрешность решения метода . Устойчивость и точность метода Холесского. Этот метод применяется для решения линейных систем с симметрической положительно определенной матрицей A. Можно показать, что в схеме Холесского коэффициент роста элементов матрицы A равен 1. Действительно, из соотношения следует, что . Если учесть, что в симметрической положительно определенной матрице A наибольший по модулю элемент положителен и находится на главной диагонали, нетрудно прийти к сформулированному заключению: g(A)=1. Таким образом, чтобы ошибки округления при реализации схемы Холесского были малыми, не требуется никаких перестановок строк либо столбцов. Наименьшая погрешность численного решения достигается при привлечении операции скалярного накопления. В этом случае и . Прямые методы для систем с разреженными матрицами. Во многих приложениях приходится решать систему с разреженной матрицей A. Назовем –матрицу A разреженной, если число ее ненулевых элементов . При решении таких систем естественно хранить в памяти ЭВМ только ненулевые элементы. Можно надеяться, что в случае разреженных матриц как вычислительная работа, так и требуемый объем памяти значительно сократятся. Проблема состоит в том, что на этапе исключения неизвестных в методе Гаусса либо на этапе построения L и U матриц в компактной схеме Гаусса могут возникать Рис. 3.3. LU-факторизация диагональной матрицы со строчным и столбцовым окаймлением новые ненулевые элементы. Обратимся к следующему простому при-меру (рис. 3.3). Хотя исходная матрица A здесь имеет ничтожное число ненулевых членов, матрицы L и U имеют такую же структуру Рис. 3.4. LU-факторизация преобразованной диагональной матрицы со строчным и столбцовым окаймлением ненулевых элементов, как и для полностью заполненной матрицы. В то же время перестановкой строк и столбцов систему можно преоб-разовать к виду, который сохраняет свойство разреженности исходной матрицы в матрицах разложения (рис. 3.4). В общем случае, конечно, не всегда удается полностью сохранить разреженность. Поэтому перестановки строк и столбцов выбирают такими, чтобы заполняемость матрицы в процессе исключения переменных была минимальной. На языке матричного описания эта задача выглядит следующим образом. Пусть P и Q – -матрицы перестановок строк и столбцов, причем . Матрицы P и Q легко получить из единичной матрицы, в которой переставлены соответствующие строки (столбцы). Например, матрица перестановки первой и третьей строк -матрицы приведена на рис. 3.5. Преобразуем систему линейных алге-браических уравнений следующим образом: или где – переупорядоченная форма A, , . К системе будем применять прямой метод решения (его эффективность зависит от выбора матриц P и Q) и далее найдем , так как . Задача осложняется тем, что перестановка строк и столбцов определяющим образом влияет не только на разреженность матриц разложения, но и на численную устойчивость прямого метода. Симметричные матрицы. Если матрица A линейной системы является симметрической положительно определенной, то процесс исключения переменных по схеме Холесского является численно устойчивым при любом выборе главных элементов из числа диагональных членов. Таким образом, если положить , т. е. при перестановке k–й и j–й строк осуществлять перестановку k–го и j–го столбцов, то выбор матрицы P можно производить только исходя из требования минимизации заполнения матрицы . Задача численного решения линейной системы сводится в этом случае к выполнению трех последовательных шагов: Шаг 1. Формирование матрицы и вектора . Шаг 2. Решение упорядоченной системы . Шаг 3. Вычисление вектора . Рассмотрим более подробно, как реализуется Шаг 1. Вообще говоря, для -матрицы A существует различных упоря-дочиваний, из которых одно или более оказывается оптимальным в смысле заполнения. Задача отыскания оптимальных упорядочиваний является чрезвычайно трудной и практически нерешенной на сегодняшний день. Существующие процедуры реализации шага 1 являются эвристическими. Они обеспечивают понижение заполнения матрицы , однако не гарантируют достижение точного минимума. Опишем в общих чертах одну из таких процедур – алгоритм минимальной степени. Основная идея алгоритма минимальной степени заключается в локальной минимизации заполнения за счет выбора главного элемента в той строке и столбце, которые обеспечивают внесение наименьшего числа ненулевых элементов в множитель . Рассмотрим, например, k–й шаг исключения прямого хода Гаусса. В поддиагональных позициях столбцов с 1–го по (k-1)–й расположены нулевые элементы. Для исключения переменной xk на этом шаге k–я строка, нормированная соответствующим образом, вычитается из всех тех строк, которые имеют ненулевые элементы в k–м столбце. Следовательно, чтобы минимизировать число новых ненулевых элементов, достаточно просмотреть активную подматрицу матрицы Ak-1, выбрать строку с наименьшим числом ненулевых элементов, назовем ее l–й строкой, и переставить как строки, так и столбцы с номерами l, k. После такой перестановки в матрицу Ak-1 вносятся новые ненулевые элементы, т. е. формируется портрет матрицы Ak , и процесс имитации исключения переменных продолжается. Выполнив имитацию (n-1)–го шага исключения, построим в итоге из единичной матрицы I искомую матрицу перестановок P. Матрицы общего вида. На k–м шаге алгоритма Гаусса при решении линейных систем с разреженной матрицей A произвольной структуры необходимо выбирать в качестве главного элемента такой элемент , который прежде всего обеспечивает устойчивость вычислительного процесса и в то же время по возможности минимизирует заполняемость матрицы Ak. Это достигается введением специального барьера – числа u≥1 и выделением множества элементов из активной подматрицы матрицы Ak-1, которые удовлетворяют условию , где – элемент множества Sk-1. Для каждого из таких элементов вводится числовая оценка , где r(l,k) – число ненулевых элементов в l–й строке активной части матрицы Ak-1, c(m,k) – число ненулевых элементов в m–м столбце этой же части матрицы Ak-1. В качестве главного элемента на k–м шаге прямого хода Гаусса выбирается такой элемент , для которого . Лекция 4 ИТЕРАЦИОННЫЕ МЕТОДЫ ЦЕЛЬ ЛЕКЦИИ: Дать общую схему линейных итерационных методов первого порядка и построить на её основе методы простой итерации, Гаусса–Зейделя, последовательной релаксации; сравнить эти методы по вычислительным затратам на итерации. В случае больших разреженных СЛАУ предпочтение отдаётся итерационным методам, т. к. они, во-первых, не приводят к появлению на итерациях новых ненулевых элементов, во-вторых, оказываются более эффективными по затратам машинного времени. Схема итерационных методов. Пусть требуется решить систему . Предположим, что Q – неособенная -матрица, рассматриваемая как апроксимация матрицы A, для которой решение системы не представляет особого труда. Это имеет место тогда, когда матрица Q является нижней треугольной, верхней треугольной, трёх-диагональной, произведением конечного числа таких простых матриц и т.п. Матрицу Q называют матрицей расщепления. Перепишем исходную систему в виде . Построим итерационный процесс Разрешим его относительно: , или . Введем обозначения: ; . Тогда итерационную схему можно представить в виде Это линейная схема первого порядка. Поскольку итерационная матрица B на итерациях является постоянной, то такую схему называют стационарной. Достаточное условие сходимости этой итерационной схемы имеет вид . Убедимся в этом. Обозначим вектор погрешности на k-м шаге через , на (k+1)-ом шаге –, точное решение – . Тогда . Учет этих соотношений в итерационной схеме позволяет записать Отсюда . При видно, что , т. е. метод сходится к решению. Построим теперь вычислительные схемы, выбирая конкретный способ формирования матрицы раcщепления. Метод простой итерации (Метод Якоби). Пусть требуется решить систему . Представим , где D – диагональная матрица с диагональными членами матрицы A; – часть матрицы A, лежащая ниже центральной диагонали, – часть матрицы A, лежащая выше центральной диагонали. Тогда , или . Запишем итерационный метод Разрешим его относительно : Эта форма удобна для реализации метода Якоби. Здесь итерационная матрица . Матрица расщепления . Следовательно, при построении метода Якоби в качестве матрицы ращепления используется диагональ матрицы СЛАУ. Запишем метод Якоби в координатной форме для системы общего вида: Нетрудно убедиться, что метод Якоби в координатной форме есть не что иное, как разрешение каждого из уравнений системы относительно одной из компонент вектора. Из первого уравнения системы выражается и его значение принимается за значение . Из второго уравнения определяется и его значение принимается за и т. д. Переменные в правой части этих соотношений при этом полагаются равными их значениям на предыдущей итерации. Рассмотрим пример системы линейных алгебраических уравнений с разреженной (пятидиагональной) матрицей A (рис. 4.1) и оценим эффективность решения такой системы методом Якоби. Будем хранить матрицу A в виде пяти одномерных массивов g, f, a, b, c. Тогда решаемая система в терминах указанных массивов примет вид . В этой записи следует учесть, что отдельные коэффициенты являются нулевыми. В частности при i=1; при ; при ; при Рис. 4.1. Система с пятидиагональной матрицей Реализация метода Якоби для такой системы линейных алгебраических уравнений осуществляется просто: Видно, что для вычисления компоненты вектора решения необходимо выполнить четыре операции умножения и сложения и одну операцию деления. Полное время одной итерации метода Якоби оценивается соотношением Метод Гаусса–Зейделя. Решаемую систему запишем в виде . Итерационная схема Гаусса–Зейделя также следует из этого пред-ставления системы: , или . Последняя форма такого итерационного метода удобна для реа-лизации. Приведем метод Гаусса–Зейделя к стандартному виду: Стандартная форма метода позволяет выписать его итерационную матрицу и провести над ней очевидные преобразования: . Матрица расщепления cодержит всю нижнюю треугольную часть матрицы A. Запишем метод Гаусса–Зейделя в координатной форме для системы общего вида: Координатная форма метода Гаусса–Зейделя отличается от координатной формы метода Якоби лишь тем, что первая сумма в правой части итерационной формулы содержит компоненты вектора решения не на k-й, а на (k+1)-й итерации. Но именно это означает, что матрица расщепления метода Гаусса–Зейделя включает не только диагональные члены матрицы A, но и все члены, расположенные ниже главной диагонали. Оценим, как и ранее, временные затраты метода при решении тестовой задачи  СЛАУ с пятидиагональной матрицей. Реализация метода Гаусса–Зейделя для такой системы: позволяет констатировать, что временные затраты на итерации метода оказываются такими же, как и в методе Якоби. Ключ к пониманию более высокой эффективности метода Гаусса–Зейделя  количество итераций, требующееся для получения решения с заданной точностью. Более высокая сходимость к решению в методе Гаусса–Зейделя достигается за счет выбора матрицы расщепления, лучше аппроксимирующей матрицу A. Метод последовательной релаксации. Воспользуемся методом Гаусса–Зейделя, записанным в форме , или . Удобная для реализации итерационная схема релаксационного метода имеет вид , где ω – некоторый параметр. Преобразуем ее к стандартной форме: . В свою очередь итерационная матрица . В итоге матрица расщепления . При указанная итерационная схема соответствует методу последовательной верхней релаксации, при – методу последовательной нижней релаксации. Координатная форма метода последовательной релаксации: Метод последовательной релаксации для системы с пятидиа-гональной матрицей: Время выполнения итерации в этом случае . Метод касательных. Метод касательных (Ньютона) применяется для уточнения действительных корней нелинейного уравнения. Имеем уравнение F(X) = 0. Если F(X) Î D[A, B], и известно F(А) F(В) < 0 (F(А) и F(В) имеют разные знаки) и F11(X) не меняет знака на интервале [А,В], то уравнение F(X) = 0 имеет решение Х* Î [A, B], F(X*) º 0. Процесс уточнения действительного корня нелинейного уравнения методом касательных заключается в том, что строится касательная к кривой F(X) и определяется точка пересечения касательной с осью абсцисс, координата этой точки используется в качестве уточнения корня. Затем к кривой F(X) в точке последнего уточнения корня строят очередную касательную и определяется точка пересечения касательной с осью Х – очередное уточнение корня уравнения. Процесс уточнения корня (построение касательных) продолжается до тех пор пока два ближайших уточнения будут отличаться на величину не более e, точность вычисления корня, | Х(i) – X(i-1) | < e где i - номер итерации.. Ограничение по применению метода. Метод касательных можно применять для уточнения действительного корня, если на интервале [А, В] функция F(X) удовлетворяет следующим свойствам: интервал [А,В] должен быть достаточно мал, чтобы на его длине график функции F(X) не имел горизонтальных участков, участков с малым наклоном касательных к кривой F(X) и экстремумов, функция F(X) должна быть монотонная, F11(А) F11(В) > 0; за начальное приближение Х0 принимается одна из границ интервала [А,В], где F(X) и F11(X) имеют одинаковые знаки, Х0 = А при условии F(А) F11(А) > 0 (Х0 = B при условии F(B) F11(B) > 0). Алгоритм метода. Исходные данные: начальное приближение корня Х, точность вычисления корня e . Организовать вычисление поправки к уточнению корня, DХ = F(X) / F1(X), Уточнить значение корня. В качестве рекуррентной формулы метода касательных используют формулу вычисления координаты пересечения касательной с осью Х, Х = Х - DХ, Проверить условие продолжения уточнения корня, | DХ | > e . Если заданное условие принимает значение «истина», то необходимо продолжить уточнение корень уравнения с пункта 2. Если заданное условие принимает значение «ложно», то корень Х найден с заданной точностью e, организовать вывод значения корня Х и прекратить уточнение корня. Часть 3. НЕЛИНЕЙНЫЕ АЛГЕБРАИЧЕСКИЕ УРАВНЕНИЯ И СИСТЕМЫ Лекция 5 ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ. МЕТОД ИТЕРАЦИй. ЦЕЛЬ ЛЕКЦИИ: Для уравнения ввести постановки задач отделения и уточнения корней, сформулировать лемму об оценке погрешности приближенного решения, построить метод простой итерации, получить для него оценку точности, доказать его сходимость. Пусть требуется решить уравнение , т. е. найти все корни , удовлетворяющие этому уравнению на отрезке . Задача численного решения уравнения сводится, во-первых, к отделению корней, во-вторых, к последующему уточнению корней. Отделение корней. Отделить корни уравнения значит заключить каждый корень в интервал , для которого выполняются условия: 1. ; 2. – знакопостоянная функция для . Первое неравенство обеспечивает наличие в интервале хотя бы одного корня, второе условие гарантирует единственность корня (см. рис. 4.1). Для отделения корней можно использовать аналитический и табличный способы. Аналитический предполагает исследование функции методами математического анализа, последующее построение графика функции, из которого и определяются интервалы, содержащие единственный корень. Недостаток аналитического способа – неалгоритмизуемость. Табличный метод предполагает составление таблицы значений, причем. Из таблицы на основании условий алгоритмически определяются искомые интервалы . Чтобы не потерять корни, интервал отделения h должен быть достаточно малым. Уточнение корней. Пусть - корень уравнения . В дальнейшем – интервал, содержащий единственный корень. Уточнить корень с заданной точностью значит найти приближенное значение такое, что . Для решения сформулированной задачи необходимо уметь производить оценку абсолютной погрешности метода, т. е. величины . Лемма об оценке погрешности приближенного решения уравнения . Лемма. Пусть уравнение на отрезке имеет корень . Пусть найдено некоторое его приближенное значение . Тогда , где . Доказательство. Вычислим по теореме о конечных приращениях: . Очевидно, что . Отсюда , или . Лемма доказана. Величину называют невязкой. Из леммы следует, что судить о величине погрешности приближенного решения только по величине невязки нельзя. Необходимо учитывать и значение первой производной. Обратимся к рис. 4.2. Из этого рисунка следует, что одна и та же невязка приводит к существенно разным погрешностям приближенного решения, если производная в окрестности решения сильно отличается. Метод итераций. Для построения метода итераций преобразуем уравнение к виду . Это можно сделать в общем случае так: , или , где . Постоянный множитель h выбирается при этом из условия сходимости метода (установим его позже). Пусть известно начальное приближение . Тогда Приведенный способ построения числовой последовательности реализуется в методе итераций: . Рассмотрим, как ведет себя погрешность решения на итерациях метода. Обозначим , где - погрешности приближенного решения на двух соседних итерациях. Подставим представленные таким образом и в итерационное правило: . Разложим функцию в ряд Тейлора в окрестности точки: . Пренебрегая остаточным членом , получим соотношение, связывающее погрешность решения метода на двух соседних итерациях: . Сделаем некоторые выводы на основании этого приближенного равенства. • Если , то можно ожидать, что и последовательность будет сходиться к решению, когда начальное приближение выбрано достаточно близким к . • Если , то скорее всего и метод будет расходиться, так как каждое последующее приближение будет отстоять от решения дальше, чем предыдущее. • При и погрешности и имеют одинаковые знаки. Сходимость будет монотонной. • При и погрешности и имеют разные знаки. Сходимость является немонотонной. Проиллюстрируем характер сходимости метода итераций графически. Образуем функции. В решении задачи значения этих функций совпадают. Первый пример (см. рис. 7.3,а) соответствует условиям и, следовательно, в этом случае наблюдается монотонная сходимость метода итераций. В следующем примере (см. рис. 4.3,б) , поэтому имеет место немонотонная сходимость. В последнем примере (см. рис. 4.3,в) . При таких значениях производной метод итераций расходится. Рис. 4.3. Иллюстрация сходимости метода итераций Теорема о сходимости и точности метода итераций. Теорема. Пусть уравнение приведено к виду таким образом, что функция дифференцируема и выполняется условие для всех . Тогда: • последовательность ; • ошибка . Доказательство. Построим два соседних приближения в методе итераций: . Очевидно, что , или . Таким образом, . Запишем это неравенство для k=1,2,…,k в следующем виде: (4.1) Построим вспомогательный ряд (4.2) Частичная сумма членов этого ряда , а значит . По этой причине . Последовательность частичных сумм ряда (4.2) совпадает с последовательностью приближений, вычисленных по методу итераций, а доказательство сходимости вспомогательного ряда эквивалентно доказательству сходимости метода итераций. Рассмотрим вспомогательный числовой ряд (4.3) Ряд (7.3) сходится абсолютно при . Но ряд (7.2) является мажорируемым к ряду (4.3) в силу неравенств (7.1). Следовательно, ряд (7.2) также сходится. Первая часть теоремы доказана. Напоминание. Ряд называется мажорируемым, если каждый его член по модулю не превосходит соответствующего члена некоторого сходящегося ряда с положительными членами. Получим теперь оценку погрешности приближенного решения. На основании леммы имеем: . Найдем . Так как , то . Теперь вычислим. По определению . Тогда . Вторая часть теоремы также доказана. Замечание 1. Величину можно использовать в качестве предельной абсолютной погрешности. Для ее расчета привлекаются два соседних приближения и минимальное по модулю значение первой производной. Предельная же абсолютная погрешность лежит в основе критерия завершения итерационного процесса. Критерий останова: , где – допустимая величина абсолютной ошибки, или . Замечание 2. Рассмотрим теперь, как выбирается при построении функции . Напомним, что . Из теоремы следует, что для сходимости метода итераций необходимо, чтобы . Значит или для. Отсюда можно сделать вывод, что • при ; • при . Лекция 6 МЕТОД НЬЮТОНА ДЛЯ СИСТЕМ И ЕГО МОДИФИКАЦИИ ЦЕЛЬ ЛЕКЦИИ: Изучить метод Ньютона для решения систем нелинейных алгебраических уравнений; обсудить основные модификации метода: метод Ньютона с кусочно-постоянной матрицей, метод Ньютона–Рафсона, методы продолжения по параметру, метод Ньютона для плохо обусловленных систем, дискретный метод Ньютона. Предполагается знакомнство с методом Ньютона из курса информатики. Метод Ньютона для решения систем нелинейных уравнений. Пусть задана система нелинейных уравнений решение которой достигается в точке пространства . Обозначим ; . Тогда исходная система запишется в виде . Предположим, что известно k-е приближение к . Построим правило Ньютона вычисления (k+1)-го приближения в форме . Разложим функцию в ряд Тейлора в окрестности точки и сохраним в разложении два члена: . Полагая, что решение системы достигается на текущей итерации, относительно поправки получим систему линейных алгебраических уравнений: . Тогда , и итерационное правило Ньютона решения системы нелинейных алгебраических уравнений запишется как Такой вид метода Ньютона неудобен на практике, потому что требует вычисления обратной матрицы, а эта операция достаточно трудоемка. На практике метод Ньютона реализуется в следующем виде: 1. Решается система линейных алгебраических уравнений и вычисляется вектор поправки: , где – матрица Якоби системы; 2. Вычисляется (k+1)-е приближение , 3. Пункты 1, 2 повторяются для k=0,1,2,… до тех пор, пока не будет достигнута требуемая точность. Критерием завершения итерационного процесса служат условия , , или в более общей форме , . Модификации метода Ньютона. Краткий обзор. Перечислим недостатки метода Ньютона, на устранение которых направлены различные его модификации: • трудность задания начального приближения, от которого метод сходится; • необходимость вычисления на каждой итерации матрицы Якоби, что может потребовать существенных вычислительных затрат; • необходимость решения на каждой итерации системы линейных алгебраических уравнений; • требование невырожденности матрицы Якоби. Рассмотрим модификации метода Ньютона, которые в той или иной мере устраняют перечисленные недостатки. Метод Ньютона с кусочно-постоянной матрицей. Чтобы уменьшить вычислительные затраты на итерации, матрица Якоби в этой модификации остается постоянной на протяжении нескольких шагов. Число шагов m, на которых J постоянна, задается в такой модификации в качестве параметра, либо момент перевычисления матрицы Якоби определяется условием , в котором , например (матрица Якоби лишь при нарушении этого условия вычисляется заново). Эффективность метода достигается в этом случае не только путем сокращения числа расчетов матрицы Якоби, но главным образом за счет того, что на m итерациях метода приходится решать линейные системы с одной и той же матрицей. Метод Ньютона–Рафсона. Чтобы обеспечить сходимость метода от выбранного начального приближения , применяется модификация, называемая методом Ньютона–Рафсона. Вычисление (k+1)-го приближения в этой модификации осуществляется по правилу , где – параметр, значение которого на k-й итерации выбирается из условия . Стратегия выбора параметра на итерации может быть такой. Вначале принимается пробное значение либо и далее это значение видоизменяется до выполнения сформулированного условия. Это условие может потребовать многократного вычисления вектора на текущей итерации. Очевидно, что при метод Ньютона–Рафсона совпадает с методом Ньютона. Методы продолжения по параметру. Эти методы позволяют обеспечить сходимость метода Ньютона от выбранного начального приближения .Сущность методов продолжения по параметру заключается в замене исходной задачи последовательностью задач, каждая последующая задача при этом незначительно отличается от предыдущей. Последовательность строится таким образом, что первая система имеет решение , а последняя система совпадает с исходной задачей. Поскольку системы отличаются незначительно, то решение предыдущей задачи окажется хорошим начальным приближением для последующей. Решая такую последовательность задач методом Ньютона, получим в итоге решение исходной системы. Рассмотрим способ построения указанной последовательности задач. Пусть при решении системы используется начальное приближение . Заменим исходное уравнение уравнением с параметром , которое при имеет решение , а при совпадает с решением исходной задачи, т. е. . В качестве можно выбрать функции либо . Разобьем отрезок точками на интервалов. Получим искомую последовательность систем: . Метод Ньютона для плохо обусловленных задач. Если матрица Якоби плохо обусловлена, погрешность решения линейной системы может оказаться значительной из-за ошибок округления. Поэтому в случае плохо обусловленных матриц при вычислении вектора поправок привлекают систему с числовым параметром , где E–единичная матрица. При модифицированная система совпадает с линейной системой стандартного метода Ньютона, при стремлении к единице число обусловленности модифицированной системы также стремится к единице. Однако скорость сходимости соответствующего метода значительно ухудшается, поскольку при метод вырождается в метод простых итераций. Дискретный метод Ньютона. Дискретный метод Ньютона базируется на аппроксимации матрицы Якоби на основе вычисленных значений функции в ряде вспомогательных точек. Построим его. Как и прежде, будем использовать векторную запись решаемой системы уравнений . Пусть известно k-е приближение к решению . Аппроксимируем функцию линейной функцией: . Для численного определения матрицы и вектора потребуем, чтобы значения функций и совпадали в (n+1) вспомогательных точках , т. е. чтобы выполнялось равенство . Вычитая из первого равенства все последующие, получим соотношения или в матричной форме , где матрицы и имеют вид Следовательно, . Условие позволяет найти вектор : . Перепишем функцию , подставив найденные соотношения для и : Примем . Будем искать из уравнения Получим итерационную формулу дискретного метода Ньютона: . Заменяя обращение матрицы решением линейной системы, придем к реализуемому на практике алгоритму дискретного метода Ньютона: 1. Вычисляется вектор , матрицы и . 2. Решается система линейных алгебраических уравнений 3. Вычисляется вектор поправки . 4. Вычисляется (k+1)-е приближение 5. Пункты 1÷4 повторяются для k=0,1,2,… до получения решения с требуемой точностью. Применение дискретного метода Ньютона предполагает хранение -матриц и . Однако на практике в качестве вектора выбирается вектор , где – диагональная матрица параметров дискретизации, – j-й столбец единичной матрицы. Элементы матрицы вычисляют по правилу , где – константа (например, 0.1). Часть 4. ИНТЕРПОЛИРОВАНИЕ ФУНКЦИЙ Лекция 7 полиномиальная интерполяция ЦЕЛЬ ЛЕКЦИИ: Пояснить сущность задачи интерполяции; дать общую схему ее решения; построить интерполяционный полином Лагранжа, ввести разделенные разности и получить на их основе интерполяционный многочлен Ньютона. Постановка задачи. Важным элементом численного анализа является задача интерполирования функций, которую приходится решать в следующих приложениях: машинная графика, автоматизация проектирования, обработка экспериментальных данных, управление, передача данных и др. При этом часто возникает задача восстановления функции на отрезке, если известны ее значения в конечном числе точек отрезка. Эти значения определяются на основании экспериментальных измерений либо в результате вычислений. Сформулируем математическую постановку задачи интерполяции. Пусть на отрезке в точках известны значения функции , равные . Требуется построить интерполянту – интерполяционную функцию , совпадающую с функцией в точках : (7.1) и в то же время аппроксимирующую ее, естественно приближенно, на всем отрезке . Основной вопрос интерполяции: как выбрать интерполянту и как оценить погрешность ? Интерполяционные функции строятся, как правило, в виде линейной комбинации некоторых линейно-независимых элементарных функций : , (7.2) где – коэффициенты, которые можно определить, используя условие (7.1). Из этого условия . (7.3) Относительно коэффициентов получили систему линейных алгебраических уравнений, порядок которой равен n+1. В силу линейной независимости элементарных функций определитель этой системы отличен от нуля, и решение для единственно. Таким образом, вычисляя из (7.3) и подставляя их в (7.2), можем получить значение в любой точке . В качестве системы линейно-независимых функций чаще всего выбирают степенные функции, например . В этом случае – полином степени n. Интерполяционная формула Лагранжа. При построении интерполяционной формулы Лагранжа в качестве используются полиномы Лагранжа степени n, удовлетворяющие условиям Для выполнения второго условия полином степени n должен иметь вид , т. е. его корнями являются все узлы интерполяции, кроме k-го. Коэффициент определим, используя первое условие: . Отсюда находим и . Решение системы (7.3) при использовании в качестве элементарных функций полиномов Лагранжа имеет вид . Таким образом, интерполяционный многочлен Лагранжа представится как . Преобразуем его к виду, используемому на практике при вычислении значений функции: . Интерполяционный многочлен совпадает с интерполируемой функцией только в точках . В остальных точках имеет место погрешность интерполяции , которая оценивается величиной , где . (Вывод оценки погрешности опустим). Погрешность интерполяции зависит от числа узлов интерполяции n, а также от их расположения на отрезке . Наилучшими узлами интерполяции следует признать те , для которых принимает наименьшее значение. Недостатком интерполяционной формулы Лагранжа является то, что каждое слагаемое зависит от всех узлов интерполяции. При добавлении узла интерполяции и, следовательно, повышении порядка полинома необходимо вычислять не только слагаемое, относящееся к новому узлу, но и перевычислять заново слагаемые, относящиеся ко всем узлам интерполяции. Понятие конечных и разделенных разностей. Пусть на равномерной сетке в узлах заданы значения функции , т. е. . Конечными разностями первого порядка называют величины . Разности конечных разностей первого порядка называют конечными разностями второго порядка: . Аналогично определяют конечные разности следующих порядков: . Схему вычисления конечных разностей удобно представить в виде следующего треугольника: Если узлы интерполяции – произвольные точки, то вводится понятие разделенных разностей. Разделенная разность первого порядка: . Разделенная разность второго порядка: . Разделенная разность k-го порядка: Отметим два важных свойства разделенных разностей: 1. Разделенная разность является симметричной функцией своих аргументов, т. е. значение разделенной разности не зависит от порядка перечисления узлов интерполяции, а зависит только от их значений; 2. Разделенная разность вычисляется и тогда, когда две или более точки интерполяции одинаковы. Разделенная разность первого порядка в случае, когда одинаковы две точки, имеет вид . Если точка среди узлов интерполяции встречается (m+1) раз, то разделенная разность m-го порядка на этих узлах . Приведем теперь схему вычисления разделенных разностей Интерполяционная формула Ньютона. Пусть . Построим первую разделенную разность . Откуда находим . Построим вторую разделенную разность и выразим из нее : . Подставим в выражение для : . Аналогично привлекаем следующие узлы интерполяции для построения интерполяционной функции. После использования всех узлов интерполяции: Таким образом, интерполяционный многочлен Ньютона имеет вид . Непосредственной проверкой нетрудно убедиться, что . Погрешность интерполяции . Более эффективное вычисление значения функции по интерполяционной формуле Ньютона можно получить, если преобразовать ее к такому виду: Интерполяционная формула Ньютона позволяет легко наращивать число узлов интерполяции, требуя при этом вычисления лишь дополнительных слагаемых. Например, добавление узла приведет к вычислению слагаемого . Лекция 8 сплайн-интерполяция ЦЕЛЬ ЛЕКЦИИ: Показать ограничения полиномиальной интерполяции, ввести понятие m-сплайна, построить вычислительную схему интерполяции на основе кубического сплайна. Ограничения многочленной интерполяции. Многочленная интерполяционная функция очень чувствительна к выбору узлов интерполяции. Рассмотрим интерполяционную формулу Лагранжа , где . Оценим норму функции , которую определим, как . С этой целью выполним очевидные преобразования: . Введем функцию Лебега: . Поскольку , получаем оценку: . Величина нормы функции зависит от распределения узлов интерполяции. Рассмотрим два случая. Случай 1. Узлы интерполяции на отрезке распределены равномерно. В этом случае (доказательство опускаем). При увеличении n многочлен может полностью оказаться непригодным для аппроксимации , т. к. неограниченно возрастает. Рассмотрим такой пример. Будем интерполировать функцию полиномом n-го порядка, выбирая узлы интерполяции равномерно распределенными. Введем норму ошибки интерполяции: и исследуем ее зависимость от порядка интерполирующего полинома. Для этого обратимся к численным результатам (см. табл.8.1). Таблица 8.1 2 0.96 8 0.25 14 1.07 4 0.71 10 0.30 16 2.10 6 0.43 12 0.56 18 4.21 Видно, что с увеличением погрешность интерполяции уменьшается вплоть до . Дальнейшее увеличение порядка полинома приводит к возрастанию . Случай 2. Узлы интерполяции являются нулями полинома Чебышева. Их нетрудно получить, если на отрезке построить полуокружность, разделить ее на равные части и спроектировать на отрезок середину каждой из них (рис. 8.1). В этом случае (доказательство опускаем). Обратимся к рассмотренному выше примеру (см. табл.8.2) Погрешность интерполяции теперь при возрастании порядка полинома монотонно падает. Рис. 8.1. Расположение нулей полинома Чебышева Таблица 8.2 2 0.93 8 0.39 14 0.12 4 0.75 10 0.27 16 0.08 6 0.56 12 0.18 18 0.06 Многочленная интерполяция в точках Чебышева позволяет увеличивать точность приближения функции посредством увеличения порядка полинома. Однако привлекать узлы Чебышева не всегда удается. (Например, в узле Чебышева функция имеет особенность.) Чтобы избежать зависимости точности аппроксимации от локальных свойств функции, переходят к кусочно-полиномиальной аппроксимации. Сплайн-интерполяция. Назовем m-сплайнами полиномы невысокого порядка m, которыми аппроксимируется функция на интервалах и которые сшиваются в узлах интерполяции на основе требования непрерывности интерполирующей функции и ее первых производных. Рассмотрим наиболее широко используемую на практике кубическую сплайн-интерполяцию. В ней искомая функция на интервале аппроксимируется полиномом третьей степени (8.1) Аналогичные полиномы можно записать для всех n интервалов, т. е. соотношение (8.1) справедливо для . Условия сшивки сплайнов в узлах интерполяции: , (8.2) , (8.3) . (8.4) Кроме того, необходимо задать граничные условия в точках и . Это можно сделать несколькими способами. Здесь будем считать, что , т. е. будем рассматривать интерполяцию так называемыми естественными сплайнами. Условие (8.2) приводит к уравнениям: (8.5) Получили уравнений относительно неизвестных , . Вычислим производные: Построим уравнения, к которым приводит условие непрерывности первой производной в узлах . При кубической сплайн-интерполяции выражения для первой производной на соседних интервалах интерполирования имеют вид Требование непрерывности первой производной в узлах (условие 8.3) приводит к уравнениям . (8.6) В свою очередь, непрерывность второй производной в этих же узлах (условие (8.4)) позволяет записать , или . (8.7) Наконец, из граничных условий , получим, что . (8.8) Соотношения (8.5), (8.6), (8.7), (8.8) составляют систему линейных алгебраических уравнений, всего уравнений, относительно неизвестных. Построим эффективный метод решения этой системы. Выразим из (8.7) : , (8.9) и подставим в (8.5): . Учтем, что : . Выразим из этого соотношения : . (8.10) Подставим теперь в формулу (8.6): Выполним очевидные преобразования и запишем результирующую систему в виде: (8.11) Получили систему линейных алгебраических уравнений относительно неизвестных . Однако, если учесть, что из граничных условий , то приходим к системе линейных алгебраических уравнений с треугольной симметрической матрицей Обратите внимание: матрица системы обладает диагональным преобладанием! Эту систему можно эффективно решать методом -факторизации. Вычислив коэффициенты , из соотношений (8.9) находим . Из условия (8.8) определяем : Коэффициенты рассчитаем из (8.10): так как . Коэффициенты , известны из (8.5). В результате для всех кубических сплайнов определены коэффициенты . Расчет значений функции методом сплайн-интерполяции осуществляется следующим образом: 1. По описанной выше методике вычисляются коэффициенты кубических сплайнов; 2. Находится интервал , которому принадлежит данное . Значение функции на этом интервале вычисляется из кубического сплайна с параметрами для интервала . Часть 5. определенные интегралы Лекция 9 методы численного интегрирования ЦЕЛЬ ЛЕКЦИИ: Дать схему построения методов численного интегрирования посредством замены подынтегральной функции интерполяционным полиномом; получить на этой основе численные методы вычисления одномерных интегралов (методы трапеций и Симпсона); распространить метод Симпсона на случай вычисления двойных интегралов; показать, как строится схема вычисления многомерных интегралов при использовании метода Монте–Карло. Постановка задачи. Речь идет о вычислении определенного интеграла в тех случаях, когда функция задана таблично, либо когда первообразная функции находится очень сложно. Общая схема построения рассматриваемых вычислительных методов расчета определенных интегралов сводится к следующему. 1. Отрезок интегрирования покрывается сеткой и вычисляются значения функции во всех узлах сетки (рис. 9.1). Рис. 9.1. Покрытие отрезка интегрирования сеткой 2. Подынтегральная функция на всем отрезке или на его отдельных частях заменяется легко интегрируемой интерполяционной функцией , для которой . В качестве интерполяционной функции чаще всего используется по- линомиальная функция, т. е. 3. Вычисляется приближенное значение определенного интеграла и оценка погрешности . Формула трапеций. Примем шаг сетки, которой покрывается отрезок интегрирования , постоянным и равным , где – число интервалов разбиения. Узлы интерполяции вычисляются в этом случае по правилу: , причем . Заменим подынтегральную функцию на отрезке интерполяционным многочленом Лагранжа первого порядка: . Тогда . Первый член в таком представлении является приближенным значением определенного интеграла, второй – погрешностью расчета. В частности Аналогично Просуммируем , где Отсюда следует, что формула трапеций имеет достаточно простой вид: . При этом ошибка вычисления определенного интеграла не превосходит величины: , где . Погрешность формулы трапеций имеет второй порядок относительно шага сетки. Точность вычисления определенного интеграла может быть повышена двумя способами: уменьшением шага сетки, увеличением точности интерполяции подынтегральной функции. Формула Симпсона. Разобьем интервал интегрирования на четное число частей . В этом случае шаг сетки (шаг интерполирования) , и сеточные узлы принимают значения , при этом . Рассмотрим простейший случай, когда сетка содержит только три узла: . Очевидно, что. Вычислим приближенное значение интеграла, заменяя подынтегральную функцию интерполяционным многочленом Лагранжа второго порядка: Первый член, составляющий приближенное значение искомого интеграла, легко интегрируется точно. В результате имеем следующее равенство: , где , – погрешность вычисления интеграла. Таким образом, формула Симпсона для случая трех узлов имеет вид . (9.1) Оценим погрешность . Так как погрешность интерполяции подынтегральной функции многочленом Лагранжа второго порядка пропорциональна третьей производной от подынтегральной функции, то соотношение (9.1) является точным для всех подынтегральных функций, описываемых полиномом второй степени. В силу симметрии эта формула является точной и для подынтегральных функций, описываемых полиномом третьей степени: , так как она точна для . В этом нетрудно убедиться, проверив справедливость равенства . Рассмотрим полином третьей степени, удовлетворяющий условиям . Он интерполирует функцию на отрезке по значениям функции в узлах и по значению ее производной в узле (узел имеет кратность два): , где – погрешность кратной интерполяции, которая равна . (Вывод погрешности кратной интерполяции опустим.) Тогда можем записать, что Найдем теперь погрешность приближения определенного интеграла: Пусть теперь сетка содержит произвольное число узлов , принадлежащих отрезку интегрирования , причем. Последовательно вычислим интеграл на отрезках длиной , интерполируя подынтегральную функцию многочленом Лагранжа второго порядка: Просуммируем левые и правые части этих соотношений: , где . Таким образом, формула Симпсона принимает следующий вид: . Погрешность вычисления определенного интеграла по формуле Симпсона имеет четвертый порядок относительно шага сетки: где . Формула Симпсона для вычисления двойных интегралов. Вычислить . Область интегрирования здесь – прямоугольник со сторонами и (рис. 9.2). Применим формулу Симпсона, вычисляя определенный интеграл сначала по , затем по . Рис. 9.2. Область интегрирования Пусть , и . Тогда Для повышения точности вычислений область покрывается сетью прямоугольников. В этом случае , и значение двойного интеграла вычисляется в виде где . Если – криволинейная область, то для применения полученной формулы Симпсона область заключают в прямоугольник и пользуются вспомогательной функцией Тогда и для вычисления последнего интеграла привлекают метод Симпсона. Вычисление многомерных интегралов методом Монте–Карло. Пусть необходимо вычислить . Заключим область интегрирования внутрь -мерного параллепипеда со сторонами , т. е. . Сделаем замену переменных . Тогда -мерный параллепипед преобразуется в -мерный единичный куб, т. к. . Область преобразуется в область , заключенную внутри -мерного единичного куба (рис. 9.3). Рис. 9.3. Преобразование области интегрирования в двумерном случае С учетом преобразования переменных где . По теореме о среднем можно положить , где – объем области интегрирования , – усредненное значение функции в области . Для вычисления и воспользуемся методом статистических испытаний (методом Монте–Карло). Пусть мы умеем строить случайные числа, распределенные в интервале по равномерному закону. Обозначим через случайную точку -мерного пространства, координаты которой являются независимыми случайными величинами, распределенными в интервале по равномерному закону. Используя датчик случайных чисел, сформируем случайных точек . Разобьем это множество точек на два подмножества: Пусть подмножество содержит элементов, т. е. точек из общего числа принадлежат области интегрирования . Тогда Следовательно, . Окончательное расчетное соотношение метода Монте–Карло для вычисления определенных интегралов принимает вид . Методом Монте-Карло пользуются в тех случаях, когда достаточен невысокий порядок точности. Часть 6. системы обыкновенных дифференциальных уравнений Лекция 10 ЧИСЛЕННОЕ РЕШЕНИЕ СИСТЕМ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ ЦЕЛЬ ЛЕКЦИИ: Дать постановку задачи численного решения системы обыкновенных дифференциальных уравнений с начальными условиями; провести классификацию численных методов по виду разностной схемы; ввести понятия локальной и интегральной точности, устойчивости по отношению к шагу интегрирования; построить простейший явный одношаговый метод, оценить его локальную погрешность, устойчивость к шагу интегрирования; пояснить свойство жесткости дифференциальных уравнений и непригодность для их решения классических явных методов. Общая характеристика численных методов. Пусть требуется решить систему обыкновенных дифференциальных уравнений (СОДУ) с начальными условиями Это значит, что необходимо найти функции на отрезке такие, что Запишем сформулированную задачу, ее называют задачей Коши для системы обыкновенных дифференциальных уравнений, в векторном виде: Решение задачи Коши заключается в нахождении такой траектории , которая бы удовлетворяла условию . Применение численных методов для решения задачи Коши предполагает приближенное вычисление значений в точках отрезка . С этой целью СОДУ заменяют разностной схемой, из которой рекуррентно вычисляют приближенные значения . Запишем общий вид разностной схемы: . Она связывает искомое решение в текущий момент времени с построенными к этому моменту времени решениями в предыдущие моменты времени. Здесь – приближенное значение , т. е. Выбор функции в разностной схеме определяет соответствующий численный метод. Приведем общую характеристику методов. • Если , т. е. , то метод является одношаговым. Он связывает решение в последующий момент времени с решением в предыдущий момент времени. Если , то метод является многошаговым. • Метод является явным, если функция не зависит от и неявным в противном случае. В случае явного одношагового метода искомое решение на текущем шаге интегрирования вычисляется тривиально: . Неявный одношаговый метод требует для расчета решать систему в общем случае нелинейных алгебраических уравнений , привлекая, например, итерационную процедуру Ньютона. • Различают локальную и интегральную погрешности метода. Локальная погрешность – ошибка на шаге интегрирования, интегральная погрешность – полная погрешность в текущий момент времени. Обратимся к графической иллюстрации погрешности. Пусть численно решается уравнение . В результате выполнения первого шага интегрирования (см. рис. 10.1) приближенное значение искомой функции в момент времени равно Рис. 10.1. Иллюстрация локальной и интегральной погрешностей . Это значение принадлежит некоторой интегральной кривой 2, являющейся решением дифференциального уравнения при другом, отличном от , начальном условии. Разница между приближенным значением и точным составляет локальную погрешность. Очевидно, что на первом шаге локальная погрешность совпадает с интегральной погрешностью. Второй шаг интегрирования приводит к приближенному значению , которое принадлежит интегральной кривой 3. Видно, что локальная погрешность на этом шаге, равная разности и значения ординаты кривой 2 в момент времени , существенно отличается от интегральной погрешности – разности между приближенным и точным значениями. Наиболее объективной характеристикой точности метода является величина интегральной погрешности. К сожалению, выполнить ее оценку крайне сложно. На практике при выборе шага интегрирования обычно используют локальную погрешность численного метода. • Если величина шага интегрирования ограничивается только допустимой локальной погрешностью метода, то такой метод является абсолютно устойчивым к шагу интегрирования. Метод, обладающий таким свойством, позволяет осуществлять выбор шага интегрирования, исходя лишь из требуемой точности. Условно устойчивый метод характеризуется тем, что шаг интегрирования ограничивается не только допустимой локальной погрешностью, но и определенными свойствами решаемой системы. Последние ограничения являются весьма обременительными для некоторых классов задач. Рассмотрим разностные схемы наиболее широко используемых на практике методов численного интегрирования. Явный метод Эйлера. Рассмотрим , где , – текущий шаг интегрирования. Разложим функцию в ряд Тейлора в окрестности точки : . Ограничившись в этом разложении двумя членами, получим разност-ную схему метода Эйлера . Локальная погрешность метода Эйлера составляет величину . В вычислительной математике численные методы решения обык-новенных дифференциальных уравнений принято характеризовать порядком точности. Определение. Если локальная погрешность численного метода ,то порядок точности такого метода равен . Метод Эйлера является методом первого порядка. Приведем геометрическую интерпретацию явного метода Эйлера для задачи Коши (см. рис. 10.2). Приращение на шаге интегрирования – катет прямоугольного треугольника, лежащий против угла, тангенс которого равен значению производной в предыдущий момент времени. Вторым катетом этого треугольника является текущий шаг интегрирования. Рис. 10.2. Геометрическая иллюстрация явного метода Эйлера Оценим устойчивость метода Эйлера по отношению к шагу ин- тегрирования. Для этого рассмотрим линейную автономную систему с отрицательно определенной матрицей простой структуры. Отрицательная определенность матрицы означает, что все собствен-ные значения матрицы действительны и отрицательны, т. е. . В этом случае все решения . Применим для решения этой системы метод Эйлера с постоянным шагом : . Здесь E – единичная матрица соответствующей размерности. Из алгебры известно, что для любой неособенной матрицы простой структуры существует такая неособенная матрица , которая преобразованием подобия приводит матрицу к диагональному виду: . Преобразуем вычислительную схему метода Эйлера следующим образом: . Введем замену переменных . Тогда , или . Запишем это соотношение для i-й компоненты вектора : , или , где, т. е. определяется начальным условием. Нетрудно видеть, что , если . Именно этим свойством обладает решение автономной системы с отрицательно определенной матрицей. Отсюда приходим к требованиям , при этом неравенство приводит к естественному условию , т. к. , а неравенство  к условию . Очевидно, чтобы , необходимо при выборе шага интегрирования выполнить условие . Таким образом, явный метод Эйлера по отношению к шагу интегрирования является условно устойчивым. Явление жесткости. Ограниченная устойчивость численного метода является серьезным недостатком при решении так называемых жестких систем. Рассмотрим дифференциальное уравнение первого порядка , для которого уравнение имеет единственное решение Рис. 10.3. К определению жесткости дифференциального уравнения . Любая интегральная кривая такого дифференциального уравнения характеризуется двумя участками с существенно различным поведением решения (см. рис. 10.3), причем первый участок значительно меньше второго. Первый участок с быстрым изменением функции отражает стремление интегральной кривой к графику функции и называется пограничным слоем. На втором участке интегральная кривая практически совпадает с графиком . Однако даже при небольшом отклонении от графика в любой его точке производная резко возрастает по сравнению с производной . Именно по этой причине малый шаг интегрирования, используемый при воспроизведении быстропротекающего участка, не может быть существенно увеличен вне пограничного слоя в явных методах численного интегрирования. Рассмотрим примеры жестких систем. 1. Простейшее дифференциальное уравнение (10.1) может быть жестким (рис. 10.4) , если интервал наблюдения зна- Рис. 10.4. Семейство решений уравнения (10.1) чительно превосходит величину . Решением такого уравнения является функция . Пунктирные линии на рис. 10.4 соответствуют различным значениям . 2. Для иллюстрации явления жесткости дифференциальных уравнений может быть полезной также любая система линейных ОДУ, матрица которой характеризуется большим разбросом собственных значений, например система (10.2) Здесь собственные числа матрицы есть . Решение этой системы (рис. 10.5) содержит быструю составляющую как для , так и для , хотя по амплитуде быстрая составляющая функции значительно превосходит аналогичную компоненту функции (ее на рис. 10.5 в таком масштабе отобразить не удалось). Рис. 10.5. Решение жесткой системы ОДУ (10.2) Завершая краткую характеристику свойства жесткости дифференциальных уравнений, отметим, что при решении научных задач жесткость уравнений является скорее правилом, чем исключением. По этой причине для таких задач разработаны специальные методы численного интегрирования. Лекция 11 одношаговые методы ЦЕЛЬ ЛЕКЦИИ: Построить простейший неявный одношаговый метод, обладающий свойством устойчивости к шагу интегрирования, оценить его локальную погрешность; дать способ Рунге–Кутта увеличения точности одношаговых методов, проиллюстрировав его методами второго и четвертого порядков точности; привести разностную схему линейных многошаговых методов, получить условия корректного выбора коэффициентов. Неявный метод Эйлера. Пусть требуется численно решить систему обыкновенных дифференциальных уравнений . Формально неявный метод Эйлера можно получить, рассматривая , где – шаг интегрирования. Разложим в ряд Тейлора в окрестности точки : . Ограничившись в разложении двумя членами, придем к разностной схеме неявного метода Эйлера: . Локальная погрешность при этом определяется отброшенными членами ряда Тейлора: . Сравнивая явный и неявный методы Эйлера между собой (см. рис. 11.1), следует отметить, что методы обладают близкой по модулю, но разной по знаку погрешностью. Рассмотрим устойчивость неявного метода Эйлера по отношению к шагу интегрирования. Применим его к системе уравнений с отрицательно определенной матрицей , полагая шаг интегрирования постоянным: . Отсюда Пусть – неособенная матрица, которая преобразованием подобия приводит матрицу к диагональному виду . Привлекая матрицу , преобразуем итерационное правило следующим образом: , или , где новая переменная . Запишем результат для -й компоненты вектора : . Отсюда следует, что при любом , поскольку все . Неявный метод Эйлера является абсолютно устойчивым по отношению к шагу интегрирования. При решении этим методом жестких систем дифференциальных уравнений шаг интегрирования выбирается только из соображений допустимой локальной погрешности. Методы Рунге–Кутта. Точность явных одношаговых методов численного решения систем обыкновенных дифференциальных уравнений вида можно повысить, сохраняя в разложении функции в ряд Тейлора большее число членов. Например, метод второго порядка имеет следующую разностную схему: , или где . Основное неудобство такой формы разностной схемы – необходимость вычисления частных производных , . Эта трудность значительно возрастает при построении методов более высокого порядка точности. В методах Рунге–Кутта функции , где – порядок точности метода, заменяются на некоторые удобно вычисляемые функции таким образом, что , где – константа, не зависящая от . В методе Рунге–Кутта второго порядка функция имеет вид . Разложим функцию в ряд Тейлора в окрестности точки : . Подставим это разложение в выражение для : Сравнивая и , нетрудно видеть, что при эти формы совпадают с точностью до члена . Если положить , то , . В результате метод Рунге–Кутта второго порядка примет вид: , где По аналогии можно построить методы Рунге–Кутта более высоких порядков. Не останавливаясь на выводе, приведем популярный на практике метод Рунге–Кутта четвертого порядка: , где На каждом шаге интегрирования в методе Рунге–Кутта четвертого порядка приходится четырежды вычислять значение функции при разных значениях аргументов. Более того, эти значения функции используются лишь однократно, что отражается на эффективности вычислений. Методы Рунге–Кутта относятся к классу явных условно устойчивых методов. По этой причине они оказываются неприемлемыми для решения жестких систем обыкновенных дифференциальных уравнений. Линейные многошаговые методы. Пусть требуется найти решение на отрезке задачи Коши . Предположим, что построены приближенные значения решения и его первой производной в моменты времени , т. е. и . Общий вид разностной схемы рассматриваемых здесь многошаговых методов имеет вид где – коэффициенты (их всего), которые должны быть определены при получении конкретного многошагового метода, – шаг интегрирования. Значения этих коэффициентов выбирают так, что если решение является полиномом степени , то разностная схема многошагового метода дает точное значение, т. е. . Поскольку полином степени имеет параметр, то разностная схема должна иметь по крайней мере коэффициент. В большинстве практических многошаговых методов и лишние коэффициенты могут быть выбраны произвольно. Получим соотношения, которым должны удовлетворять все коэффициента разностной схемы в предположении, что метод дает точное решение для задачи Коши, точным решением которой является полином степени . Поскольку полином -й степени включает в себя все полиномы степени ниже , то разностная схема должна также давать точное решение для всех задач Коши, имеющих полиномиальное решение степени меньшей, чем . В частности: 1. , . Класс задач с таким решением задается уравнением . Поэтому , , . Подставив эти значения в разностную схему, получим первое условие, которому должны удовлетворять коэффициенты : . 2. , . Класс задач с таким решением задается в виде . Для удобства выберем . Тогда , , , . В этом случае . Подставим их в разностную схему: . Если учесть первое условие, то это соотношение преобразуется к виду . Наконец, разделив левую и правую части на , получим условие корректности для полиномиальных решений первой степени: . 3. , . Класс задач с таким решением . Полагаем, как и прежде, и находим, что Перепишем с учетом этих соотношений разностную схему многошагового метода или Разделим левую и правую части этого соотношения на . Условие корректности для полиномиальных решений второй степени примет следующий вид: . 4. Общий случай: . Класс задач с таким решением . Условие корректности для полиномиальных решений степени : . Анализ выписанных условий корректности для полиномиальных решений до степени включительно свидетельствует, что они имеют одинаковую форму, а именно: Этим соотношениям должны удовлетворять все коэффициента разностной схемы линейного многошагового метода. Лекция 12 многошаговые методы ЦЕЛЬ ЛЕКЦИИ: На основе общей разностной схемы линейных многошаговых методов и условий корректного выбора коэффициентов построить явные многошаговые методы Адамса, неявные многошаговые методы Адамса, неявные многошаговые методы Гира различных порядков точности. Явные многошаговые методы Адамса. Этот класс методов легко получить из общей разностной схемы линейных многошаговых методов и условий корректного выбора коэффициентов Полагая , , получим разностную схему явного многошагового метода порядка , где коэффициенты , их всего , однозначно определяются из условия корректности. Из первого условия корректности находим . Обозначим . Эта запись подчеркивает, что значение коэффи-циента зависит от порядка метода. Чтобы определить , подставим значения в оставшиеся условий корректности: Получили систему из линейных алгебраических уравнений относи-тельно : Запишем явные методы Адамса различных порядков. 1. . Коэффициент и . Разностная схема явного метода Адамса первого порядка совпадает с разностной схемой явного метода Эйлера. 2. . Соотношения для являются системой из двух линейных алгебраических уравнений Следовательно , и разностная схема явного метода Адамса второго порядка принимает вид . 3. . Относительно , запишем систему линейных алгебраических уравнений третьего порядка , из которой следует, что . В результате разностная схема явного метода Адамса третьего порядка записывается следующим образом: . Этот процесс записи разностных схем явных многошаговых методов более высоких порядков точности можно продолжить и далее. Не останавливаясь на выводе, приведем выражение для локальной погрешности явных многошаговых методов Адамса порядка : , где – постоянный коэффициент, величина которого зависит от метода, в частности – значение -й производной функции в точке . Неявные многошаговые методы Адамса. Неявный метод Адамса -го порядка получается из общей разностной схемы линейных многошаговых методов при условии , т. е. . Здесь , их всего , определяются из условий корректности полиномиальных решений -го порядка. Как и в явных методах Адамса, . Обозначим , , подчеркивая тем самым зависимость значений коэффициентов от порядка метода. Подставим выбранные значения параметров в условие корректности: . Относительно неизвестных коэффициентов получим систему из линейных алгебраических уравнений: . Решение этой системы однозначно определяет коэффициенты неявного метода Адамса -го порядка. Запишем неявные методы Адамса первого, второго и третьего порядков. 1. . Коэффициент и . Разностная схема неявного метода Адамса первого порядка совпадает с неявным методом Эйлера. 2. . В этом случае , и искомые коэффициенты имеют следующие значения: , . Неявный метод Адамса второго порядка , как легко видеть, является методом трапеций. 3. . Коэффициенты вычисляются в этом случае из системы . Ее решение: . Неявный метод Адамса третьего порядка принимает вид: . Запись неявных методов Адамса более высоких порядков точности можно продолжить по аналогии. Приведем теперь без вывода соотношение для расчета локальной погрешности неявного метода Адамса -го порядка: , где значения коэффициента для рассмотренных выше разностных схем соответственно равны: . Сравнивая явные и неявные методы Адамса одного порядка между собой, можно отметить, что методы имеют одинаковую по порядку, но разную по знаку локальную погрешность. Неявные многошаговые методы Гира. Неявный многошаговый метод Гира -го порядка получается из общей разностной схемы многошаговых методов при следующем выборе параметров: . Его разностная схема . Коэффициенты должны быть определены таким образом, чтобы выполнялись условия корректности полиномиальных решений, которые, с учетом выбранных значений , принимают вид Запись подчеркивает зависимость значений этих коэффициентов от порядка метода. Запишем условия корректности полиномиальных решений многошаговых методов Гира в развернутом виде: . Решение этой системы линейных алгебраических уравнений единственным образом определяет коэффициенты метода Гира -го порядка. Приведем методы Гира первого, второго и третьего порядков. 1. . Коэффициенты определяются системой , из которой находим , и разностная схема метода Гира записывается следующим образом: . Разностная схема метода Гира первого порядка совпадает с неявным методом Эйлера. 2. . В этом случае неизвестные коэффициенты , являются решением системы . Видно, что . Разностная схема метода Гира второго порядка имеет вид . 3. . При этом система линейных алгебраических уравнений определяет следующие значения неизвестных коэффициентов: . Подставив эти коэффициенты в общую формулу методов Гира, получим метод Гира третьего порядка: . По аналогии можно построить многошаговые методы Гира и более высоких порядков. Приведем в заключение оценку локальной погрешности (без вывода) метода Гира -го порядка: , где константа зависит от порядка метода и для приведенных схем соответственно равна: . Лекция 13 устойчивость многошаговых методов ЦЕЛЬ ЛЕКЦИИ: Дать общую схему исследования устойчивости линейных многошаговых методов к шагу интегрирования; построить области устойчивости явных и неявных методов Адамса, неявных методов Гира, порядок точности которых изменяется в диапазоне 16; ввести понятия -устойчивости, жесткой устойчивости; сравнить методы по критерию устойчивости; показать, как осуществляются начальные вычисления в многошаговых методах. Выбор тестовой задачи. При анализе устойчивости многошаговых методов обращаются к простейшей задаче Коши , где – в общем случае комплексная величина. Такой выбор тестовой задачи объясняется тем, что устойчивая линейная система (ее решение стремится к нулю) общего вида с действительной -матрицей простой структуры преобразованием приводится к виду , где – неособенная матрица преобразования, – собственные значения матрицы (в общем случае комплексные величины). Таким образом, каждое уравнение преобразованной устойчивой линейной системы является простейшим дифференциальным уравнением, совпадающим по форме с тестовым уравнением привлекаемой задачи Коши. Точное решение тестовой задачи имеет вид . Если , то при любом значении . Применяя численный метод с постоянным шагом к тестовому уравнению, получим последовательность . Если , то такой метод является устойчивым. Условие устойчивости линейного многошагового метода. Будем полагать теперь, что для решения тестового уравнения используется многошаговый метод. В этом случае . Преобразуем это выражение следующим образом: , где .Введем подстановку . Тогда разностное уравнение приводится к виду , или . Каждый корень полиномиального уравнения , их всего для каждого фиксированного , порождает частное решение разностного уравнения многошагового метода. Так как разностное уравнение многошагового метода линейно и для него справедлив принцип суперпозиции, общее решение в предположении, что все корни различны, можно представить в виде , где – постоянные коэффициенты. Если полиномиальное уравнение содержит кратный корень кратности , то соответствующий член в общем решении будет таким : . Отсюда следует, что , если . Таким образом, многошаговый метод является численно устойчивым для тех значений , для которых корень полиномиального уравнения лежит внутри единичной окружности . Множество всех значений , для которых многошаговый метод является численно устойчивым, называют областью устойчивости метода. Построение области устойчивости. Рассмотрим, как строится область устойчивости многошагового метода. Преобразуем полиномиальное уравнение следующим образом: , где Будем считать, что и вычислим соответствующее . Так как любое на единичной окружности в комплексной плоскости можно выразить в виде то представляет множество точек границы области устойчивости на комплексной плоскости (геометрическая иллюстрация такого преобразования приведена на рис. 13.1). Нетрудно проверить, что кривая симметрична относительно оси . Поэтому преобразование выполняют только для углов . Границу области устойчивости метода для углов восстанавливают путем симметричного Рис. 13.1. Геометрическая иллюстрация построения области устойчивости в плоскости σ отображения кривой относительно оси . Покажем теперь области устойчивости многошаговых методов различных порядков точности (от первого до шестого), построенные в плоскости по описанной методике (см. рис. 13.2, 13.3, 13.4). Из рисунков следует, что размер области устойчивости линейных многошаговых методов (как явных, так и неявных) уменьшается с увеличением порядка точности. Наибольшая область устойчивости характерна для методов первого порядка точности, наименьшая – для методов шестого порядка точности. Условие устойчивости явных многошаговых методов Адамса, впрочем, как и неявных методов Адамса выше второго порядка точности, накладывает существенные ограничения на величину шага интегрирования. По этой причине эти методы не подходят для решения жестких задач. Жесткие дифференциальные уравнения целесообразно ин-тегрировать неявными методами Адамса первого либо второго порядков точности, а также неявными методами Гира. Рис. 13.2. Области устойчивости явных многошаговых методов Адамса Рис. 13.3. Области устойчивости неявных многошаговых методов Адамса Рис. 13.4. Области устойчивости неявных многошаговых методов Гира Определение 1. Метод называется -устойчивым, если его область устойчивости содержит всю правую полуплоскость комплексной плоскости . Свойство -устойчивости метода позволяет выбирать шаг интегрирования исходя только из требования точности. -устойчивыми являются неявные методы Адамса и Гира первого и второго порядков точности. Определение 2. Метод называется жестко устойчивым, если его область устойчивости содержит две подобласти и (см. рис. 13.5) комплексной плоскости , где определяется условиями , – условием , причем в области метод обеспечивает требуемую точность. Сущность требования заключается в том, чтобы при малых , когда , обеспечить требуемую точность быстрых составляющих решения, а при больших гарантировать затухание быстрых составляющих, не учитывая погрешность их воспроизведения. Жестко устойчивые методы численного интегрирования позволяют эффективно решать жесткие системы обыкновенных дифференциальных уравнений, так как шаг интегрирования в этих методах ограничивается только условием допустимой погрешности. Анализ областей устойчивости методов Гира от третьего до шестого порядков включительно свидетельствует о том, что эти методы относятся к классу жестко устойчивых (методы Гира более высоких порядков не являются жестко устойчивыми). В частности, для метода третьего порядка , для четвертого порядка , для пятого порядка , для метода шестого порядка . Начальные вычисления в многошаговых методах. В противоположность одношаговым методам многошаговые методы численного интегрирования не являются самоначинающимися. Так, явный метод Адамса -го порядка является -шаговым. Для вычисления на его основе необходимо знание решения в предыдущих временных точках, т. е. ..., . В начале вычислений известно только . По этой причине первый шаг интегрирования может быть осуществлен только одношаговым методом, второй шаг – одношаговым либо двухшаговым методом, третий шаг – одношаговым, двухшаговым либо трехшаговым методом и т. д. Такая же стратегия вычислений реализуется и при использовании неявных -шаговых методов Адамса и Гира. Лекция 14 изменение шага в многошаговых методах ЦЕЛЬ ЛЕКЦИИ: Показать, как можно на основе так называемого вектора Нордсика, не прибегая к задаче интерполяции, эффективно изменять шаг интегрирования в процессе численного решения обыкновенных дифференциальных уравнений. Сущность проблемы изменения шага в многошаговых методах сводится к следующему. В текущий момент времени в памяти компьютера хранятся вычисленные значения решения и умноженные на значения его производной в моменты времени . Обозначим их вектором . Если в момент времени принимается решение об изменении шага интегрирования, т. е. , где , то для выполнения шага интегрирования многошаговым методом необходимы предшествующие значения решения и значения его производной в моменты времени . Эти значения образуют вектор , причем только . Остальные компоненты вектора должны быть перевычислены. Использование для этой цели интерполяционной процедуры может потребовать значительных вычислительных усилий. Имеется другой путь. Введем так называемый вектор Нордсика, который для значений временного шага и имеет вид , . Заметим, что , так как эти величины относятся к одному и тому же моменту времени . Следовательно, значения векторов и связаны простым соотношением , или в компактной форме , где – диагональная матрица. С другой стороны, вектор Нордсика может быть определен посредством следующего линейного преобразования переменных: , где – невырожденная матрица преобразования переменных линейного многошагового метода. Объединение введенных матричных преобразований приводит к линейной системе , из которой для заданного вектора можно найти искомый вектор . Так как матрица имеет невысокий порядок, то такая процедура перехода от к оказывается достаточно эффективной. Рассмотрим, как строится матрица преобразования переменных многошагового метода на примере явного метода Адамса третьего порядка. В этом случае и . При построении матрицы воспользуемся тем свойством, что явный метод Адамса третьего порядка позволяет найти точное решение, описываемое полиномом третьей степени . Очевидно, что . Полагая величину постоянной на трех шагах, выберем , , и вычислим из соотношений для , зна-чения : . Решим эту систему относительно коэффициентов : Подставляя значения коэффициентов в соотношения для , нетрудно найти, что В итоге связь вектора Нордсика явного метода Адамса третьего порядка с вектором имеет следующий вид: . Это преобразование переменных имеет место также для неявного метода Адамса четвертого порядка. Действительно, неявный метод Адамса четвертого порядка является трехшаговым методом и для него векторы и такие же, как и для явного метода Адамса третьего порядка. Наконец, в случае метода Гира третьего порядка (он является трехшаговым) векторы и соответственно равны Действуя по аналогии, нетрудно получить, что . В заключение подчеркнем, что для любого многошагового метода можно построить матрицу преобразования переменных. В свою очередь, привлечение матриц и позволяет эффективно перевычислять вектор при изменении шага интегрирования. Сделаем ряд практических замечаний по проблеме выбора шага. Замечание 1. При оценке погрешности интегрирования многошагового метода порядка по формуле необходимо вычислить значение . Это легко сделать, используя -ю компоненту вектора Нордсика на двух соседних шагах: . Погрешность расчета находится следующим образом: . Замечание 2. При изменении шага интегрирования локальная погрешность оценивается по формуле . Из этого соотношения легко найти параметр : , где – допустимая величина локальной погрешности.