Математическое моделирование и вычислительный эксперимент

Введение Математическое моделирование и вычислительный эксперимент 1. Схема вычислительного эксперимента. Эффективное решение крупных естественнонаучных и народнохозяйственных задач сейчас невозможно без применения быстродействующих электронно-вычислительных машин (ЭВМ). В настоящее время выработалась технология исследования сложных проблем, основанная на построении и анализе с помощью ЭВМ математических моделей изучаемого объекта. Такой метод исследования называют вычислительным экспериментом. Пусть, например, требуется исследовать какой-то физический объект, явление, процесс. Тогда схема вычислительного эксперимента выглядит так, как показано на рисунке 1. Формулируются основные законы, управляющие данным объектом исследования (I) и строится соответствующая математическая модель (II), представляющая обычно запись этих законов в форме системы уравнений (алгебраических, дифференциальных, интегральных и т. д.). Рисунок 1 – Этапы построения и анализа с помощью ЭВМ математической модели объекта При выборе физической и, следовательно, математической модели мы пренебрегаем факторами, не оказывающими существенного влияния на ход изучаемого процесса. Типичные математические модели, соответствующие физическим явлениям, формулируются в виде уравнений математической физики. Большинство реальных процессов описывается нелинейными уравнениями и лишь в первом приближении (при малых значениях параметров, малых отклонениях от равновесия и др.) эти уравнения можно заменить линейными. После того как задача сформулирована в математической форме, необходимо найти ее решение. Но что значит решить математическую задачу? Только в исключительных случаях удается найти решение в явном виде, например в виде ряда. Иногда утверждение «задача решена» означает, что доказано существование и единственность решения. Ясно, что этого недостаточно для практических приложений. Необходимо еще изучить качественное поведение решения и найти те или иные количественные характеристики. Именно на этом этапе требуется привлечение ЭВМ и, как следствие, развитие численных методов (см. III на рис. 1). Под численным методом здесь понимается такая интерпретация математической модели («дискретная модель»), которая доступна для реализации на ЭВМ. Например, если математическая модель представляет собой дифференциальное уравнение, то численным методом может быть аппроксимирующее его разностное уравнение совместно с алгоритмом, позволяющим отыскать решение этого разностного уравнения. Результатом реализации численного метода на ЭВМ является число или таблица чисел. Отметим, что в настоящее время помимо собственно численных методов имеются также методы, которые позволяют проводить на ЭВМ аналитические выкладки. Однако аналитические методы для ЭВМ не получили пока достаточно широкого распространения. Чтобы реализовать численный метод, необходимо составить программу для ЭВМ (см. IV на рис. 1) или воспользоваться готовой программой. После отладки программы наступает этап проведения вычислений и анализа результатов (V). Полученные результаты изучаются с точки зрения их соответствия исследуемому явлению и, при необходимости, вносятся исправления в численный метод и уточняется математическая модель. Такова в общих чертах схема вычислительного экспе­римента. Его основу составляет триада: модель — метод (алгоритм) — программа. Опыт решения крупных задач показывает, что метод математического моделирования и вычислительный эксперимент соединяют в себе преимущества традиционных теоретических и экспериментальных методов исследования. Можно указать такие крупные области применения вычислительного эксперимента, как энергетика, аэрокосмическая техника, обработка данных натурного эксперимента, совершенствование технологических процессов. 2. Вычислительный алгоритм. Предметом данной книги является изложение вопросов, отражающих этапы III, IV, V вычислительного эксперимента. Таким образом, здесь не обсуждаются исходные задачи и их математиче­ская постановка. Необходимо подчеркнуть, что процесс исследования исходного объекта методом математического моделирования и вычислительного эксперимента неизбежно носит приближенный характер, потому что на каждом этапе вносятся те или иные погрешности. Так, построение математической модели связано с упрощением исходного явления, недостаточно точным заданием коэффициентов уравнения и других входных данных. По отношению к численному методу, реализующему данную математическую модель, указанные погрешности являются неустранимыми, поскольку они неизбежны в рамках данной модели. При переходе от математической модели к численному методу возникают погрешности, называемые погрешностями метода. Они связаны с тем, что всякий численный метод воспроизводит исходную математическую модель приближенно. Наиболее типичными погрешностями метода являются погрешность дискретизации и погрешность округления. Поясним причины возникновения таких погреш­ностей. Обычно построение численного метода для заданной математической модели разбивается на два этапа: а) фор­мулирование дискретной задачи, б) разработка вычисли­тельного алгоритма, позволяющего отыскать решение дискретной задачи. Например, если исходная математическая задача сформулирована в виде системы дифференциальных уравнений, то для численного решения необходимо заменить ее системой конечного, может быть, очень большого числа линейных или разностных алгебраических уравнений. В этом случае говорят, что проведена дискретизация исходной математической задачи. Простейшим примером дискретизации является построение разностной схемы, путем замены дифференциальных выражений конечно-разностными отношениями. В общем случае дискретную модель можно рассматривать как конечномерный аналог исходной математической задачи. Ясно, что решение дискретизированной задачи отличается от решения исходной задачи. Разность соответствующих решений и называется погрешностью дискретизации. Обычно дискретная модель зависит от некоторого параметра (или множества параметров) дискретизации, при стремлении которого к нулю должна стремиться к нулю и погрешность дискретизации. При этом число алгебраических уравнений, составляющих дискретную модель, неограниченно возрастает. В случае разностных методов таким параметром является шаг сетки. Как уже отмечалось, дискретная модель представляет собой систему большого числа алгебраических уравнений. Невозможно найти решение такой системы точно и в явном виде. Поэтому приходится использовать тот или иной численный алгоритм решения системы алгебраических уравнений. Входные данные этой системы, а именно коэффициенты и правые части, задаются в ЭВМ не точно, а с округлением. В процессе работы алгоритма погрешности округления обычно накапливаются, и в результате решение, полученное на ЭВМ, будет отличаться от точного решения дискретизированной задачи. Результирующая погрешность называется погрешностью округления (иногда ее называют вычислительной погрешностью). Величина этой погрешности определяется двумя факторами: точностью представления вещественных чисел в ЭВМ и чувствительностью данного алгоритма к погрешностям округления. Алгоритм называется устойчивым, если в процессе его работы вычислительные погрешности возрастают незначительно, и неустойчивым — в противоположном случае. При использовании неустойчивых вычислительных алгоритмов накопление погрешностей округления приводит в процессе счета к переполнению арифметического устройства ЭВМ. Итак, следует различать погрешности модели, метода и вычислительную. Какая же из этих трех погрешностей является преобладающей? Ответ здесь неоднозначен. Ви­димо, типичной является ситуация, возникающая при ре­шении задач математической физики, когда погрешность модели значительно превышает погрешность метода, а погрешностью округления в случае устойчивых алгоритмов можно пренебречь по сравнению с погрешностью ме­тода. С другой стороны, при решении, например, систем обыкновенных дифференциальных уравнений возможно применение столь точных методов, что их погрешность будет сравнима с погрешностью округления. В общем случае нужно стремиться, чтобы все указанные погрешности имели один и тот же порядок. Например, нецелесообразно пользоваться разностными схемами, имеющими точность , если коэффициенты исходных уравнений задаются с точностью . 3. Требования к вычислительным методам. Одной и той же математической задаче можно поставить в соответствие множество различных дискретных моделей. Однако далеко не все из них пригодны для практической реализации. Вычислительные алгоритмы, предназначенные для быстродействующих ЭВМ, должны удовлетворять много­образным и зачастую противоречивым требованиям. Попытаемся здесь сформулировать основные из этих требований в общих чертах. Можно выделить две группы требований к численным методам. Первая группа связана с адекватностью дискретной модели исходной математической задаче, и вторая группа – с реализуемостью численного метода на ЭВМ. К первой группе относятся такие требования, как сходимость численного метода, выполнение дискретных аналогов законов сохранения, качественно правильное поведение решения дискретной задачи. Поясним эти требования. Предположим, что дискретная модель математической задачи представляет собой систему большого, но конечного числа алгебраических уравнений. Обычно, чем точнее мы хотим получить решение, тем больше уравнений приходится брать. Говорят, что численный метод сходится, если при неограниченном увеличении числа уравнений решение дискретной задачи стремится к решению исходной задачи. Поскольку реальная ЭВМ может оперировать лишь с конечным числом уравнений, на практике сходимость, как правило, не достигается. Поэтому важно уметь оценивать погрешность метода в зависимости от числа уравнений, составляющих дискретную модель. По этой же причине стараются строить дискретную модель таким образом, чтобы она правильно отражала качественное поведение решения исходной задачи даже при сравнительно небольшом числе уравнений. Например, дискретной моделью задачи математической физики может быть разностная схема. Для ее построения область изменения независимых переменных заменяется дискретным множеством точек – сеткой, а входящие в исходное уравнение производные за­меняются, на сетке, конечно-разностными отношениями. В результате получаем систему алгебраических уравнений относительно значений искомой функции в точках сетки. Число уравнений этой системы равно числу точек сетки. Известно, что дифференциальные уравнения математической физики являются следствиями интегральных законов сохранения. Поэтому естественно требовать, чтобы для разностной схемы выполнялись аналоги таких законов сохранения. Разностные схемы, удовлетворяющие этому требованию, называются консервативными. Оказалось, что при одном и том же числе точек сетки консервативные разностные схемы более правильно отражают поведение решения исходной задачи, чем неконсервативные схемы. Сходимость численного метода тесно связана с его корректностью. Предположим, что исходная математическая задача поставлена корректно, т.е. ее решение существует, единственно и непрерывно зависит от входных данных. Тогда дискретная модель этой задачи должна быть построена таким образом, чтобы свойство корректности сохранилось. Таким образом, в понятие кор­ректности численного метода включаются свойства однозначной разрешимости соответствующей системы уравнений и ее устойчивости по входным данным. Под устойчивостью понимается непрерывная зависимость решения от входных данных, равномерная относительно числа уравнений, составляющих дискретную модель. Вторая группа требований, предъявляемых к численным методам, связана с возможностью реализации данной дискретной модели на данной ЭВМ, т. е. с возможностью получить на ЭВМ решение соответствующей системы алгебраических уравнений за приемлемое время. Основным препятствием для реализации корректно поставленного алгоритма является ограниченный объем оперативной памяти ЭВМ и ограниченные ресурсы времени счета. Реальные вычислительные алгоритмы должны учитывать эти обстоятельства, т. е. они должны быть экономичными как по числу арифметических действий, так и по требуемому объему памяти. Численные методы алгебры и анализа 1 Решение систем линейных алгебраических уравнений Рассмотрим систему линейных алгебраических урав­нений: (1.1) или в матричной форме: Aх=b, (1.2) где: A={aij} квадратная матрица размерности (mm,); х=(х1,…., хm)T; T – операция транспонирования; b=(b1,….,bm)T; detA0. Предположим, что определитель матрицы A не равен нулю. Тогда решение х существует и единственно. На практике встречаются системы, имеющие большой поря­док. Методы решения системы (1.1) делятся на две группы: 1) прямые (точные методы); 2) итерационные методы (приближенные). 1.1 Точные методы В точных методах решение х находится за конечное число действий, но из-за погрешности округления и их накопления прямые методы можно назвать точными, только отвлекаясь от погрешностей округления. 1.1.1 Метод Гаусса Вычисления с помощью метода Гаусса (который называют также методом последовательного исключения неизвестных) состоят из двух основных этапов: прямого хода и обратного хода. Прямой ход метода заключается в последовательном исключении неизвестных из системы для преобразования ее к эквивалентной системе с треугольной матрицей. На этапе обратного хода производят вычисления значений неизвестных. Рассмотрим простейший вариант метода Гаусса, называемый схемой единственного деления. Прямой ход метода 1-й шаг. Предположим, что а110. Поделим первое уравнение на этот элемент, который назовем ведущим элементом первого шага : . (1.3) Остальные уравнения системы (1.1) запишем в виде , (1.4) где i=. Уравнение (1.3) умножаем на ai1 и вычитаем из i-го уравнения системы (1.4). Это позволит обратить в нуль коэффициенты при x1 во всех уравнениях, кроме первого. Получим эквивалентную систему вида: . (1.5) , где i,j=. Система (1.5) имеет матрицу вида: . Дальше работаем с укороченной системой, т.к. х1 входит только в 1-ое уравнение . 2-й шаг. На этом шаге исключаем неизвестное х2 из уравнений с номерами i=3,4,…,m. Если ведущий элемент второго шага , то из укороченной системы ана­логично исключаем неизвестное x2 и получаем матрицу коэффициентов такого вида: . Аналогично повторяем указанные действия для неиз­вестных х3,х4,...,хm-1 и приходим к системе : . (1.6) Эта система с верхней треугольной матрицей: . Обратный ход метода. Из последнего уравнения сис­темы (1.6) находим хm, из предпоследнего хm-1, ..., из пер­вого уравнения – х1. Общая формула для вычислений: xm=ym/cmm, , (i=m-1,…,1). Для реализации метода Гаусса требуется примерно (1/3)m3 арифметических операций, причем большинство из них приходится на прямой ход. Ограничение метода единственного деления заключа­ется в том, что ведущие элементы на k-ом шаге ис­ключения не равны нулю, т.е. ≠0. Но если ведущий элемент близок к нулю, то в процессе вычисления может накапливаться погрешность. В этом случае на каждом шаге исключают не хk, a хj (при jk). Такой подход называется методом выбора главного элемента. Для этого выбирают неизвестные xj с наибольшим по абсолютной величине коэффициентом либо в строке, либо в столбце, либо во всей матрице. Для его реализации требуется − арифметических действий. 1.1.2 Связь метода Гаусса с разложением матрицы на множители. Теорема об LU разложении. Пусть дана система Aх=b (1.1), которая при прямом ходе преобразуется в эквивалентную систему (1.6) и запишем ее в виде Cх=y, (1.6*) где С – верхняя треугольная матрица с единицами на главной диагонали, полученная из (1.6) делением последнего уравнения системы на сmm. Как связаны в системе (1.1) элементы b и элементы y из (1.6*)? Если внимательно посмотреть на прямой ход метода Гаусса, то можно увидеть, что . Для произвольного j имеем , (1.7) где j=, dji – числовые коэффициенты: . (1.8) Можно записать систему: b=Dy, где D−нижняя треугольная матрица с элементами на главной диагонали (j=,). В связи с тем, что в методе Гаусса угловые коэффици­енты не равны нулю , то на главной диагонали матрицы D стоят не нулевые элементы. Следовательно, эта матрица имеет обратную, тогда y=D-1b, Сx= D-1b. Тогда DCx=b. (1.9) В результате использования метода Гаусса, получили разложение матрицы А на произведение двух матриц A = DC, где D − нижняя треугольная матрица, у которой элементы на главной диагонали не равны нулю, а C – верхняя тре­угольная матрица с единичной диагональю. Таким образом, если задана матрица A и вектор b, то в методе Гаусса сначала производится разложение этой матрицы А на произведение D и C, а затем последова­тельно решаются две системы: Dy=b, Cx=y. (1.10) Из последней системы находят искомый вектор x. При этом разложение матрицы А на произведение СD – есть прямой ход метода Гаусса, а решение систем (1.10) обратный ход. Обозначим нижнюю треугольную матрицу через L, верхнюю треугольную матрицу − U. Теорема об LU разложении Введем обозначения: − угловой минор порядка j матрицы А, т.е. Теорема. Пусть все угловые миноры матрицы А не равны нулю (Δj0 для j=). Тогда матрицу А можно представить единственным образом в виде произведе­ния А=L*U. Идея доказательства. Рассмотрим матрицу А второго порядка и будем искать разложение этой матрицы в виде L и U. . Сопоставляя эти два равенства, определяем элементы матриц L и U (перемножим и приравняем неизвестные). Система имеет единственное решение. Методом математической индукции сказанное можно обобщить для матрицы размерности mm. Следствие. Метод Гаусса (схему единственного деле­ния) можно применять только в том случае, когда угловые миноры матрицы А не равны нулю. 1.1.3 Метод Гаусса с выбором главного элемента Может оказаться так, что система (1.1) имеет единст­венное решение, хотя какой либо из миноров матрицы А равен нулю. Заранее неизвестно, что все угловые миноры матрицы А не равны нулю. В этом случае можно использовать метод Гаусса с выбором главного элемента. 1. Выбор главного элемента по столбцу, когда на k-ом шаге исключения в качестве главного элемента выбирают максимальный по модулю коэффициент при неизвестном хk в уравнениях с номерами i=k,k+1,…,m. Затем уравнение, соответствующее выбранному коэффициенту, меняют местами с k-ым уравнением системы, чтобы ведущий элемент занял место коэффициента .После перестановки исключение неизвестного хk выполняют как в схеме единственного деления. ПРИМЕР. Пусть дана система второго порядка Предположим, чтоа21>а11, тогда переставим уравнения и применяем первый шаг прямого хода метода Гаусса. В этом случае имеет место перенумерация строк. 2. Выбор главного элемента по строке, т.е. производится перенумерация неизвестных системы. При а12> а11, на первом шаге вместо неизвестного х1 исключают х2: . К этой системе применяем первый шаг прямого хода метода Гаусса. 3. Поиск главного элемента по всей матрице заклю­чается в совместном применении методов 1 и 2. Всё это приводит к уменьшению вычислительной погрешности, но может замедлить процесс решения задачи. 1.1.4 Метод Холецкого (метод квадратных корней) Пусть дана система Ах=b, (1.11) где А − симметричная положительно определенная матрица. Тогда решение системы (1.11) проводится в два этапа: 1. Симметричная матрица А представляется как произведение двух матриц А = L∙ LТ. Рассмотрим метод квадратных корней на примере сис­темы 4-го порядка: . Перемножаем матрицы в правой части разложения и сравниваем с элементами в левой части: , , , , , , , . 2. Решаем последовательно две системы Ly=b, LTx=у. Замечания 1) Под квадратным корнем может получиться отрица­тельное число, следовательно в программе необходимо предусмотреть использование правил действия с комплексными числами. 2) Возможно переполнение, если угловые элементы близки к нулю. 1.2 Итерационные методы решений систем алгеб­раических уравнений Итерационные методы обычно применяются для реше­ния систем большой размерности и они требуют приведе­ния исходной системы к специальному виду. Суть итерационных методов заключается в том, что решение х системы (1.1) находится как предел последова­тельности . Так как за конечное число итераций предел не может быть достигнут, то задаётся малое число  − точность, и последовательные приближения вычисляют до тех пор, пока не будет выполнено неравенство , где n=n() – функция , ||x|| − норма вектора. Определения основных норм в пространстве векторов и матриц. Для вектора x=( x1,x2,…,xn)T нормы вычисляются по следующим формулам: ; ; . Согласованные с ними нормы в пространстве матриц: ; ; − величина, называемая евклидовой нормой матрицы A. Прямые методы рассчитаны для решения систем, порядок которых не больше 100, иначе для практических вычислений используются итерационные методы. 1.2.1 Метод Якоби (простых итераций) Исходную систему (1.11) Ах=b преобразуем к виду: (1.12) где i=1,2,...,m; aii0. Первая сумма равна нулю, если верхний предел сум­мирования меньше нижнего. Так (1.12) при i=1 имеет вид По методу Якоби (n+1 приближение хi) ищем по формуле (1.13) где n – номер итерации (0,1,…); i=. Итерационный процесс (1.13) начинается с начальных значений , которые в общем случае задаются произ­вольно, но предпочтительнее за взять свободные члены исходной системы. Условие окончания счета: , где i=. 1.2.2 Метод Зейделя Система (1.11) преобразуется к виду (1.12) и органи­зуется итерационная процедура, где неизвестные хi на n+1 шаге определяются по формулам (1.14) Например, (1.15) (1.16) и так далее. Итерационные процессы (1.13) и (1.14) сходятся, если норма матрицы А (А − матрица коэффициентов при неиз­вестных в правой части систем (1.13) и (1.14)) удовлетво­ряет условию: ||A||<1. 1.2.3 Матричная запись методов Якоби и Зейделя Исходную матрицу системы (1.11) представим в виде суммы трёх матриц A=A1+D+A2, где D − диагональная матрица; D =diаg[а11а22…аmm]; A1 − нижняя треугольная матрица; A2 − верхняя треугольная матрица. Пример: Дана матрица размерности (33): . А1 А2 D Тогда исходную систему (1.11) можно записать в виде x=–D-1A1 x – D-1A2 x+D-1 b. Тогда метод Якоби можно записать в виде: или . (1.17) В матричной форме метод Зейделя будет выглядеть: или . (1.18) Преобразуем формулы (1.17) и (1.18): , (1.19) . (1.20) Из (1.19) и (1.20) видно, что если итерационный метод сходится, то он сходится к точному решению. Иногда при решении задач большой размерности, в итерационные методы вводятся числовые параметры, которые могут зависеть от номера итерации. Пример для метода Якоби. , где t – числовой параметр. Возникают вопросы: 1) При каких значениях t сходимость будет наиболее быстрой? 2) При каких значениях t метод сходится? На примере двух методов просматривается вывод о том, что одни и те же методы можно записывать несколькими способами. Поэтому вводят каноническую (стандартную) форму записи: . (1.21) Формула (1.21) получена путем объединения (1.19) и (1.20). Матрица Dn+1 здесь задает тот или иной метод. Если существует обратная матрица к этой матрице, то из по­следней системы мы можем найти все неизвестные. 1. Метод (1.21) – явный, если матрица Dn совпадает с единичной матрицей и неявный – в противном случае. 2. Метод (1.21) – стационарный, если матрица Dn+1=D, и параметр t не зависит от номера итерации и нестацио­нарный – в противном случае. 1.2.4 Метод Ричардсона Явный метод с переменным параметром t: , (1.21а) называется методом Ричардсона. 1.2.5 Метод верхней релаксации (обобщённый метод Зейделя) , (1.21б) где  – числовой параметр. Если матрица А − симметричная и положительно определена, то последний метод сходится при (0 <  < 2). Последнюю формулу запишем в следующем виде: , (1.22) где Е – единичная матрица. Тогда для вычисления неизвестных хi (i=) можно записать итерационную процедуру в виде: . (1.23) Например, для х1 это будет такое выражение: . 1.2.6 Сходимость итерационных методов Рассмотрим систему Ax=b, где А – невырожденная действительная матрица. Для решения системы рассмотрим одношаговый стационарный метод , (1.24) при n=0,1,2…. Предположим, что задан начальный вектор решения. Тогда метод (1.24) сходится, если норма вектора Теорема. Условие сходимости итерационного ме­тода. Пусть А – симметричная положительно определенная матрица и выполнено условие D - 0.5tA > 0 (где t > 0). Тогда итерационный метод (1.24) сходится. Следствие 1. Пусть А – симметричная и положительно определенная матрица с диагональным преобладанием, то есть , при j=1,2,…,m. Тогда метод Якоби сходится. Следствие 2. Пусть А – симметричная и положительно определенная матрица с диагональным преобладанием, тогда метод верхней релаксации сходится при (0< <2). Проверяется, при каком значении  метод достигает заданной точности быстрее. В частности, при  = 1 метод верхней релаксации превращается в метод Зейделя, следовательно, при  = 1 метод Зейделя сходится. Теорема. Итерационный метод (1.24) сходится при любом начальном векторе x0 тогда и только тогда, когда все собственные значения матрицы по модулю меньше единицы. 2 Плохо обусловленные системы линейных алгебраических уравнений Дана система линейных алгебраических уравнений Ах=b (2.1) Если система плохо обусловлена, то это значит, что погрешности коэффициентов матрицы А и правых частей b или же погрешности их округления сильно искажают решение системы. Для оценки обусловленности системы вводят число обусловленности MА . Чем больше значение MА,тем система хуже обусловлена. Свойства числа обусловленности: 1) МЕ =1; 2) MА ³ 1; 3) MА³, где lмах, lmin – соответственно максимальное и минимальное собственные числа матрицы А; 4) MАВ £ MА * MВ; 5) Число обусловленности матрицы А не меняется при умножении матрицы на произвольное число a¹0. Найдем выражение для полной оценки погрешности решения системы. Пусть в системе (2.1) возмущены коэффициенты матрицы А и правая часть b, т.е. , , . Теорема. Пусть матрица A имеет обратную матрицу, и выполняется условие . Тогда матрица имеет обратную и справедлива следующая оценка относительной погрешности: . х1 а) х1 б) х1 в) х2 х2 Рисунок 2 – а) система имеет единственное решение; б) система не имеет решения; в) система плохо обусловлена. В случае в) малейшее возмущение системы сильно меняет положение точки пересечения прямых. В качестве примера рассмотрим систему Решение этой системы x1 » 1.981 x2 » 0.4735. Оценим влияние погрешности правых частей на результат. Рассмотрим “возмущенную” систему с правой частью b* = (2.505 , 2.415) и решим эту систему: x1* » 2.877 x2* » -0.4629. Относительная погрешность правой части d (b) = 0.005/2.51 » 0.28% привела к относительной погрешности решения d (x*) =0.9364/1.981 » 47.3%. Погрешность возросла примерно в 237 раз. Число обусловленности системы (2.1) приблизительно равно 237. Подобные системы называются плохо обусловленными. Возникает вопрос: какими методами можно решать такие системы? 2.1 Метод регуляризации для решения плохо обусловленных систем Рассмотрим систему Ах=b. (2.1) Для краткости перепишем эту систему в эквивалентный форме . (2.2) Для примера рассмотрим систему . Тогда ее можно представить как (2х1-х2-1)2+(х1-2х2-2)2=0. (2.2*) Решение системы (2.2) совпадает с решением системы (2.2*). Если коэффициенты A или b известны неточно, то решение также является не точным, поэтому вместо равенства можем потребовать приближенного выполнения равенства и в этом виде задача становится не определенной и нужно добавить дополнительные условия. В качестве дополнительного условия вводят требование, чтобы решение как можно меньше отклонялось от заданного х0 т.е. (х-х0, х-х0) было минимальным. Следовательно, приходим к регуляризованной задаче вида (Ах-b, Ax-b)+(x-x0, x-x0)=min, (2.3) где  >0. Используя свойства скалярного произведения, выражение (2.3) перепишем в виде (x,ATAx)-2(x,ATb)+(b,b)+[(x,x)-2(x,x0)+(x0,x0)]=min. (2.4) Варьируя x в уравнении (2.4), получим уравнение вида (ATА+E)x=ATb+x0. (2.5) Система (2.5) – система линейных алгебраических уравнений, эквивалентная системе (2.1). Систему (2.5) решаем с помощью метода Гаусса или с помощью метода квадратных корней. Решая систему (2.5) найдем решение, которое зависит от числа . Выбор управляющего параметра . Если =0, то система (2.5) перейдет в плохо обусловленную систему (2.1). Если же – велико, то система (2.5) переходит в хорошо обусловленную систему и решение этой системы может сильно отличаться от решения системы (2.1). Оптимальное значение  – это такое число, при котором система (2.5) удовлетворительно обусловлена. На практике пользуются невязкой вида , и эту невязку сравнивают по норме с известной погрешностью правых частей и с влиянием погрешности коэффициентов матрицы . Если – слишком велико, то или . Если – мало, то или . Поэтому проводят серию расчетов, при различных  и в качестве оптимального значения выбирают то значение , когда выполнено следующее условие . Для выбора вектора х0 нужно знать приближенное решение или же, если приближенное решение трудно определить, то х0 =0. 2.2 Метод вращения (Гивенса) Метод Гивенса, как и метод Гаусса состоит из прямого и обратного ходов. Прямой ход метода. Исключаем неизвестное х1 из всех уравнений, кроме первого. Для исключения х1 из 2-го уравнения вычисляют числа , , где  и  такие, что , . Первое уравнение системы заменяем линейной комбинацией первого и второго уравнений с коэффициентами и , а второе уравнение такой же комбинацией с и –. В результате получим систему . (2.6) Здесь где j=. Преобразование системы (2.1) к системе (2.6) эквивалентно умножению слева матрицы A и вектора b на матрицу С12 вида . Аналогично для исключения из третьего уравнения вычисляем числа и , такие, что . Затем первое уравнение системы (2.6) заменяем линейной комбинацией первого и третьего уравнений с коэффициентами ,, а третье уравнение системы (2.6) заменяем линейной комбинацией тех же уравнений, но с коэффициентами и –. Это преобразование эквивалентно умножению слева на матрицу . Исключая неизвестное х1 из всех последующих уравнений получим систему А(1) х=b(1), где матрица на первом шаге A(1)=C1m…C13C12A, , а вектор правых частей . Здесь и далее через Сkj обозначена матрица элементарного преобразования, отличающаяся от единичной матрицы Е только четырьмя элементами. Действие матрицы Сkj на вектор x эквивалентно повороту вектора x вокруг оси, перпендикулярной плоскости OXkXj на угол φkj такой, что , . Операцию умножения на матрицу Сkj называют плоским вращением или преобразованием Гивенса. Первый этап состоит из m-1 шагов, в результате чего получается система (2.7) В матричной форме получаем А(1) х=b(1). На втором этапе, состоящем из m-2 шагов, из уравнений системы (2.7) с номерами 3,4,…,m исключают неизвестное х2. B результате получим систему В матричной форме получаем , где , . После завершения (m-1)-го шага придем к системе с верхней треугольной матрицей вида , где , . Обратный ход метода вращений проводится точно так же, как и для метода Гаусса. 3 Решение нелинейных уравнений и систем нелинейных уравнений Рассмотрим систему нелинейных уравнений с m неизвестными вида . (3.1) Задача решения такой системы является более сложной, чем нахождение корней одного нелинейного уравнения, и чем задача решения линейных алгебраических уравнений. В отличие от систем линейных уравнений здесь использование прямых методов исключено и решение находится с использованием итерационных методов, т.е. находится приближенное решение x* = (x1*, ... , xm*), удовлетворяющее при заданном  > 0 условию . Задача (3.1) совсем может не иметь решения или же число решений может быть произвольным. Введем векторную запись решения задачи: x=(x1,…,xm)T, f=(f1,…,fm)T, f(x)=0. (3.2) Будем считать, что функции fi непрерывно дифференцируемы в некоторой окрестности точки х. Введем матрицу Якоби . Как и в случае решения одного уравнения начинаем с этапа локализации решения (отделения корней). Пример. Дана система 2-х уравнений с двумя неизвестными . Найдем на плоскости место расположения решения. Строим графики уравнений этой системы: а) – график 1-го уравнения, б) – график 2-го уравнения, в) – совмещенные графики. а) б) в) Рисунок 3 – Графики уравнений системы Определяем границы координат пересечения графиков. Данная система имеет три решения. Координаты точек (B,C,A): B: x1=4, x2=4 C: 3,5 < x1 < 4; 1,5 < x2 < 2,5. Точки А и С симметричны относительно прямой х1=х2. Координаты точки С определим приближенно: x13,8, x22. Обусловленность и корректность решения системы (3.1). Предположим что система (3.1) имеет решение х и в некоторой окрестности этого решения матрица Якоби не вырождена. Это означает, что в указанной окрестности нет других решений системы. В одномерном случае нахождение корня нелинейного уравнения приводит к определению интервала неопределенности (х*–, х*+). Так как значения функции f(x) чаще всего вычисляются на ЭВМ с использованием приближенных методов нельзя ожидать, что в окрестности корня относительная погрешность окажется малой. Сама погрешность корня ведет себя крайне нерегулярно и в первом приближении может восприниматься как некоторая случайная величина. На рисунке 4а представлена идеальная ситуация, отвечающая исходной математической постановке задачи, а на рисунке 4б – реальная, соответствующая вычислениям значений функции f на ЭВМ. а) б) y y x* а x* b x х х*–d х*+d Рисунок 4 – Графическое изображение определения интервала неопределенности В этом случае мы не можем определить, какая же точка в интервале неопределённости является решением. Радиус интервала неопределенности d прямо пропорционален погрешности вычисления значения f . Кроме того, d возрастает (обусловленность задачи ухудшается) с уменьшением . Оценить величину d довольно сложно, но выполнить это необходимо по следующим причинам: • не имеет смысла ставить задачу о вычислении корня с точностью ε < d; • после попадания очередного приближения в интервал неопределенности или близко от него, вычисления следует прекратить (этот момент для итерационных методов определяется крайне нерегулярным поведением приближений). Если случай многомерный, то получаем некоторую область неопределённости D, и можем получить оценку радиуса  этой области: Роль абсолютного числа обусловленности играет норма матрицы, обратной матрице Якоби . Чем число обусловленности больше, тем хуже эта система обусловлена. 3.1 Метод простых итераций Систему (3.1) преобразуем к следующему эквивалентному виду: . (3.3) Или в векторной форме х=φ(х) (3.4) Пусть задано начальное приближение . Подставляем его в правую часть системы (3.4) и получаем x(1)= (x(0)), продолжая подстановку, находим х(2) и т.д. Получим последовательность точек , которая приближается к искомому решению х. 3.1.1 Условия сходимости метода. Пусть '(x) – матрица Якоби, соответствующая системе (3.4) и в некоторой -окрестности решения х функции (i=1,2,…,m) дифференцируемы и выполнено неравенство вида: , где (0  q < 1), q − постоянная. Тогда независимо от выбора х(0) из -окрестности корня итерационная последовательность {хk} не выходит за пределы данной окрестности, метод сходится со скоростью геометрической прогрессии и справедлива оценка погрешности . 3.1.2 Оценка погрешности. В данной окрестности решения системы, производные функции i(x) (i=1,…,m) должны быть достаточно малы по абсолютной величине. Таким образом, если неравенство не выполнено, то исходную систему (3.1) следует преобразовать к виду (3.3). Пример. Рассмотрим предыдущий пример и приведем систему к удобному для итераций виду Проверяем условие сходимости вблизи точки С. Вычислим матрицу Якоби . Так как x13,8, x22, то при этих значениях вычисляем норму матрицы ||||||(3,8 ; 2) || 0,815. Запишем итерационную процедуру Следовательно, метод простых итераций будет сходиться со скоростью геометрической прогрессии, знаменатель которой q0,815. Вычисления поместим в таблице 1. Таблица 1 Решение системы нелинейных уравнений k 1 … 8 9 3,80000 3,75155 …. 3,77440 x1=3,77418 2,00000 2,03895 … 2,07732 x2=2,07712 При k=9 критерий окончания счета выполняется при =10-3 и можно положить x1=3,7740,001 и x2=2,077 0,001. 3.2 Метод Ньютона Суть метода состоит в том, что система нелинейных уравнений сводится к решению систем линейных алгебраических уравнений. Пусть дана система (3.1) и задано начальное приближение x(0). Приближение к решению х строим в виде последовательности x(0), x(1), …, x(n). В исходной системе (3.1) каждую функцию где i=, раскладывают в ряд Тейлора в точке х(n) и заменяют линейной частью её разложения . В результате получим систему линейных алгебраических уравнений . (3.5) В матричной форме f(x(n))+ f'(x(n))*(x−x(n))=0, (3.6) где f ' − матрица Якоби. Предположим, что матрица невырожденная, то есть существует обратная матрица . Тогда система (3.6) имеет единственное решение, которое и принимается за очередное приближение x(n+1). Отсюда выражаем решение x(n+1) по итерационной формуле: . (3.7) Формула (3.7) и есть итерационная формула метода Ньютона для приближенного решения системы нелинейных уравнений. Замечание. В таком виде формула (3.7) используется редко в виду того, что на каждой итерации нужно находить обратную матрицу. Поэтому поступают следующим образом: вместо системы (3.6) решают эквиваленту ей систему линейных алгебраических уравнений вида f(x(n))*x(n+1) =−f(x(n)). (3.8) Это система линейных алгебраических уравнений относительно поправки x(n+1)= x(n+1) − x(n). Затем полагают x(n+1) =x(n) +x(n+1). (3.9) 3.2.1 Сходимость метода Теорема. Пусть в некоторой окрестности решения х системы (3.1) функции fi (при i=) дважды непрерывно дифференцируемы и матрица Якоби не вырождена. Тогда найдется такая малая  окрестность вокруг решения х, что при выборе начального приближения x0 из этой окрестности итерационный метод (3.7) не выйдет за пределы этой окрестности решения и справедлива оценка вида , где n − номер итерации. Метод Ньютона сходится с квадратичной скоростью. На практике используется следующий критерий остановки: . 4 Решение проблемы собственных значений Пусть дана квадратная матрица A размерностью (m*m) и существует такое число , что выполняется равенство A∙x=λ∙x, x≠0, тогда такое число  называется собственным значением матрицы А, а x– соответствующим ему собственным вектором. Перепишем это равенство в эквивалентной форме (A − E)x = 0 . (4.1) Система (4.1) – однородная система линейных алгебраических уравнений. Для существования нетривиального решения системы (4.1) должно выполняться условие det(A − E) = 0 . (4.2) Определитель в левой части уравнения является многочленом m-ой степени относительно , его называют − характеристическим определителем (характеристическим многочленом). Следовательно, уравнение (4.2) имеет m корней или m собственных значений. Среди них могут быть как действительные, так и комплексные корни. Задача вычисления собственных значений сводится к нахождению корней характеристического многочлена (4.2). Корни могут быть найдены одним из итерационных методов (в частности методом Ньютона). Если найдено некоторое собственное значение матрицы A, то, подставив это число в систему (4.1) и решив эту систему однородных уравнений, находим собственный вектор х, соответствующий данному собственному значению. Собственные вектора будем при нахождении нормировать (вектор х умножаем на ||х||-1, и таким образом они будут иметь единичную длину), нахождение собственных значений матрицы A и соответствующих им собственных векторов и есть полное решение проблемы собственных значений. А нахождение отдельных собственных значений и соответствующих им векторов − называется решением частной проблемы собственных значений. Эта проблема имеет самостоятельное значение на практике. Например, в электрических и механических системах собственные значения отвечают собственным частотам колебаний, а собственные вектора характеризуют соответствующие формы колебаний. Эта задача легко решается для некоторых видов матриц: диагональных, треугольных и трехдиагональных матриц. К примеру, определитель треугольной или диагональной матрицы равен произведению диагональных элементов, тогда и собственные числа равны диагональным элементам. Пример. Матрица А – диагональная . Тогда det(А−Е)=, а характеристическое уравнение имеет трехкратный корень =а. Собственными векторами для матрицы А будут единичные векторы Пример. Найдем собственные числа матрицы . Составим характеристический многочлен Используя метод Ньютона, определим один из корней уравнения Р3()=0, а именно -7,87279. Разделив многочлен P3(λ) на (−1) получим многочлен второй степени: P2(λ)=2+3,02711+2,66765. Решив квадратное уравнение, находим оставшиеся два корня: 2,3-1,513560,613841*i (комплексные сопряженные корни). Существуют прямые методы нахождения собственных значений и итерационные методы. Прямые методы неудобны для нахождения собственных значений для матриц высокого порядка. В таких случаях с учетом возможностей компьютера более удобны итерационные методы. 4.1 Прямые методы нахождения собственных значений 4.1.1 Метод Леверрье Метод разделяется на две стадии: - раскрытие характеристического уравнения, - нахождение корней многочлена. Пусть det(A-E) – есть характеристический многочлен матрицы А={aij} (i,j=1,2,…,m), т.е., и 1,2,…,m – есть полная совокупность корней этого многочлена (полный спектр собственных значений). Рассмотрим суммы вида (k=1,2,…,m), т.е. , (4.3) где − след матрицы. В этом случае при km справедливы формулы Ньютона для всех (1k m) , (4.4) Откуда получаем при k=1 р1 = -S1, при k=2 р2 = -1/2∙ (S2 + р1∙S1), . . . . . . . . . . . . . . при k=m рm = -1/n∙ (Sm + р1∙Sm-1 + р2∙Sm-2 +. . .+ + рm-1∙S1). (4.5) Следовательно, коэффициенты характеристического многочлена рi можно определить, если известны суммы S1,S2, ... ,Sm. Тогда схема алгоритма раскрытия характеристического определителя методом Леверрье будет следующей: 1) вычисляют степень матрицы: Аk=Аk-1∙А для k=1,…,m; 2) определяют Sk − суммы элементов, стоящих на главной диагонали матриц Аk; 3) по формулам (4.5) находят коэффициенты характеристического уравнения рi(i=1,2,…,m). 4.1.2 Усовершенствованный метод Фадеева Алгоритм метода: 1) вычисляют элементы матриц A1,A2 ,.., Am: (в конце подсчета Bm – нулевая матрица для контроля); 2) определяют коэффициенты характеристического уравнения рi: q1 = –р1, q2 = –р2,..., qm = –рm. Существуют и другие методы раскрытия характеристического определителя: метод Крылова, Данилевского и др. 4.1.3 Метод Данилевского Две матрицы A и B называются подобными, если одна получается из другой путем преобразования с помощью некоторой не вырожденной матрицы S: B=S-1∙AS, если это равенство справедливо, то матрицы A и B подобны, а само преобразование называется преобразованием подобия (переход к новому базису в пространстве m - мерных векторов). Пусть y − результат применения матрицы A к вектору х y=A∙х. Сделаем замену переменных: x=S∙x' , y=S∙y'. Тогда равенство y=A∙х преобразуется к виду y'=S-1∙A∙S∙x'. В этом случае матрица B и матрица A имеют одни и те же собственные числа. Это можно легко увидеть раскрыв определитель . Следовательно, матрицы A и B − подобные, имеют одни и те же собственные значения. Но собственные векторы х и не совпадают, они связаны между собой простым соотношением х = S∙х'. Такую матрицу A с помощью преобразования подобия или же последовательности таких преобразований можно привести к матрице Фробениуса вида: . Детерминант матрицы F – det(F) можно разложить по элементам первой строки: . Тогда коэффициенты характеристического многочлена матрицы А будут р1 = f11 , p2 = f12,…, pm = f1m.. Второй случай. Матрицу А преобразованием подобия можно привести к матрице В верхнего треугольного вида . Тогда собственными числами будут диагональные элементы матрицы B: . Третий случай. Матрицу A с помощью преобразования подобия можно привести к Жордановой форме , где i – собственные числа матрицы A; Si – константы (0 или 1); если Si=1, то i=i+1. К четвёртому случаю относятся матрицы, которые с помощью преобразования подобия можно привести к диагональному виду (матрица простой структуры): , у которой, как известно, собственными числами являются диагональные элементы. 4.1.4 Метод итераций определения первого собственного числа матрицы. Пусть дано характеристическое уравнение: det(A−∙E) = 0, где 1, 2,..., n − собственные значения матрицы А. Предположим, что |1|>|2||3|…|m|, т.е. 1 – наибольшее по модулю собственное число. Тогда для нахождения приближенного значения 1 используется следующая схема: 1) выбирают произвольно начальный вектор у(0); 2) строят последовательность итераций вида: 3) выбирают , тогда или , где и – соответствующие координаты векторов y(m) и y(m+1). Возникает вопрос выбора начального вектора у(0). При неудачном выборе можем не получить значения нужного корня, или же предела может не существовать. Этот факт при вычислении можно заметить по прыгающим значениям этого отношения, следовательно, нужно изменить у(0). В качестве первого собственного вектора можно взять вектор у(n+1) и пронормировать его. Пример. Найти наибольшее по модулю собственное значение и соответствующий ему собственный вектор матрицы А . 1) Выбираем начальный вектор . 2) Вычисляем последовательно векторы y(1), y(2), …, y(10). Вычисления помещаем в таблицу 2. Таблица 2 – Вычисление векторов y(n+1) y(0) А∙y(0) А2∙y(0) А3∙y(0) …….. А9∙y(0) А10∙y(0) 1 4 17 69 243569 941370 1 5 18 67 210663 812585 1 2 7 25 73845 284508 3) Вычисляем отношения координат векторов 4)Вычисляем l1 как среднее арифметическое . . 5) Определим соответствующий числу l1собственный вектор: Нормируем вектор y(10), разделив на его длину получим вектор Далее можем определить второе собственное число , где i=1,2,…,n. При вычислении собственных чисел подобным образом, будет накапливаться ошибка. Данная методика позволяет приближенно оценить собственные значения матрицы. 5 Задача приближения функции Постановка задачи. Пусть на отрезке [a,b] функция у=f(x) задана таблицей своих значений . Допустим, что вид функции f(x) неизвестен. На практике часто встречается задача вычисления значений функции у=f(x) в точках х, отличных от . Кроме того, в некоторых случаях, не смотря на то, что аналитическое выражение у=f(x) известно, оно может быть слишком громоздким и неудобным для математических преобразований (например, специальные функции). Кроме этого значения yi могут содержать ошибки эксперимента. Определение . Точки называются узлами интерполяции. Требуется найти аналитическое выражение функции F(x), совпадающей в узлах интерполяции со значениями данной функции, т.е. Определение. Процесс вычисления значений функции F(x) в точках отличных от узлов интерполирования называется интерполированием функции f(x). Если , то задача вычисления приближенного значения функции в т. х называется интерполированием, иначе – экстраполированием. Геометрически задача интерполирования функции одной переменной означает построение кривой, проходящей через заданные точки (рисунок 5). То есть задача в такой постановке может иметь бесконечное число решений. у x x x x х Рисунок 5 − Геометрическая иллюстрация задачи интерполирования функции Задача становится однозначной, если в качестве F(x) выбрать многочлен степени не выше n, такой что: Fn (x0)=y0, Fn (x1)=y1, ..., Fn (xn)=yn. Определение. Многочлен Fn(x), отвечающий вышеназванным условиям, называется интерполяционным многочленом. Знание свойств функции f позволяет осознанно выбирать класс G аппроксимирующих функций. Широко используется класс функций вида (5.1) являющихся линейными комбинациями некоторых базисных функций 0(x), ..., m(x). Будем искать приближающую функцию в виде многочлена степени m, с коэффициентами с0, ..., сm, которые находятся в зависимости от вида приближения. Функцию Фm(х) называют обобщенным многочленом по системе функций 0(х), 1(х), …,m(х), а число m – его степенью. Назовем обобщенный многочлен Фm(х) интерполяционным, если он удовлетворяет условию Фm(хi)=yi, (i=0,1,…,n). (5.2) Покажем, что условие (5.2) позволяет найти приближающую функцию единственным образом (5.3) Система (5.3) есть система линейных алгебраических уравнений относительно коэффициентов с0,с1,…,сm. Эта система n линейных уравнений имеет единственное решение, если выполняется условие m=n и определитель квадратной матрицы Р Определение. Система функций (x),...,n(x) называется Чебышевской системой функций на [a,b],если определитель матрицы отличен от нуля detP≠0 при любом расположении узлов xi[a,b], i=0,1,…n, когда среди этих узлов нет совпадающих. Если мы имеем такую систему функций, то можно утверждать, что существует единственный для данной системы функций интерполяционный многочлен Фm(х), коэффициенты которого определяются единственным образом из системы (5.3). Пример. При mn система функций 1, х, х2,…, хm линейно независима в точках х0, х1, …, хn, если они попарно различны. 5.1 Интерполяционный многочлен Лагранжа Рассмотрим случай, когда узлы интерполирования не равноотстоят друг от друга на отрезке [a,b]. Тогда шаг h=xi+1 −x≠const. Задача имеет единственное решение, если в качестве интерполирующей функции F(x) взять алгебраический многочлен Ln(x)=a0+a1 x+a2 x2+…+anxn, где анеизвестные постоянные коэффициенты. Используя условие (5.2) можем записать (5.4) Запишем это в виде: (5.5) Эта система однозначно разрешима, так как система функций 1,х,х2,…,хn линейно независима в точках х0,х1,…,хn. Однозначная разрешимость следует из того факта, что определитель этой системы (определитель Вандермонда) Без вывода приведем одну из форм записи интерполяционного многочлена Лагранжа (5.6) Определение. Этот многочлен называется интерполяционным многочленом Лагранжа и сокращенно записывается в виде Ln(x)= (5.7) На практике часто пользуются линейной и квадратичной интерполяцией. В этом случае формула Лагранжа имеет вид L1(x)= − при линейной интерполяции; L2(x)= − при квадратичной интерполяции. Рассмотрим теперь случай с равноотстоящими узлами. Тогда интерполяционная формула Лагранжа заметно упрощается. В этом случае шаг h=xi+1-xi=const. Введем в рассмотрение многочлен вида Введем обозначение q= , отсюда следует, что = , = , . . . . . . . . . . . = = Тогда многочлен Qi примет вид Произведя простейшие преобразования, получим выражение вида: Qi(q)== , где – число сочетаний из n элементов по i Тогда интерполяционный многочлен Лагранжа для равноотстоящих узлов имеет вид: . 5.1.1 Оценка погрешности интерполяционного многочлена Оценить погрешность интерполяционной формулы Лагранжа можно только тогда, когда известно аналитическое выражение интерполируемой функции, а точнее, если известно максимальное значение (n+1)-ой производной функции f(x) на отрезке [a,b]. Пусть |R(x)| =| f(x) −L(x)|, где Rn(x) –погрешность; f(x) − точное значение функции в точке х; Ln(x) − приближенное значение, полученное по полиному Лагранжа. Если обозначить через M=, где ξ[a,b], причем х0=а, хn= b, то . 5.2 Интерполяционные полиномы Ньютона 5.2.1 Интерполяционный многочлен Ньютона для равноотстоящих узлов Вычисление значений функции для значений аргумента, лежащих в начале таблицы удобно проводить, пользуясь первой интерполяционной формулой Ньютона. Для этого введем понятие конечной разности. Определение. Конечной разностью перового порядка называется разность между значениями функции в соседних узлах интерполяции. Тогда конечные разности в точках х0,х1,…,хn-1 , , . . . . . . . . . . . . Конечная разность второго порядка имеет вид: Рассмотрим некоторые свойства конечных разностей. Вторая конечная разность в точке хi . Аналогично третья конечная разность . Общее выражение для конечной разности n-го порядка имеет вид а вообще, конечная разность порядка m от конечной разности порядка n . Конечные разности n-го порядка от многочлена степени n – есть величины постоянные, а конечные разности n+1-го порядка равны нулю. Для вычисления значений функции в начале таблицы требуется построить интерполяционный многочлен степени n такой, что выполнены условия интерполяции . В силу единственности многочлена степени n, построенного по n+1 значениям функции f(x) многочлен , в конечном счете, совпадает с многочленом Лагранжа. Найдем этот многочлен в виде: , где аi(i=0,1,…, n) – неизвестные коэффициенты. Для нахождения а0 положим . Тогда , отсюда а0=у0. Для вычисления рассмотрим первую конечную разность для многочлена Pn(x) в точке х. В результате преобразований получим Вычислим первую конечную разность многочлена Pn(x) в точке х0 , но откуда . Чтобы определить коэффициент а2, составим конечную разность второго порядка Отсюда после преобразования получим . Вычисляя конечные разности более высоких порядков и полагая х=х0, придем к общей формуле для определения коэффициентов: (i=0,1,2,…,n). Подставим значения ai в многочлен, в результате получим первую интерполяционную формулу Ньютона: Первую интерполяционную формулу можно записать в том виде, в котором ее удобнее использовать для интерполирования в начале таблицы. Для этого введем переменную q=(x-x0)/h, где h– шаг интерполирования. Тогда первая формула примет вид 5.2.2 Вторая интерполяционная формула Ньютона Эта формула используется для интерполирования в конце таблицы. Построим интерполяционный многочлен вида Неизвестные коэффициенты а0,а1,…,аn подберем так, чтобы были выполнены равенства . Для этого необходимо и достаточно, чтобы (i=0,1,…,n). В случае, если положить x=xn, то сразу определяется коэффициент а0 . Из выражения для первой конечной разности найдем: Отсюда, полагая х=хn-1, получим . Из выражения для второй конечной разности найдем а2: . Общая формула для коэффициента аi имеет вид . Подставим эти коэффициенты в формулу многочлена и получим вторую интерполяционную формулу Ньютона: На практике используют формулу Ньютона в другом виде. Положим q=(x-xn)/h. Тогда . 5.3 Интерполирование сплайнами Многочлен Лагранжа или Ньютона на всем отрезке [a,b] с использованием большого числа узлов интерполирования часто приводит к плохому приближению, что объясняется накоплением погрешностей в ходе вычислений. Кроме того, из-за расходимости процесса интерполирования увеличение числа узлов не обязательно приводит к повышению точности вычислений. Поэтому построим такой вид приближения, который: • позволяет получить функцию, совпадающую с табличной функцией в узлах; • приближающая функция в узлах таблицы имеет непрерывную производную до нужного порядка; В силу вышесказанного на практике весь отрезок [a,b] разбивается на частичные интервалы и на каждом из них приближающая функция f(x) заменяется многочленом невысокой степени. Такая интерполяция называется кусочно-полиномиальной интерполяцией. Определение. Сплайн − функцией называют кусочно-полиномиальную функцию, определенную на отрезке [a,b] и имеющую на этом отрезке некоторое число непрерывных производных. Слово сплайн означает гибкую линейку, которую используют для проведения гладких кривых через определенное число точек на плоскости. Преимущество сплайнов – сходимость и устойчивость процесса вычисления. Рассмотрим частный случай (часто используемый на практике), когда сплайн определяется многочленом третьей степени. 5.3.1 Построение кубического сплайна Пусть на отрезке [a,b] в узлах сетки заданы значения некоторой функции f(x), т.е. , (i= 0,1,…, n). Сплайном, соответствующим этим узлам функции называется функция S(х), которая: 1) на каждом частичном отрезке является многочленом третьей степени; 2) функция S(x) и ее две первые производные непрерывны на [a,b]; 3) . На каждом частичном отрезке будем искать сплайн , где многочлен третьей степени . (5.8) То есть для нужно построить такую функцию, где подлежат определению. Для всего отрезка интерполирования [a,b], таким образом, необходимо определить 4n неизвестных коэффициента. Доопределим . Требование непрерывности функции S(x) приводит к условия (i=0,1,…,n-1). Отсюда из (5.8) получаем следующие уравнения: (i= 1,2,…,n-1). Введем шаг интерполирования . Тогда последнее равенство можно переписать в виде (i= 1,2,…,n). Из непрерывности первой производной следует (i=2,3,…,n), а из непрерывности второй производной hidi=ci−ci−1 (i=2,3,…,n). Объединив все три вида уравнений, получим систему из 3n-2 уравнений относительно 3n неизвестных . Два недостающих уравнения получим, задав граничные условия для функции S(x). Для этого воспользуемся граничными условиями для сплайн-функции в виде (концы гибкой линейки свободны). Тогда получим систему уравнений (5.9) Решая систему методом подстановки (исключаем из (5.9) неизвестные bi,di), получим систему: (5.10) Система (5.10) имеет трехдиагональную матрицу. Эта система может быть решена методом прогонки или Гаусса. Метод прогонки рассматривается в пункте 6.7.2. После решения системы коэффициенты сплайна определим через коэффициенты сi с помощью явных формул , (i= 1,2,…,n). Существуют специальные виды записи сплайнов на каждом из промежутков [xi ,xi+1] [9], которые позволяют уменьшить число неизвестных коэффициентов сплайна. Вводятся обозначения На отрезке [xi ,xi+1] кубический сплайн записывается в виде Кубический сплайн, записанный в таком виде, на каждом из промежутков [xi ,xi+1] непрерывен вместе со своей первой производной на [a,b]. Выберем mi таким образом, чтобы и вторая производная была непрерывна во всех внутренних узлах. Отсюда получим систему уравнений: где К этим уравнениям добавим уравнения, полученные из граничных условий В результате получаем аналогичную систему с трехдиагональной матрицей. Решаем систему линейных уравнений относительно коэффициентов mi методом пргонки. 5.3.2 Сходимость процесса интерполирования кубическими сплайнами Доказывается, что при неограниченном увеличении числа узлов на одном и том же отрезке [a,b] . Оценка погрешности интерполяции зависит от выбора сетки и степени гладкости функции f(x). При равномерной сетке (i=0,1,…,n) где. Другие постановки задачи интерполирования функций. 1. Если функция периодическая, то используется тригонометрическая интерполяция с периодом l, которая строится с помощью тригонометрического многочлена , коэффициенты которого находятся из системы уравнений (i= 1,2,…, 2n+1). 2. Выделяют приближение функций рациональными, дробно – рациональными и другими функциями. В данном пособии эти вопросы не рассматриваются. 5.4 Аппроксимация функций методом наименьших квадратов К такой задаче приходят при статистической обработке экспериментальных данных с помощью регрессионного анализа. Пусть в результате исследования некоторой величины x значениям поставлены в соответствие значения некоторой величины у. Требуется подобрать вид аппроксимирующей зависимости y=f(x), связывающей переменные х и у. Здесь могут иметь место следующие случаи. Во-первых: значения функции f(x) могут быть заданы в достаточно большом количестве узлов; во-вторых: значения таблично заданной функции отягощены погрешностями. Тогда проводить приближения функции с помощью интерполяционного многочлена нецелесообразно, т.к. - число узлов велико и пришлось бы строить несколько интерполяционных многочленов; - построив интерполяционные многочлены, мы повторили бы те же самые ошибки, которые присущи таблице. Будем искать приближающую функцию из следующих соображений: 1) приближающая функция не проходит через узлы таблицы и не повторяет ошибки табличной функции; 2) чтобы сумма квадратов отклонений приближающей функции от таблично заданной в узлах таблицы была минимальной. у уn yn-1 отклонения y1 y0 х0 х1 … хn-1 хn х Рисунок 6 – Графическое изображение отклонений приближающей функции от таблично заданной Рассмотрим линейную задачу наименьших квадратов. Определение. Уровень погрешности, допускаемый при снятии характеристики измеряемой величины, называется шумом таблицы. Пусть функция y=f(x) задана таблицей приближенных значений yif(xi), i=0,1,…,n, полученных с ошибками , где . Пусть даны функции , назовем их базисными функциями. Будем искать приближающую (аппроксимирующую) функцию в виде линейной комбинации базисных функций . (5.11) Такая аппроксимация называется линейной, а Фm(х) – обобщенным многочленом. Будем определять коэффициенты обобщенного многочлена c0,…,cm используя критерий метода наименьших квадратов. Согласно этому критерию вычислим сумму квадратов отклонений таблично заданной функции от искомого многочлена в узлах: . (5.12) Выражение для можно рассматривать как функцию от неизвестных c0,…, cm. Нас интересует, при каких значениях c0,…, cm, значение будет минимально. Для этого воспользуемся условием существования экстремума функции, а именно, найдем частные производные от по всем переменным c0,…, cm и приравняем их к нулю. Получим систему вида: . (5.13) Система (5.13) – система линейных уравнений относительно c0,…, cm. Чтобы систему (5.13) записать компактно, ведем определение. Определение. Скалярным произведением функций f на g на множестве точек называется выражение . Тогда систему (5.13) можно записать в виде: (5.13а) Системы (5.13) или (5.13а) будем называть нормальной системой уравнений. Решив ее, мы найдем коэффициенты c0,…,cm и следовательно, найдем вид аппроксимирующего многочлена. Напомним, что это возможно, если базисные функции линейно независимы, а все узлы различны. Осталось определить степень многочлена m. Прямому вычислению поддаются только значения среднеквадратичного отклонения , анализируя которые будем выбирать степень многочлена. Алгоритм выбора степени многочлена m. В случае, когда m=n мы получим интерполяционный многочлен, поэтому выберем m<0 и 2>0 должны быть такими, чтобы находилось между ними; 2) первоначально m выбирают произвольно, но учитывая условие, что m< 1, то степень необходимо уменьшить хотя бы на единицу; б) если <2, то степень необходимо увеличить хотя бы на единицу. 5) затем строим приближающую функцию . Очень часто для приближения по методу наименьших квадратов используются алгебраические многочлены степени mn, т.е. . Тогда нормальная система (5.13) принимает следующий вид: , (k= 0,1,…,m). (5.14) Запишем систему (5.14) в развернутом виде в двух наиболее простых случаях m =1 и m =2. В случае многочлена первой степени P1(x)=c0+c1x, нормальная система имеет вид (5.15) Для многочлена второй степени P2(x)=c0+c1x+c2x2, нормальная система имеет вид (5.16) 6 Численные методы решения задачи Коши для обыкновенных дифференциальных уравнений и систем дифференциальных уравнений Будем рассматривать задачу Коши для системы обыкновенных дифференциальных уравнений(ОДУ).Запишем систему в векторной форме , (6.1) где: u – искомая вектор-функция; t– независимая переменная; ; , m– порядок системы; координаты; t0; . Запишем систему (6.1) в развернутом виде , (6.2) где: i=1,...,m; . В случае i=1 – это будет ОДУ 1-го порядка, а при i=2 – система из двух уравнений первого порядка. В случае i=1 решение задачи Коши предполагает нахождение интегральной кривой, проходящей через заданную точку и удовлетворяющую заданному начальному условию. Задача состоит в том, чтобы найти искомую вектор-функцию u, удовлетворяющую (6.1) и заданным начальным условиям. Известны условия, гарантирующие существование и единственность решения (6.1) или (6.2). Предположим, что функции fi(i=1, …, m) непрерывны по всем аргументам в некоторой замкнутой области D={t}, где a,b– известные константы. Из непрерывности функций следует их ограниченность, т.е. функции сверху ограничены некоторой константой М: | fi |0. По оси t введем равномерную сетку с шагом  , т.е. рассмотрим систему точек . Обозначим через u(t) точное решение (6.4) , а через приближенные значения функций u в заданной системе точек. Приближенное решение является сеточной функцией, т.е. определено только в точках сетки . 6.1 Семейство одношаговых методов решения задачи Коши 6.1.1 Метод Эйлера (частный случай метода Рунге-Кутта) Уравнение (6.4) заменяется разностным уравнением , n=0,1,2,…,. В окончательной форме значения можно определить по явной формуле . (6.5) Вследствие систематического накопления ошибок метод используется редко или используется только для оценки вида интегральной кривой. Определение 1. Метод сходится к точному решению в некоторой точке t , если при , . Метод сходится на интервале (0,t], если он сходится в каждой точке этого интервала. Определение 2. Метод имеет р-й порядок точности, если существует такое число р>0, для которого при , где:  – шаг интегрирования; O – малая величина порядка τp. Так как , то метод Эйлера имеет первый порядок точности. Порядок точности разностного метода совпадает с порядком аппроксимации исходного дифференциального уравнения. 6.1.2 Методы Рунге-Кутта Метод Рунге-Кутта второго порядка точности Отличительная особенность методов Рунге-Кутта от метода (6.5) заключается в том, что значение правой части уравнения вычисляется не только в точках сетки, но и также в середине отрезков(промежуточных точках). Предположим, что приближенное значение yn решения задачи в точке t=tn уже известно. Для нахождения yn+1 поступают следующим образом: 1) используют схему Эйлера в таком виде (6.6) и отсюда вычисляют промежуточное значение yn+1/2; 2) воспользуемся разностным уравнением вида , (6.7) откуда найдем значение yn+1. Далее подставим значение в уравнение (6.7). Тогда получим разностное уравнение , (6.8) где . Можно показать, что метод (6.8) имеет второй порядок точности, т.е. . Реализация метода (6.8) в виде двух этапов (6.6),(6.7) называется методом предиктор–корректор (прогноза и коррекции) в том смысле, что на первом этапе (6.6) приближенное значение предсказывается с невысокой точностью , а на втором этапе (6.7) это предсказанное значение исправляется и результирующая погрешность имеет второй порядок . Будем рассматривать явные методы. Задаем числовые коэффициенты ai, bij, i=2,...,m; j=1,2,...,m-1 и σi, i=1,2,...,m. Последовательно вычисляем функции ; ; ; …………………………………………….. . Затем из формулы находим новое значение =y(tn+1). Здесь ,–числовые параметры, которые определяются или выбираются из соображений точности вычислений. Чтобы это уравнение аппроксимировало уравнение (6.4), необходимо потребовать выполнения условия . При m=1 получается метод Эйлера, при m=2 получаем семейство методов , (6.9) где: ; ; y0=u0. Семейство определяет явные методы Рунге-Кутта. Подставив нужные 1 и 2, получаем окончательную формулу. Точность этих методов совпадает с точностью аппроксимирующего метода и равна . Невязкой, или погрешностью аппроксимации метода (6.9) называется величина , полученная заменой в (6.9) приближенного решения точным решением. При 1+2=1 получим первый порядок точности. Если же потребовать дополнительно , то получим методы второго порядка точности вида при . Приведем один из методов Рунге-Кутта третьего порядка точности , где: ; ; . Метод Рунге-Кутта 4-го порядка точности , где: Методы Рунге-Кутта сходятся и порядок их точности совпадает с порядком аппроксимации разностным отношением. Теорема. Пусть правая часть уравнения (6.4) удовлетворяет условию Липшица по аргументу u с константой L, и пусть – невязка метода Рунге-Кутта. Тогда для погрешности метода при справедлива оценка , где: ; ; . На практике обычно пользуются правилом Рунге. Для этого сначала проводят вычисления с шагом τ , затем – τ/2. Если – решение при шаге τ, а – при шаге , p – порядок точности метода, то справедлива оценка. Тогда за оценку погрешности для метода четвертого порядка точности при шаге принимают величину . 6.2 Многошаговые разностные методы решения задачи Коши для обыкновенных дифференциальных уравнений Рассмотрим задачу Коши для обыкновенного дифференциального уравнения , (6.10) где . Для решения задачи Коши для уравнения (6.10) при t>0 введем равномерную сетку с постоянным шагом  . Введем понятие линейного m шагового разностного метода для решения задачи (6.10). Линейным m-шаговым разностным методом называется система разностных уравнений , (6.11) где: n=m,m+1...; – числовые коэффициенты, не зависящие от n; k=0,1,…,m, причем a0≠0. Систему (6.11) будем рассматривать как рекуррентные соотношения, выражающие новые значения через ранее найденные значения , причем расчет начинают с индекса n=m, т.е. с уравнения . Отсюда следует, что для начала расчета по формулам (6.11) надо знать m предыдущих значений функции y, причем y0=u0 определяется исходной задачей (6.10). Эти предыдущие m значений могут быть найдены одним из одношаговых методов Рунге-Кутта. Отличие от одношаговых методов состоит в том, что по формулам (6.11) расчет ведется только в точках сетки. Определение. Метод (6.11) называется явным, если коэффициент b0=0. Тогда значение легко выражается через . В противном случае метод называется неявным, и для нахождения yn придется решать нелинейное уравнение вида . (6.12) Обычно это уравнение решают методом Ньютона при начальном значении . Коэффициенты уравнения (6.11) определены с точностью до множителя, тогда, чтобы устранить этот произвол, вводят условие , с тем условием, что правая часть (6.11) аппроксимирует правую часть дифференциального уравнения (6.10). На практике используют частный случай методов (6.11), т.н. методы Адамса, когда производная аппроксимируется разностным отношением, включающим две соседние точки и . Тогда ; =0, k=2,...,m и . (6.13) Это и есть методы Адамса. При b0=0 метод будет явным, в противном случае – неявным. 6.2.1 Задача подбора числовых коэффициентов aк , bк Выясним, как влияют коэффициенты ak, bk на погрешность аппроксимации уравнения (6.11), на устойчивость и сходимость. Определение. Невязкой, или погрешностью аппроксимации методов (6.11) называется функция , (6.14) которая получается в результате подстановки точного решения u(t) дифференциального уравнения (6.10) в разностное уравнение (6.11) . Если разложить функции в ряд Тейлора в точках равномерной сетки, окончательно получим функцию . (6.15) Из вида функции следует, что порядок аппроксимации будет равен p, если выполнены условия (6.16) где l=1,...,p. Условия (6.16) представляют собой систему линейных уравнений из p+2 уравнений относительно неизвестных a0,…,am, b0,…,bm. Их количество равно 2(m+1). Решив систему (6.16), получаем неизвестные числовые коэффициенты. Для неявных m-шаговых методов наивысшим достижимым порядком аппроксимации будет p=2m, а для явных – p=2m-1. Запишем систему (6.16) для методов Адамса (6.17) где l=2,...,p . Отсюда видно, что наивысший порядок аппроксимации для неявного m-шагового метода Адамса p=m+1, для явного (b0=0) – p=m. 6.2.2 Устойчивость и сходимость многошаговых разностных методов Наряду с системами уравнений (6.11) будем рассматривать т.н. однородные разностные уравнения вида , (6.18) где n=m,m+1,... . Будем искать его решение в виде функции , где q – число, подлежащее определению. Подставив в (6.18) получаем уравнение для нахождения q . (6.19) Уравнение (6.19) принято называть характеристическим уравнением для разностных методов (6.11). Говорят, что разностный метод (6.11) удовлетворяет условию корней, если все корни уравнения (6.19) q1, …, qm лежат внутри или на границе единичного круга комплексной плоскости, причем на границе нет кратных корней. Разностный метод (6.11), удовлетворяющий условию корней, называется устойчивым методом. Теорема. Пусть разностный метод (6.11) удовлетворяет условию корней и выполнено условие при . Тогда при , nm и достаточно малых  будет выполнена оценка , (6.20) где: rk− погрешность аппроксимации; − погрешность в задании начальных условий; M –константа , зависящая от L,T и не зависящая от n. Методы Адамса удовлетворяют условию корней, т.к. для них a0=-a1=1, следовательно, q=q1=1. 6.2.3 Примеры m-шаговых разностных методов Адамса Явные методы. При m=1 порядок точности p=1. Тогда метод описывается формулой . В этом случае получаем метод Эйлера. При m=2 порядок точности p=2. Тогда метод описывается формулой . При m=3 порядок точности p=3. Тогда метод описывается формулой . При m=4 порядок точности p=4. Метод описывается формулой . Неявные m-шаговые методы Адамса: m=1, p=2, – метод трапеций; m=2, p=3, ; m=3, p=4, . Неявные методы содержат искомое значение yn нелинейным образом, поэтому для его нахождения применяют итерационные методы решения нелинейных уравнений. 6.3 Численное интегрирование жестких систем обыкновенных дифференциальных уравнений Жесткие системы можно сравнить с плохо обусловленными системами алгебраических уравнений. Рассмотрим систему дифференциальных уравнений(ДУ) , (6.21) где . Для решения (6.21) рассмотрим разностные методы вида , (6.22) где n= m, m+1, m+2,…. Устойчивость и сходимость методов (6.22) определяется расположением корней характеристического уравнения, т.е. |q|1 – корни принадлежат единичному кругу. Среди методов (6.22) выделим те, которые позволяют получить асимптотически устойчивые решения. Пример. В качестве частного случая (6.21) рассмотрим уравнение вида , (6.23) где: ; <0; – решение ДУ. При  <0 решение есть монотонно убывающая функция при t. Для этого решения можно записать при любом шаге  >0 , (6.24) что означает устойчивость решения u(t). Рассмотрим для задачи (6.23) метод Эйлера , где: n=0,1,2,…, , q – промежуточный параметр, равный 1+. Оценка (6.24) для метода Эйлера будет выполнена тогда и только тогда, когда |q|. Шаг  лежит в интервале 0 <  ≤ . Метод Эйлера для задачи (6.23) устойчив при выполнении этого условия. Определение 1. Разностный метод (6.22) называется абсолютно устойчивым, если он устойчив при любом  >0. Определение 2. Разностный метод называется условно устойчивым, если он устойчив при некоторых ограничениях на шаг . Например, метод Эйлера для (6.23) условно устойчив, при 0 <  ≤ . Примером абсолютно устойчивого метода является неявный метод Эйлера , т.к. в этом случае, при любых >0. Замечание. Условная устойчивость является недостатком явных методов в связи с тем, что приходится выбирать мелкий шаг интегрирования. Пример для задачи (6.23). Если =−200, тогда  0.01. Если мы рассмотрим интервал (0,1], то необходимо будет 100 шагов. Неявные методы со своей стороны приводят к решению на каждом шаге нелинейного уравнения, но это уже недостаток неявного метода. 6.3.1 Понятие жесткой системы ОДУ Замечание. Все вышерассмотренные методы легко реализуются на примере одного уравнения и легко переносятся на системы ДУ, но при решении систем возникают дополнительные трудности, связанные с разномасштабностью описанных процессов. Рассмотрим пример системы двух уравнений: , где: t >0; a1,a2>0. Эта система однородных независимых ДУ имеет решение . Это решение монотонно убывает с ростом t. Пусть коэффициент а2 на порядок больше а1, т.е. а2>>a1. В этом случае компонента u2 затухает гораздо быстрее u1, и тогда, начиная с некоторого момента времени t , решение задачи u(t) почти полностью будет определяться поведением компоненты u1. Однако при численном решении данной задачи шаг интегрирования  будет определяться, как правило, компонентой u2, не существенной с точки зрения поведения решения системы. Рассмотрим метод Эйлера для решения данной системы . Он будет устойчив, если на шаг  наложены ограничения . Учитывая, что >>>0, получаем окончательное ограничение на  . Такие трудности могут возникнуть при решении любых систем ОДУ. Рассмотрим в качестве примера систему , (6.25) где А– квадратная матрица m*m. Если матрица А имеет большой разброс собственных чисел, то возникают проблемы с разномасштабностью описываемых системой процессов. Допустим, что матрица А постоянна (т.е. не зависит от t). Тогда система (6.21) будет называться жесткой, если: 1) вещественные части собственных чисел для всех k, где k=1,…,m; 2) число велико (десятки и сотни), и число S называется числом жесткости системы. Если же матрица А зависит от t, то и собственные числа зависят от t и зависят от t. Решение жесткой системы (6.25) содержит как быстро убывающие, так и медленно убывающие составляющие. Начиная с некоторого t >0 решение системы определяется медленно убывающей составляющей. При использовании явных разностных методов быстро убывающая составляющая отрицательно влияет на устойчивость, поэтому приходится брать шаг интегрирования  слишком мелким. 6.3.2 Некоторые сведения о других методах решения жестких систем На практике разностные методы (6.22) для решения жестких систем используются в виде методов Гира (неявный разностный метод) и метода матричной экспоненты (метод Ракитского). 6.3.2.1 Методы Гира Это частный случай методов (6.22), когда коэффициент , . Запишем числовые коэффициенты, которые определяются из условия p-го порядка точности аппроксимации системы разностными методами ; ; (6.26) где l=2,...,p. Решив систему линейных уравнений (6.26) с учетом предыдущих условий, получаем все нужные коэффициенты. Трехшаговый метод Гира (частный случай методов (6.22) с учетом условий (6.26)) имеет вид (6.27) Метод имеет третий порядок точности. При m=4, получаем четырехшаговый метод Гира (6.28) Запишем систему (6.26) в виде . (6.29) Решив (6.29) для каждого случая можем найти коэффициенты ak, k=1,2,…,т. Чисто неявные разностные методы обладают хорошими свойствами устойчивости, поэтому используются для решения жестких систем уравнений. 6.3.2.2 Метод Ракитского (матричной экспоненты) решения систем ОДУ , (6.30) где: ; ; А– матрица размерности n*n. Допустим, что матрица А – постоянная, т.е. ее элементы не зависят от времени. Система (6.30) – однородная, с постоянными коэффициентами. Запишем аналитическое решение (6.30) , (6.31) где – матричная экспонента и +…. (6.32) Проинтегрируем уравнение (6.30) при значениях t=, 2, 3,…. Если точно знать матрицу , то точное решение в указанных точках можно получить по формуле (6.31), т.е. решение можно записать Таким образом, задача сводится к тому, чтобы достаточно точно знать матрицу. На практике поступают следующим образом: при больших  нельзя воспользоваться рядом Тейлора в связи с его бесконечностью, т.е. для удовлетворительной точности пришлось бы взять много членов ряда, что трудно. Поэтому поступают так: отрезок [0,] разбивают на k частей, чтобы длина h=/k удовлетворяла условию ||A∙h||<0.1. Тогда запишем по схеме Горнера . Каждый столбец матрицы− вычисляют по формуле , где – вектор столбец, в i-ой строке которого 1, а в остальных – нули. Если эта матрица найдена, то решение находится по (6.31). Для исследования разностных методов при решении жестких систем рассматривают модельное уравнение , (6.33) где – произвольное комплексное число. Для того, чтобы уравнение (6.33) моделировало исходную систему (6.30) его нужно рассматривать при таких значениях , которые являются собственными числами матрицы А. Многошаговые разностные методы (6.31) имеют вид , (6.34) где: n=m, m+1…; ∙ . Если решение уравнения (6.34) искать в виде , то для нахождения числа q получим характеристическое уравнение вида . Для устойчивости метода достаточно выполнения условия корней . В случае жестких систем используются более узкие определения устойчивости. Предварительные сведения. Областью устойчивости разностных методов называется множество всех точек комплексной плоскости ∙, для которых разностный метод применительно к уравнению (6.33) устойчив. Определение 1. Разностный метод называется А-устойчивым, если область его устойчивости содержит левую полуплоскость Re<0. Замечание. Решение модельного уравнения (6.33) асимптотически устойчиво при значениях Re<0, поэтому сущность А-устойчивого метода заключается в том, что А-устойчивый разностный метод является абсолютно устойчивым, если устойчиво решение исходного дифференциального уравнения. Так как класс А-устойчивых методов узок, то пользуются А()-устойчивым методом. Определение 2. Разностный метод (6.31) называется А()-устойчивым, если область его устойчивости содержит угол меньший , т.е. |arg(-)|<, где ∙. Исходя из этого, определяется, что при  А() устойчивость совпадает с определением А-устойчивого метода. 6.4 Краевые задачи для обыкновенных дифференциальных уравнений Постановка краевой задачи. Рассматриваем дифференциальное уравнение порядка n (n2) . (6.35) Если сделаем замену переменных вида (6.36) где: ; k=0,1,…,(n-1), то задача (6.35) сводится к задаче Коши для нормальной системы ОДУ порядка n. Типовые примеры краевых задач. Рассмотрим дифференциальное уравнение . (6.37) Для уравнения (6.37) краевая задача формулируется следующим образом: найти решение y=y(x), удовлетворяющее уравнению (6.37), для которой значения ее производных в заданной системе точек удовлетворяют n независимым краевым условиям, в общем виде нелинейным. Эти краевые условия связывают значения искомой функции y и ее производных до (n-1) порядка на границах заданного отрезка. 1. Рассмотрим уравнение второго порядка . Необходимо найти решение уравнения, удовлетворяющее заданным краевым условиям: y(a)=A, y(b)=B, т.е. необходимо найти интегральную кривую, проходящую через две заданные точки (рисунок 7). 2. Рассмотрим уравнение с краевыми условиями , . Из графика на рисунке 8 видно, что tg()=A1, tg()=B1. Здесь интегральная кривая пересекает прямые x=a и x=b под заданными углами  и  соответственно. 3. Смешанная краевая задача включает оба случая. Рассматривается то же самое уравнение , но с краевыми условиями y(a)=A1,(b)=B1 . Геометрическую иллюстрацию этих краевых условий легко представить, используя рисунки 7 и 8. Замечание. Краевая задача для уравнения (6.37) в общем случае может не иметь решений, иметь единственное решение, иметь несколько решений или бесконечное множество решений. 4. Поражение заданной цели баллистическим снарядом. Дифференциальные уравнения движения снаряда с учетом сопротивления воздуха имеют вид , где: – вторая производная по времени; E=E(y,v) – известная функция высоты и скорости; ; g=g(y) – ускорение силы тяжести; – угол наклона к горизонту касательной к траектории движения снаряда; . Предполагая, что при t=t0 снаряд выпущен из точки, совпадающей с началом координат с начальной скоростью v0 под углом 0 , а в момент t=t1 он поразит неподвижную мишень в точке получаем краевые условия Здесь неизвестны значения 0 и t1. Решив данную краевую задачу, можем найти начальный угол , где: ; 0 – угол, при котором поражается цель в точке M. Аналогично ставятся краевые задачи для систем дифференциальных уравнений. 6.5 Решение линейной краевой задачи Рассмотрим важный частный случай решения краевой задачи, когда дифференциальное уравнение и краевые условия линейны. Для этого рассмотрим уравнение , (6.38) где: и f(x) известные непрерывные функции на отрезке [a, b]. Предположим, что в краевые условия входят две абсциссы x=a, x=b. Это двухточечные краевые задачи. Краевые условия называются линейными, если они имеют вид =, (6.39) где:    − заданные константы. Причем они одновременно не равны нулю, т.е. , при v=1,2,…,n. Например, краевые условия во всех трех рассмотренных ранее задачах линейны, т.к. их можно записать в виде , причем – для первой задачи. 6.6 Решение двухточечной краевой задачи для линейного уравнения второго порядка сведением к задаче Коши Запишем линейное уравнение второго порядка в виде , (6.40) где: p, q, f − известные непрерывные функции на некотором отрезке [a,b]. Требуется найти решение уравнения (6.40), удовлетворяющее заданным краевым условиям (6.41) Причем константы  и  одновременно не равны нулю Решение задачи (6.40), (6.41) будем искать в виде линейной комбинации y=Cu+V, где С – константа, u – общее решение соответствующего однородного уравнения , (6.42) а V– некоторое частное решение неоднородного уравнения . (6.43) Потребуем ,чтобы первое краевое условие было выполнено при любом C, , откуда следует, что , . Тогда , (6.44) где k– некоторая константа, не равная нулю. Значение функции V и ее производная в точке а могут быть, например,выбраны равными , (6.45) если коэффициент α0 и , (6.46) если коэффициент α1. Из этих рассуждений следует, что функция u – есть решение задачи Коши для однородного уравнения (6.42) с начальными условиями (6.44), а функция V – есть решение задачи Коши для неоднородного уравнения (6.43) с начальными условиями (6.45) или (6.46) в зависимости от условий. Константу C надо подобрать так, чтобы выполнялись условия (6.41) (вторая строчка) в точке x=b . Отсюда следует, что , где знаменатель не должен быть равен нулю, т. е. . (6.47) Если условие (6.47) выполнено, то краевая задача (6.35), (6.36) имеет единственное решение. Если же (6.47) не выполняется, то краевая задача (6.35), (6.36) либо не имеет решения, либо имеет бесконечное множество решений. 6.7 Методы численного решения двухточечной краевой задачи для линейного уравнения второго порядка 6.7.1 Метод конечных разностей Рассмотрим линейное дифференциальное уравнение (6.48) с двухточечными краевыми условиями ; , (6.49) где: p, q, f – известные непрерывные функции на некотором отрезке [a,b]. Одним из наиболее простых методов решения этой краевой задачи является сведение ее к системе конечно-разностных уравнений. Основной отрезок [a,b] делим на n равных частей с шагом h=(b-a)/n, то есть рассматриваем равномерную сетку , i=0,1,…,n. В каждом внутреннем узле дифференциальное уравнение (6.48) аппроксимируем, используя формулы численного дифференцирования второго порядка точности (6.50) где i=1,..., n-1. Для граничных точек и , чтобы не выходить за границы отрезка, используем формулы численного дифференцирования первого порядка точности . (6.51) Используя отношение (6.50) исходное дифференциальное уравнение (6.48) аппроксимируем конечно-разностными уравнениями , (6.52) где i=1,...,n-1. Учитывая краевые условия, получим еще два уравнения . (6.53) Таким образом, получена линейная система n+1 уравнений с n+1 неизвестными , представляющими собой значения искомой функции в точках . Разностная схема (6.52)-(6.53) аппроксимируеткраевую задачу (6.48)- (6.49) с порядком 1 по h за счет краевых условий. Решив эту систему, получим таблицу значений искомой функции y. 6.7.2 Метод прогонки (одна из модификаций метода Гаусса) При применении метода конечных разностей к краевым задачам для дифференциальных уравнений второго порядка получается система линейных алгебраических уравнений с трехдиагональной матрицей, т.е. каждое уравнение системы содержит три соседних неизвестных. Для решения таких систем разработан специальный метод – «метод прогонки». Для этого систему (6.52) перепишем в виде (6.54) для внутренних точек (i=1,…,n-1), где: ; ; . На концах отрезка x0=a и xn=b производные заменяем разностными отношениями и . Учитывая эту замену, получим еще два уравнения (6.55) Обратим внимание на внешний вид записи системы (6.54), (6.55). В каждом уравнении системы присутствует три ненулевых элемента. В первом и последнем – по два ненулевых коэффициента. Разрешая уравнение (6.54) относительно yi, получим . (6.56) Предположим, что с помощью полной системы (6.54), (6.55) из уравнения (6.56) исключена неизвестная yi-1. Тогда это уравнение примет вид , (6.57) где: – некоторые коэффициенты; i=1,2,…,n-1. Отсюда . Подставляя это выражение в уравнение (6.54), получим , а отсюда . (6.58) Сравнивая (6.57) и (6.58), получим для определения ci и di рекуррентные формулы ; i=1,…,n-1. (6.59) Из первого краевого условия (6.55) и из формулы (6.57) при i=0 находим ;. (6.60) На основании формул (6.59), (6.60) последовательно определяются коэффициенты ci, di (i=1,…,n-1) до cn-1 и dn-1 включительно (прямой ход). Обратный ход начинается с определения . Для этого из второго краевого условия (6.55) и из формулы (6.57) при i=n-1 найдем . (6.61) Далее по формуле (6.57) последовательно находим . Заметим, что метод прогонки обладает устойчивым вычислительным алгоритмом. 7 Приближенное решение дифференциальных уравнений в частных производных В реальных физических процессах искомая функция зависит от нескольких переменных, а это приводит к уравнениям в частных производных от искомой функции. Как и для обыкновенных дифференциальных уравнений(ОДУ), в этом случае для выбора одного конкретного решения, удовлетворяющего уравнению в частных производных, кроме начальных условий, необходимо задавать дополнительные условия (т.е. краевые условия ). Чаще всего такие задачи на практике не имеют аналитического решения и приходится использовать численные методы их решения, в том числе метод сеток, метод конечных разностей и так далее. Мы будем рассматривать класс линейных уравнений в частных производных второго порядка. В общем виде в случае двух переменных эти уравнения записываются в виде , (7.1) где: A, B, C, a, b, c − заданные непрерывные функции двух переменных, имеющие непрерывные частные производные, u – искомая функция. Для сокращения записи введем обозначения ; ; ; ; . Будем рассматривать упрощенную форму записи (7.1) вида (7.2) и рассмотрим частный случай (7.2), когда a=b=c=F0, т.е. . (7.3) Путем преобразований уравнение (7.3) может быть приведено к каноническому виду (к одному из трех стандартных канонических форм) эллиптическому типу, гиперболическому типу, параболическому типу. Причем тип уравнения будет определяться коэффициентами А, В, С, а именно – знаком дискриминанта D=B2−4AC. Если D <0, то имеем уравнение эллиптического типа в точке с координатами x, y; если D=0, то (7.3) − параболического типа; если D>0, то (7.3) − гиперболического типа; если D не сохраняет постоянного знака, то (7.3) − смешанного типа. Замечание. Если А, В, С − константы, тогда каноническое уравнение (7.3) называется полностью эллиптического, параболического, гиперболического типа. Введем понятие оператора Лапласа для сокращенной записи канонических уравнений вида . Используя это определение, запишем сокращенные канонические уравнения всех трех типов 1. u=0. Это уравнение эллиптического типа, так называемое уравнение Лапласа. В механике это уравнение описывает стационарные тепловые поля, установившееся течение жидкости и т.д. 2.u=−f , где f − заданная непрерывная функция. Это уравнение Пуассона имеет эллиптический тип и описывает процесс теплопередачи с внутренним источником тепла. 3. , где a − константа. Не во всех уравнениях в качестве переменных будут выступать стандартные переменные x, y. Может быть также переменная времени. Это уравнение диффузии описывает процесс теплопроводности и является уравнением параболического типа. 4. , а– константа. Это уравнение гиперболического типа − так называемое волновое уравнение и оно описывает процесс распространения волн. 7.1 Метод сеток для решения смешанной задачи для уравнения параболического типа (уравнения теплопроводности) Смешанная задача означает, что следует найти искомую функцию, удовлетворяющую заданному уравнению в частных производных, краевым, а так же начальным условиям. Различить эти условия можно в том случае, если одна из независимых переменных − время, а другая − пространственная координата. При этом условия, относящиеся к начальному моменту времени, называются начальными, а условия, относящиеся к фиксированным значениям координат − краевыми. Рассмотрим смешанную задачу для однородного уравнения теплопроводности , k=const>0. (7.4) Задано начальное условие (7.5) и заданы краевые условия первого рода (7.6) Требуется найти функцию u(x,t), удовлетворяющую в области D (00. Есть и явная схема (рисунок 11), но она устойчива только при , т.е. при . Вычисления по этой схеме придется вести с малым шагом по , что приводит к большим затратам машинного времени. Рисунок 11 – Четырехточечный шаблон явной схемы 7.2 Решение задачи Дирихле для уравнения Лапласа методом сеток Рассмотрим уравнение Лапласа . (7.8) Будем рассматривать уравнение Лапласа в прямоугольной области с краевыми условиями ; ; ; , где − заданные функции. Заметим, что чаще всего область бывает не прямоугольной. Введем обозначения uij=u(xi,yj). Накладываем на прямоугольную область сетку ; i=0,1,…,n; ; j=0,1,…,m. Тогда , =b. Частные производные аппроксимируем по формулам и заменим уравнение Лап­ласа конечно-разностным уравнением Рисунок 12 – Схема узлов «крест» , (7.9) где: i=1,…,n-1, j=1,...,m-1 (т.е. для внутренних узлов). Погрешность замены дифференциального уравнения разностным составляет величину О(h2+l2). Уравнения (7.9) и значения ui,j в граничных узлах образуют систему линейных алгебраических уравнений относительно приближенных значений функции u(x,y) в узлах сетки. Выразим ui,j при h=l, и заменим систему (7.10) Систему (7.10) линейных алгебраических уравнений можно решить любым итерационным методом (Зейделя, простых итераций и т.д.). При построении системы использовалась схема типа «крест» (рисунок 12). Строим последовательность итераций по методу Зейделя , где s– текущая итерация. Условие окончания итерационного процесса . (7.11) Условие (7.11) ненадежно и на практике используют другой критерий где . Схема «крест»– явная устойчивая схема ( малое изменение входных данных ведет к малому изменению выходных данных). 7.3 Решение смешанной задачи для уравнения гиперболического типа методом сеток Рассмотрим уравнение колебания однородной и ограниченной струны. Задача состоит в отыскании функции u(x,t) при t>0, удовлетворяющей уравнению гиперболического типа , (7.12) где: 0