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

Введение Математическое моделирование и вычислительный эксперимент 1. Схема вычислительного эксперимента. Эффективное решение крупных естественнонаучных и народнохозя йственных задач сейчас невозможно без применения быс тродействующих электронно-вычислительных машин (ЭВМ). В настоящее время выработалась технология и сследования сложных проблем, основанная на построении и анализе с помощью ЭВМ математических моделей изуча емого объекта. Такой метод исследования называют вычислительным экспериментом. Пусть, например, требуется исследовать какой -то физический объект, явление, процесс. Тогда схема вычисл ительного эксперимента выглядит так, как показано на р исунке 1. Формулируются основные законы, управляющие данным объектом исследования ( I) и строится соответствующая математическая модель (II), представляющая обычно запись этих законов в форме системы уравнений (алгебраических, дифференциальных, интегральных и т. д.). I II Объект исследования Проведение вычислений и анализ результатов Математическая модель V Численные III методы(дискретная модель вычислительного алгоритма) IV Программирование на ЭВМ Рисунок 1 – Этапы построения и анализа с помощью ЭВМ математической модели объекта 1 При выборе физической и, следовательно, математич еской модели мы пренебрегаем факторами, не оказывающ ими существенного влияния на ход изучаемого процесса. Типичные математические модели, соответствующие физ ическим явлениям, формулируются в виде уравнений математической физики. Большинство реальных процессов оп исывается нелинейными уравнениями и лишь в первом пр иближении (при малых значениях параметров, малых откл онениях от равновесия и др.) эти уравнения можно заменить линейными. После того как задача сформулирована в математич еской форме, необходимо найти ее решение. Но что значит решить математическую задачу? Только в исключител ьных случаях удается найти решение в явном виде, напр имер в виде ряда. Иногда утверждение «задача решена» означ ает, что доказано существование и единственность решения. Ясно, что этого недостаточно для практических прилож ений. Необходимо еще изучить качественное поведение р ешения и найти те или иные количественные характерист ики. Именно на этом этапе требуется прив лечение ЭВМ и, как следствие, развитие численных методов (см. III на рис. 1). Под численным методом здесь понимается такая инте рпретация математической модели («дискретная м одель»), которая доступна для реализации на ЭВМ. Н апример, если математическая модель представляет собой дифференц иальное уравнение, то численным методом может быть а ппроксимирующее его разностное уравнение совместно с а лгоритмом, позволяющим отыскать решение этого разнос тного уравнения. Результатом реализации численного мет ода на ЭВМ является число или таблица чисел. Отметим, что в настоящее время помимо собственно численных м етодов имеются также методы, которые позволяют пров одить на ЭВМ аналитические выкладки. Однако аналитич еские методы для ЭВМ не получили пока достаточно шир окого распространения. Чтобы реализовать численный метод, необходимо с оставить программу для ЭВМ (см. IV на рис. 1) или воспользоваться готовой программой. 2 После отладки программы наступает этап проведения вычислений и анализа результатов ( V). Полученные результаты изучаются с точки зрения их соответствия иссл едуемому явлению и, при необходимости, вносятся испра вления в численный метод и уточняется математическая м одель. Такова в общих чертах схема вычислительного экспе римента. Его основу составляет триада: модель — метод (алгоритм) — программа. Опыт решения крупных задач показывает, что метод математического моделирования и вычислительный эксперимент соединяют в себе преимущества традиционных теоретических и экспериментальных методов исследования. Можно указать так ие крупные области применения вычислительного эксперимента, как энергетика, аэрокосмическая техника, обработка дан ных натурного эксперимента, совершенствование технологич еских процессов. 2. Вычислительный алгоритм. Предметом данной книги является изложение вопросов, отражающих этапы III, IV, V вычислительного эксперимента. Т аким образом, здесь не обсуждаются исходные задачи и их математиче ская постановка. Необходимо подчеркнуть, что процесс исследования исходного объекта методом математического моделиров ания и вычислительного эксперимента неизбежно носит приближенный характер, потому что на каждом этапе вн осятся те или иные погрешности. Так, построение математ ической модели связано с упрощением исход ного явления, недостаточно точным заданием коэффициентов уравнения и других входных данных. По отношению к численному методу, реализующему данную математическую модель, указанные погрешности являются неустранимыми, поскольку они неизбежны в рамках данной м одели. При переходе от математической модели к численн ому методу возникают погрешности, называемые погрешностями метода. Они связаны с тем, что всякий числе нный метод воспроизводит исходную математическую модель приближенно. 3 Наиболее типичными погрешностями метода являются погрешность дискретизации и погрешность округления. Поясним причины возникновения таких погрешностей. Обычно построение численного метода для заданной математической модели разбивается на два этапа: а) фор мулирование дискретной задачи, б) разработка вычисли тельного алгоритма, позволяющего отыскать решение дискретной задачи. Например, если исходная математ ическая задача сформулирована в виде системы дифференциальных уравнений, то для численного решения необходимо зам енить ее системой конечного, может быть, очень большого числа линейных или разностных алгебраических уравнений. В этом случае говорят, что проведена дискретизация исходной математической задачи. Простейшим примером дискретизации является построение разностной схемы, путем замены дифференциальных выражений конечноразностными отношениями. В общем случае дискретную модель можно рассматривать как конечномерный аналог исходной математической задачи. Ясно, что ре шение дискретизированной задачи отличае тся от решения исходной задачи. Разность соответствующих решений и называется погрешностью дискретизации. Обычно дискретная модель зависит от некоторого параметра (или множества параме тров) дискретизации, при стремлении которого к нулю должна стремиться к нулю и погрешность дискретизации. При этом число алгебраических уравнений, составляющих дискретную модель, неограниченно возрастает. В случае разностных методов таким параметром я вляется шаг сетки. Как уже отмечалось, дискретная модель представляет собой систему большого числа алгебраических уравн ений. Невозможно найти решение такой системы точно и в явном виде. Поэтому приходится использовать тот или иной чи сленный алгоритм решения системы алгебраических уравн ений. Входные данные этой систе мы, а именно коэффициенты и правые части, задаются в ЭВМ не точно, а с округ лением. В процессе работы алгоритма погрешности округл ения обычно накапливаются, и в результате решение, по лученное на ЭВМ, будет отличаться от точного реш ения дискре4 тизированной задачи. Результирующая погрешность наз ывается погрешностью округления (иногда ее называют вычислительной погрешностью). Величина этой погрешности определяется двумя факт орами: точностью представления вещественных чисел в ЭВМ и чувствительностью данного алгоритма к погрешн остям округления. Алгоритм называется устойчивым, если в процессе его работы вычислительные погрешности возрастают незначительно, и неустойчивым — в противоположном случае. При использовании неустойчивых вычислительных алгоритмов накопление погрешностей округления прив одит в процессе счета к переполнению арифметического устройства ЭВМ. Итак, следует различать погрешности модели, м етода и вычислительную. Какая же из этих трех погрешностей я вляется преобладающей? Ответ здесь неоднозначен. Ви димо, типичной является ситуация, возникающая при ре шении задач математической физи ки, когда погрешность модели значительно превышает погреш ность метода, а погрешностью округления в случае устойчивых алгори тмов можно пренебречь по сравнению с погрешностью метода. С другой стороны, при решении, например, с истем обыкновенных дифференциальных уравнений возможн о применение столь точных методов, что их погрешность будет сра внима с погрешностью округления. В общем случае нужно стремиться, чтобы все указанные погре шности имели один и тот же порядок. Например, нецелес ообразно пользоваться разностными схемами, имеющими точность 10 6 , если коэффициенты исходных уравнений задаются с точностью 10 2 . 3. Требования к вычислительным методам. Одной и той же математической задаче можно поставить в соотве тствие множество различных дискретн ых моделей. Однако далеко не все из них пригодны для практической реализ ации. Вычислительные алгоритмы, предназначенные для б ыстродействующих ЭВМ, должны удовлетворять многообразным и зачастую противоречивым требованиям. 5 Попытаемся здесь сформулировать ос новные из этих требований в общих чертах. Можно выделить две группы требований к численным методам. Первая группа связана с адекватностью дискре тной модели исходной математической задаче, и вторая группа – с реализуемостью численного метода на ЭВМ. К первой группе относятся такие тре бования, как сходимость численного метода, выполнение дискрет ных аналогов законов сохранения, качественно правильное пове дение решения дискретной задачи. Поясним эти требования. Предположим, что дискре тная модель математической задачи представляет собой систему большого, но конечного числа алгебраических уравнений. Обычно, чем точнее мы хотим получить решение, тем больше уравнений приходится брать. Говорят, что числе нный метод сходится, если при неограниченном увеличении числа уравнений решение дискретной задачи стремится к решению исходной задачи. Поскольку реальная ЭВМ может оперировать лишь с конечным числом уравнений, на практике сходимость, как правило, не достигается. Поэтому важно уметь оценивать погрешность метода в зависимости от числа уравнений, с оставляющих дискретную модель. По этой же причине стараются строить дискретную модель таким образом, чтобы она правильно отражала качественное пов едение решения исходной задачи даже при сравнительно небольшом числе уравнений. Например, дискретной моделью задачи математич еской физики может быть разностная схема. Для ее построения область изменения независимых переменных заменяется дискретным множеством точек – сеткой, а входящие в исходное уравнение производные заменяются, на сетке, конечно-разностными отношениями. В результате получаем систему алгебраических уравнений относительно значений искомой функции в точках сетки. Число уравнений этой системы равно числу точек се тки. Известно, что дифференциальные уравнения математической физики являются следствиями интегральных зак онов сохранения. Поэтому естественно требовать, чтобы для 6 разностной схемы выполнялись ан алоги таких законов сохранения. Разностные схемы, удовлетворяю щие этому требованию, называются консервативными. Оказалось, что при одном и том же числе точек сетки консервативные ра зностные схемы более правильно отражают поведение р ешения исходной задачи, чем неконсервати вные схемы. Сходимость численного метода тесно связана с его ко рректностью. Предположим, что исходная мате матическая задача поставлена корректно, т.е. ее решение сущ ествует, единственно и непрерывно зависит от входных данных. Т огда дискретная модель этой задачи должна быть пост роена таким образом, чтобы свойство корректности сохранилось. Таким образом, в понятие корректности численного метода включаются свойства однозначной разрешимости соо тветствующей системы уравнений и ее устойчи вости по входным данным. Под устойчивостью понимается непрерывная зависимость решения от входных данных, равн омерная относительно числа уравнений, составляющих ди скретную модель. Вторая группа требований, предъявляемых к числе нным методам, связана с возможностью реализации данной ди скретной модели на данной ЭВМ, т. е. с возможностью п олучить на ЭВМ решение соответствующей системы алг ебраических уравнений за приемлемое время. Основным пр епятствием для реализации корректно поставленного алгоритма является ограниченный объем оперативной п амяти ЭВМ и ограниченные ресурсы времени счета. Реальные в ычислительные алгоритмы должны учи тывать эти обстоятельства, т. е. они должны быть экономичными как по чи слу арифметических действий, так и по требуемому объему памяти. 7 Численные методы алгебры и анализа 1 Решение систем линейных алгебраических уравнений Рассмотрим систему линейных алгебраических уравнений: а 11 х1 а 21 х1 ..... а m1 х1 а 12 х 2 ... а 1m х m b1 а 22 х 2 ... а 2 m х m b2 .... а m 2 х 2 ... а mm х m bm (1.1) или в матричной форме: Aх=b, (1.2) где: A={a ij } квадратная матрица размерности (m m,); х=(х 1 ,…., х m ) T ; T – операция транспонирования; b=(b 1 ,….,b m ) T ; detA 0. Предположим, что определитель матрицы A не равен нулю. Тогда решение х существует и единственно. На практике встречаются системы, имеющие большой поря док. Методы решения системы (1.1) делятся на две группы: 1) прямые (точные методы); 2) итерационные методы (приближенные). 1.1 Точные методы В точных методах решение х находится за конечное число действий, но из-за погрешности округления и их накопления прямые методы можно назвать точными, только отвлекаясь от погрешностей округления. 1.1.1 Метод Гаусса Вычисления с помощью метода Гаусса (который наз ывают также методом последовательного исключения неи звестных) состоят из двух основных этапов: прямого хода и обратного хода. Прямой ход метода заключается в последовательном исключении неизвестных из системы для пр еобразования ее к эквивалентной системе с треугольной 8 матрицей. На этапе обратного хода производят вычисления значений неизвестных. Рассмотрим простейший вариант метода Гаусса, называемый схемой единс твенного деления. Прямой ход метода 1-й шаг. Предположим, что а 1 1 0. Поделим первое уравнение на этот элемент, который назовем ведущим элементом первого шага : x1 c12 x 2 ... c1m x m y1 . (1.3) Остальные уравнения системы (1.1) запишем в виде a i1 x1 a i 2 x 2 ... a im x m bi , (1.4) где i= 2, m . Уравнение (1.3) умножаем на a i1 и вычитаем из i-го уравнения системы (1.4). Это позволит обратить в нуль к оэффициенты при x 1 во всех уравнениях, кроме первого. Получим эквивалентную систему вида: x 1 c12 x 2 ... c1m x m y1 (1) (1) ( 1) a 22 x 2 ... a 2 m x1 b 2 . . . . . . . . . . ( 1) a m( 12) x 2 ... a mm x1 b m( 1 ) (1.5) a ij(1) a ij c1 j a i1 , bi(1) bi y1 a i1 где i,j= 2, m . Система (1.5) имеет матрицу вида: 1 x ... x x ... x ... ... ... ... x ... . x Дальше работаем с укороченной системой, т.к. х 1 входит только в 1-ое уравнение 9 (1) (1) (1) a 22 x 2 ... a 2 m x m b2 . . . (1) am2 x2 . . . . (1) ... a mm xm b m(1) . 2-й шаг. На этом шаге исключаем неизвестное х 2 из уравнений с номерами i=3,4,…,m. Если ведущий элемент (1) второго шага а 22 0 , то из укороченной системы ана логично исключаем неизвестное x 2 и получаем матрицу коэффициентов такого вида: 1 x x ... x 1 x ... x x ... x . ...  ...  ... x ... x Аналогично повторяем указанные действия для неи звестных х 3 ,х 4 ,...,х m-1 и приходим к системе : x1 c12 x2 ... c1m xm x2 ... c2m xm . . xm 1 . cm . . 1, m xm cm m xm y1 y2 . .. ym 1 ym (1.6) Эта система с верхней треугольной матрицей: 1 x x ... x x 1 x ... x x 1 x ... x       ... 1 x ... ... x . Обратный ход метода. Из последнего уравнения сис темы (1.6) находим х m , из предпоследнего х m-1 , ..., из первого уравнения – х 1 . 10 Общая формула для вычислений: x m =y m /c mm , m x i y i c x , (i=m-1,…,1). ij j j i 1 Для реализации метода Гаусса требуется примерно (1/3)m 3 арифметических операций, причем больши нство из них приходится на прямой ход. Ограничение метода единственного делен ия заключается в том, что ведущие элементы на k-ом шаге исключения k 1 не равны нулю, т.е. аkk ≠0. Но если ведущий элемент близок к нулю, то в проце ссе вычисления может накапливаться погрешность. В этом случае на каждом шаге исключают не х k , a х j (при j k). Такой подход называется методом выбора главного эл емента. Для этого выбирают неизвестные x j с наибольшим по абсолютной величине коэффициентом либо в строке, либо в столбце, либо во всей матрице. Для его реализации треб уется m(m 2 3m 1) − арифметических действий. 3 1.1.2 Связь метода Гаусса с разложением матрицы на множители. Теорема об LU разложении. Пусть дана система Aх=b (1.1), которая при прямом ходе преобразуется в эквивалентную систему (1.6) и з апишем ее в виде Cх=y, (1.6*) где С – верхняя треугольная матрица с единицами на гла вной диагонали, полученная из (1.6) делением последнего уравнения системы на с mm . Как связаны в системе (1.1) элементы b и элементы y из (1.6*)? Если внимательно посмотреть на прямой ход метода Гаусса, то можно увидеть, что b1 a11 y1 . (1) b2 a 21 y1 a 22 y2 11 Для произвольного j имеем bj d j1 y1 d j 2 y2 ... d jj yi , (1.7) где j= 1, m , d ji – числовые коэффициенты: d j j a (j jj 1) . (1.8) Можно записать систему: b=Dy, где D−нижняя треугольная матрица с элементами a (jjj 1) на (0) а11 ). главной диагонали (j= 1, m , а11 В связи с тем, что в методе Гаусса угловые коэ ффици( j 1) 0 , то на главной диагонали матенты не равны нулю a jj рицы D стоят не нулевые элементы. Следовательно, эта матрица имеет обратную, тогда y=D -1 b, Сx= D -1 b. Тогда 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. 12 Теорема об LU разложении Введем обозначения: j − угловой минор порядка j матрицы А, т.е. . 1 a11 , 2 det . m . a11 a12 a 21 a 22 . . , . det ( A ). Теорема. Пусть все угловые миноры матрицы А не равны нулю (Δ j 0 для j= 1, m ). Тогда матрицу А можно представить единственным образом в виде произведе ния А=L*U. Идея доказательства. Рассмотрим матрицу А второго порядка и будем искать разложение этой матрицы в виде L и U. A a11 a12 l11 a 21 a 22 l 21 l 22 * 1 u12 1 . Сопоставляя эти два равенства, определяем элементы матриц L и U (перемножим и приравняем неизвестные). Система имеет единственное решение. Методом математ ической индукции сказанное можно обобщить для ма трицы размерности m m. Следствие. Метод Гаусса (схему единственного деления) можно применять только в том случае, когда угл овые миноры матрицы А не равны нулю. 1.1.3 Метод Гаусса с выбором главного элемента Может оказаться так, что система (1.1) имеет еди нственное решение, хотя какой либо из миноров матрицы А равен нулю. Заранее неизвестно, что все угловые миноры матрицы А не равны нулю. В этом случае можно использовать метод Гаусса с выбором главного элемента. 1. Выбор главного элемента по столбцу, когда на k-ом шаге исключения в качестве главного элемента выбирают 13 максимальный по модулю коэффициент при неизвестном х k в уравнениях с номерами i=k,k+1,…,m. Затем уравнение, соответствующее выбранному коэффициенту, меняют ме стами с k-ым уравнением системы, чтобы ведущий элемент k 1 занял место коэффициента аkk .После перестановки исключение неизвестного х k выполняют как в схеме единственн ого деления. ПРИМЕР. Пусть дана система второго порядка a11 x1 a12 x 2 b1 a 21 x1 a 22 x 2 b2 Предположим, что а 2 1 > а 1 1 , тогда переставим уравнения a 21 x1 a 22 x 2 a11 x1 a12 x 2 b2 b1 и применяем первый шаг прямого хода метода Гаусса. В этом случае имеет место перенумерация строк. 2. Выбор главного элемента по строке, т.е. производится перенумерация неизвестных системы. При а 12 > а 11 , на первом шаге вместо неизвестного х 1 исключают х 2 : a12 x 2 a11 x1 b1 a 22 x 2 a 21 x1 b2 . К этой системе применяем первый шаг прямого хода метода Гаусса. 3. Поиск главного элемента по всей матрице з аключается в совместном применении методов 1 и 2. Всѐ это приводит к уменьшению вычислительной п огрешности, но может замедлить процесс решения задачи. 1.1.4 Метод Холецкого (метод квадратных корней) Пусть дана система Ах=b, (1.11) где А − симметричная положительно определенная матрица. 14 Тогда решение системы (1.11) проводится в два эт апа: 1. Симметричная матрица А представляется как произведение двух матриц А = L∙ L Т . Рассмотрим метод квадратных корней на примере сис темы 4-го порядка: a11 ... a14 ... ... ... a 41 ... a 44 l11 l 21 l 22 l31 l32 l33 l 41 l 42 l 43 l 44 * l11 l 21 l31 l 41 l 22 l32 l 42 l33 l34 l 44 . Перемножаем матрицы в правой части разложения и сравниваем с элементами в левой части: l11 l32 a11 , l 21 a12 a 32 l 21l31 l 22 a11 , l31 , l33 a13 a11 , l 41 a14 a11 a33 l312 l322 , l 44 , l 22 a 22 2 l 21 , 2 2 2 . a 44 l 41 l 42 l 43 2. Решаем последовательно две системы Ly=b, L T x=у. Замечания 1) Под квадратным корнем может получиться отрицательное число, следовательно в программе необходимо предусмотреть использование правил действия с комплексными числами. 2) Возможно переполнение, если угловые элементы близки к нулю. 1.2 Итерационные методы решений систем а лгебраических уравнений Итерационные методы обычно применяются для реш ения систем большой размерности и они требуют приведе ния исходной системы к специальному виду. 15 Суть итерационных методов заключается в том, что р ешение х системы (1.1) находится как предел последов ательности lim x(n) . n Так как за конечное число итераций пред ел не может быть достигнут, то задаѐтся малое число − точность, и последовательные приближения вычисляют до тех пор, п ока не будет выполнено неравенс тво xn xn 1 , где n=n( ) – функция , ||x|| − норма вектора. Определения основных норм в пространстве векторов и матриц. Для вектора x=( x 1 ,x 2 ,…,x n ) T нормы вычисляются по следующим формулам: x 1 max xi ; 1 i n n x 2 xi ; i 1 n x xi 3 2 . i 1 Согласованные с ними нормы в пространстве матриц: n A1 aij ; max 1 i n j 1 n A 2 aij max 1 j n n A 3 ; i 1 aij 2 − величина, называемая евклидовой нор- i, j 1 мой матрицы A. Прямые методы рассчитаны для решения систем, порядок которых не больше 100, иначе для практических вычислений используются итерационные методы. 1.2.1 Метод Якоби (простых итераций) Исходную систему (1.11) 16 Ах=b преобразуем к виду: i 1a ij xi j a 1 ii m xj j i aij xj a ii 1 bi , aii (1.12) где i=1,2,...,m; a ii 0. Первая сумма равна нулю, если верхний предел сум мирования меньше нижнего. Так (1.12) при i=1 имеет вид m x1 j По методу Якоби xin формуле xin 1 i 1a ij 1 j a 1 ii a1 j xj a 2 11 b1 . a11 (n+1 приближение х i ) ищем по m x nj j i aij n xj a 1 ii bi . aii (1.13) где n – номер итерации (0,1,…); i= 1, m . Итерационный процесс (1.13) начинается с начальных значений xi0 , которые в общем случае зад аются произвольно, но предпочтительнее за xi0 взять свободные члены исходной системы. Условие окончания счета: max x n 1 i i xn , i где i= 1, m . 1.2.2 Метод Зейделя Система (1.11) преобразуется к виду (1.12) и органи зуется итерационная процедура , где неизвестные х i на n+1 шаге определяются по формулам xin 1 i 1a ij a j 1 ii x nj 1 m j i aij n xj a 1 ii bi . aii (1.14) Например, 17 m x1n 1 x2n a1 j j 2 a11 m 1 j x nj a2 j n xj a 3 22 b1 a11 (1.15) , b2 a22 a21 n 1 x , a22 1 (1.16) и так далее. Итерационные процессы (1.13) и (1.14) сходятся, если норма матрицы А (А − матрица коэффициентов при неизвестных в правой части систем (1.13) и (1.14)) удовлетв оряет условию: ||A||<1. 1.2.3 Матричная запись методов Якоби и Зейделя Исходную матрицу системы (1.11) представим в виде суммы трѐх матриц A=A 1 +D+A 2 , где D − диагональная матрица; D =diаg[а 11 а 2 2 …а mm ]; A 1 − нижняя треугольная матрица; A 2 − верхняя треугольная матрица. Пример: Дана матрица размерности (3 3): 0 а12 а13 0 0 а23 0 0 0 0 0 0 а21 0 0 а31 а32 0 А1 а11 0 0 0 а22 0 0 0 а33 А2 А. D Тогда исходную систему (1.11) можно записать в виде x=–D -1 A 1 x – D - 1 A 2 x+D -1 b. Тогда метод Якоби можно записать в виде: xn 1 D 1 A1 x n D 1 A2 x n D 1b или D хn 1 ( A1 A2 ) x n b. (1.17) 18 В матричной форме метод Зейделя будет выгл ядеть: xn 1 D 1 A1 x n 1 D 1 A2 x n 1 А2 х n b . D 1b или ( D А1 ) х n (1.18) Преобразуем формулы (1.17) и (1.18): D( х n 1 x n ) Ах n b , ( D А1 )( х n 1 (1.19) x n ) Ax n b . (1.20) Из (1.19) и (1.20) видно, что если итерационный метод сходится, то он сходится к точному решению. Иногда при решении задач большой размерности, в итерационные м етоды вводятся числовые параметры, которые могут зав исеть от номера итерации. Пример для метода Якоби. D хn 1 tn xn Ax n b , 1 где t – числовой параметр. Возникают вопросы: 1) При каких значениях t сходимость будет наиболее быстрой? 2) При каких значениях t метод сходится? На примере двух методов просматривается вывод о том , что одни и те же методы можно записывать несколькими способами. Поэтому вводят каноническую (ста ндартную) форму записи: Dn xn 1 1 tn xn Ax n b . (1.21) 1 Формула (1.21) получена путем объединения (1.19) и (1.20). Матрица D n+1 здесь задает тот или иной метод. Если существует обратная матрица к этой матрице, то из п оследней системы мы можем найти все неизвестные. 1. Метод (1.21) – явный, если матрица D n совпадает с единичной матрицей и неявный – в противном случае. 19 2. Метод (1.21) – стационарный, если матрица D n +1 =D, и параметр t не зависит от номера итерации и нестационарный – в противном случае. 1.2.4 Метод Ричардсона Явный метод с переменным параметром t: xn 1 tn xn (1.21а) Ax n b , 1 называется методом Ричардсона. 1.2.5 Метод верхней релаксации (обобщѐнный метод Зейделя) ( D ω A1 ) xn 1 xn ω (1.21б) Ax n b , где – числовой параметр. Если матрица А − симметричная и положительно определена, то последний метод сходится при (0 < < 2). Последнюю формулу запишем в следующем в иде: ( E ω D 1 A1 ) x n 1 ((1 ω ) E ω D 1 A2 ) x n ω D 1b , (1.22) где Е – единичная матрица. Тогда для вычисления неизвестных х i (i= 1, m ) можно записать итерационную процедуру в виде: xin 1 ω i 1 j 1 a ij aii x nj 1 (1 ω ) xin ω m a ij j i 1 a ii x nj ω bi a ii . (1.23) Например, для х 1 это будет такое выражение: x1n 1 (1 ω ) x1n ω m j 2 a1 j a11 x nj ω b1 a11 . 1.2.6 Сходимость итерационных методов Рассмотрим систему Ax=b, где А – невырожденная действительная матрица. 20 Для решения системы рассмотрим одношаговый ст ационарный метод D xn 1 xn Ax n b , t (1.24) при n=0,1,2…. Предположим, что задан начальный вектор решения. Тогда метод (1.24) сходится, если норма вектора x xn 0. n Теорема. Условие сходимости итерационного метода. Пусть А – симметричная положительно определе нная матрица и выполнено условие D - 0.5tA > 0 (где t > 0). Тогда итерационный метод (1.24) сходится. Следствие 1. Пусть А – симметричная и положительно определенная матрица с диагональным преобладан ием, то есть m | aij | , | a jj | i 1 i j при j=1,2,…,m. Тогда метод Якоби сходится. Следствие 2. Пусть А – симметричная и положительно определенная матрица с диагональным преобладанием, т огда метод верхней релаксации сходится при (0< <2). Проверяется, при каком значении метод достигает заданной точности быстрее. В частности, при = 1 метод верхней релаксации превращается в метод Зейделя, следовательно, при = 1 метод Зейделя сходится. Теорема. Итерационный метод (1.24) сходится при любом начальном векторе x 0 тогда и только тогда, когда все собственные значения матрицы S E tD 1 A по модулю меньше единицы. 21 2 Плохо обусловленные системы линейных алге браических уравнений Дана система линейных алгебраических уравнений Ах=b (2.1) Если система плохо обусловлена, то это значит, что п огрешности коэффициентов матрицы А и правых частей b или же погрешности их округления сильно искажают р ешение системы. Для оценки обусловленности системы вводят число обусловленности M А МА А1 А . на. Чем больше значение M А ,тем система хуже обусловле- Свойства числа обусловленности: 1) М Е =1; 2) M А 1; 3) M А max | / | min , где м ах , min – соответственно максимальное и минимальное собственные числа матрицы А; 4) M АВ M А * M В ; 5) Число обусловленности матрицы А не меняется при умножении матрицы на произвольное число 0. Найдем выражение для полной оценки погрешности решения системы. Пусть в системе (2.1) возмущены коэффициенты матр ицы А и правая часть b, т.е. ~ ~ δ A A A , δ b b b , δ x ~x x . Теорема. Пусть матрица A имеет обратную матрицу, и выполняется условие A A 1 1 ~ . Тогда матрица А δ A A имеет обратную и справедлива следующая оценка относ ительной погрешности: δx x MA 1 MA δA δA δb A b . A 22 х1 а) х1 б) х2 х1 х2 в) х2 Рисунок 2 – а) система имеет единственное решение; б) система не имеет решения; в) система плохо обусловлена. В случае в) малейшее возмущение системы сильно м еняет положение точки пересечения прямых. В качестве примера рассмотрим систему 1 .03 x 1 0 .991 x 2 2 .51 0 .991 x 1 0 .943 x 2 2 .41 . Решение этой системы x 1 1.981 x 2 0.4735. Оценим влияние погрешности правых частей на результат. Рассмотрим ―возмущенную‖ систему с правой частью b * = (2.505 , 2.415) и решим эту систему: x1* x2 * 2.877 -0.4629. Относительная погрешность правой части (b) = 0.005/2.51 0.28% привела к относительной п огрешности решения (x * ) =0.9364/1.981 47.3%. Погрешность возросла примерно в 237 раз. Число об условленности системы (2.1) приблизительно равно 237. Подобные системы называются плохо обусловленн ыми. Возникает вопрос: какими методами можно решать такие системы? 23 2.1 Метод регуляризации для решения плохо обусловленных систем Рассмотрим систему Ах=b. (2.1) Для краткости перепишем эту систему в эквивален тный форме ( Ax b , Ax b ) 0 . (2.2) Для примера рассмотрим систему 2 x1 x1 1 . x1 2 x2 2 Тогда ее можно представить как (2х 1 -х 2 -1) 2 +(х 1 -2х 2 -2) 2 =0. (2.2*) Решение системы (2.2) совпадает с решением систем ы (2.2*). Если коэффициенты A или b известны неточно, то решение также является не точным, поэтому вместо равенс тва ( Ax b , Ax b ) 0 можем потребовать приближенного в ыполнения равенства ( Ax b , Ax b ) 0 и в этом виде задача становится не определенной и нужно добавить дополн ительные условия. В качестве дополнительного условия вводят требов ание, чтобы решение как можно меньше отклонялось от з аданного х 0 т.е. (х-х 0 , х-х 0 ) было минимальным. Следовательно, приходим к регуляризованной задаче вида (Ах-b, Ax-b)+ (x-x 0 , x-x 0 )=min, (2.3) где >0. Используя свойства скалярного произведения, выраж ение (2.3) перепишем в виде (x,A T Ax)-2(x,A T b)+(b,b)+ [(x,x)-2(x,x 0 )+(x 0 ,x 0 )]=min. (2.4) Варьируя x в уравнении (2.4), получим уравнение вида (A T А+ E)x=A T b+ x 0 . (2.5) 24 Система (2.5) – система линейных алгебраических уравнений, эквивалентная системе (2.1). Систему (2.5) р ешаем с помощью метода Гаусса или с помощью метода квадратных корней. Решая систему (2.5) найд ем решение, которое зависит от числа . Выбор управляющего параметра . Если =0, то система (2.5) перейдет в плохо обусловленную систему (2.1). Если же – велико, то система (2.5) переходит в х орошо обусловленную систему и решение этой системы может сильно отличаться от решения системы (2.1). Оптимальное значение – это такое число, при кот ором система (2.5) удовлетворительно обусловл ена. На практике пользуются невязкой вида r Ax α b , и эту невязку сравнивают по норме с известной погрешнос тью правых частей δ b и с влиянием погрешности коэффицие нтов матрицы δ А . δ b или δ А . Если – Если – слишком велико, то r δ b или δ А . мало, то r Поэтому проводят серию расчетов, при различных и в качестве оптимального значения выбирают то знач ение , когда выполнено следующее условие r δb δ A x . Для выбора вектора х 0 нужно знать приближенное решение или же, если приближенное решение трудно определить, то х 0 =0. 2.2 Метод вращения (Гивенса) Метод Гивенса, как и метод Гаусса состоит из прям ого и обратного ходов. Прямой ход метода. Исключаем неизвестное х 1 из всех уравнений, кроме первого. Для исключения х 1 из 2-го уравнения вычисляют числа α 12 a11 2 a11 2 a 21 , β 12 a 21 2 a11 2 a 21 , 25 2 2 и такие, что α 12 β 12 1 , β 12 а11 α 12 а 21 0 . Первое уравнение системы заменяем линейной комб инацией первого и второго уравнений с коэффициент ами α 12 и β 12 , а второе уравнение такой же комбинацией с α 12 и – β 12 . В результате получим систему где a11(1) x1 a12(1) x 2 ... a1(1m) x m b1(1) (1) a 22 x 2 ... a 2(1m) x m b2(1) a 31 x1 a 32 x 2 ... a 3 m x m b3 . (2.6) . . . . . . . . . a m1 x1 a m 2 x 2 ... a mm x m bm Здесь a1(1j) α12 a1 j b2(1) α12b2 β12 a2 j , a2(1)j α12 a2 j b1(1) β12 a1 j , α12b1 β12b2, β12b1, где j= 1, m . Преобразование системы (2.1) к системе (2.6) эквив алентно умножению слева матрицы A и вектора b на матрицу С 1 2 вида 12 12 C 12 12 ... 12 ... 1 ... 0 . . . . . . ... 1 Аналогично для исключения x1 из третьего уравнения вычисляем числа a11(1) 13 (1) 2 11 (a ) (1) 2 31 (a ) и a 31 13 (1) 2 11 (a ) (1) 2 31 , (a ) 2 2 (1) 0. такие, что 13 13 1, 13 a 31 13 a11 Затем первое уравнение системы (2.6) заменяем лине йной комбинацией первого и третьего уравнений с коэфф ициентами 13 , 13 , а третье уравнение системы (2.6) заменяем линейной комбинацией тех же уравнений, но с коэффициентами 13 и – 13 . Это преобразование эквивалентно умножению слева на матрицу 26 13 1 13 C 13 13 13 ... ... ... 1 ... . . . . . . ... 1 . Исключая неизвестное х 1 из всех последующих уравнений получим систему А ( 1 ) х=b (1 ) , где матрица на первом шаге A (1) =C 1 m …C 13 C 12 A, , а вектор правых частей b (1) C1m ...C13 C12 b . Здесь и далее через С k j обозначена матрица элемента рного преобразования, отличающаяся от ед иничной матрицы Е только четырьмя элементами. Действие матрицы С k j на вектор x эквивалентно повороту вектора x вокруг оси, перпендикулярной плоскости OX k X j на угол φ k j такой, что sin φkj . kj cos φkj , kj Операцию умножения на матрицу С k j называют плоским вращением или преобразованием Гиве нса. Первый этап состоит из m-1 шагов, в результате чего получается система (m a11 . 1) (m a12 x1 . . . 1) x2... a1(mm (1) a22 x2 ... . am(1)2 x2 . . ... 1) b1(m xm a2(1m) xm . . (1) amm xm 1) b2(1) . (2.7) . bm(1). В матричной форме получаем А ( 1) х=b (1) . На втором этапе, состоящем из m-2 шагов, из уравнений системы (2.7) с номерами 3,4,…,m исключают неизвестное х 2 . B результате получим систему 27 (m a11 . 1) x1 . . (m a12 1) x2 (m a13 1) x3... a1(mm 1) xm b1(m 1) (m a22 1) x2 (m a23 1) x3... a2(mm 1) xm b2(m 1) (2) a33 x3 ... . . . . . am(23) x3 В . ... . . (2) amm xm bm(2) . . bm(2). A(2) x b (2) , получаем где (2) (1) (2) (1) A C 2 m ...C 24 C 23 A , b C 2 m ...C 24 C 23 b . После завершения (m-1)-го шага придем к системе с верхней треугольной матрицей вида A ( m 1) x b ( m 1) , где матричной . a3(2m) xm A(m 1) Cm 1,m форме A(m 2) , b(m 1) Cm 1, m b (m 2) . Обратный ход метода вращений проводится точно так же, как и для метода Гаусса. 28 3 Решение нелинейных уравнений и систем нелинейных уравнений Рассмотрим систему нелинейных уравнений с m неизвестными вида f1 ( x1 ,..., x m ) 0 f 2 ( x1 ,..., x m ) 0 . (3.1) . . . . . f m ( x1 ,..., x m ) 0 Задача решения такой системы является более сло жной, чем нахождение корней одного нелинейного уравн ения, и чем задача решения линейных алгебраических ура внений. В отличие от систем линейных уравнений здесь использ ование прямых методов исключено и реш ение находится с использованием итерационных методов, т.е. находится приближенное решение x * = (x 1 * , ... , x m * ), удовлетворяющее при заданном > 0 условию х х . Задача (3.1) совсем может не иметь решения или же число решений может быть произвольным. Введем векто рную запись решения задачи: x=(x 1 ,…,x m ) T , f=(f 1 ,…,f m ) T , f(x)=0. ( 3.2) Будем считать, что функции f i непрерывно дифференцируемы в некоторой окрестности точки х. Введем матрицу Якоби f ( x) f1 f1 x1 f2 x2 f2 x1 . fm x2 . fm x1 x2 ... ... . ... f1 xm f2 xm . fm . xm 29 Как и в случае решения одного уравнения начинаем с этапа локализации решения (отделения корней). Пример. Дана система 2-х уравнений с двумя неизвестными х13 х23 х1 ln x2 8 х1х2 . x2 ln x1 Найдем на плоскости место расположения решения. Строим графики уравнений этой системы: а) – график 1-го уравнения, б) – график 2-го уравнения, в) – совмещенные графики. а) б) в) Рисунок 3 – Графики уравнений системы Определяем границы координат пересечения граф иков. Данная система имеет три решения. Координаты точек (B,C,A): B: x 1 =4, x 2 =4 C: 3,5 < x 1 < 4; 1,5 < x 2 < 2,5. 30 Точки А и С симметричны относительно прямой х 1 =х 2 . Координаты точки С определим приближенно: x 1 3,8, x 2 2. Обусловленность и корректность решения системы (3.1). Предположим что система (3.1) имеет р ешение х и в некоторой окрестности этого решения матрица Якоби не вырождена. Это означает, что в указанной окрестности нет других решений системы. В одномерном случае нахождение корня нелинейного уравнения приводит к определению интервала неопред еленности (х * – , х * + ). Так как значения функции f(x) чаще всего вычисляются на ЭВМ с использованием приб лиженных методов нельзя ожидать, что в окрестности корня относительная погрешность окажется малой. Сама погрешность корня ведет себя кра йне нерегулярно и в первом приближении может восприниматься как некоторая сл учайная величина. На рисунке 4а представлена идеальная ситуация, отвечающая исходной математической пост ановке задачи, а на рисунке 4б – реальная, соответствующая вычислениям значений функции f на ЭВМ. а) б) y y x* а x* b х x х*– х*+ Рисунок 4 – Графическое изображение определения интервала неопределенности В этом случае мы не можем определить, какая же то чка в интервале неопределѐнности является решением. Радиус интервала неопределенности прямо пропорционален погрешности вычисления значения f . Кроме того, возрастает (обусловленность задачи ухудшается) с уменьшением 31 f ( x* ) . Оценить величину довольно сложно, но выпол- нить это необходимо по следующим причинам: не имеет смысла ставить задачу о вычислении корня с точностью ε < ; после попадания очередного приближения в интервал неопределенности или близко от него, вычисления сл едует прекратить (этот момент для итерационных методов определяется крайне нерегулярным поведением приближений). Если случай многомерный, то получаем некоторую о бласть неопределѐнности D, и можем получить оценку радиуса этой области: ( f ( x )) f x 1 f ( x* ) (f) ( f ). Роль абсолютного числа обусловленности играет но рма матрицы, обратной матрице Якоби f ( x ) . Чем число обусловленности больше, тем хуже эта система обусловлена. 3.1 Метод простых итераций Систему (3.1) преобразуем к следующему эквивален тному виду: x1 φ1 ( x1 ,...,x m ) x 2 φ2 ( x1 ,...,x m ) . . . . . . x m φm ( x1 ,...,x m ) . (3.3) Или в векторной форме х=φ(х) ( 3.4) Пусть задано начальное приближение x ( 0 ) ( x1( 0 ) ..., xm( 0 ) )T . Подставляем его в правую часть системы (3.4) и получ аем x ( 1) = (x (0) ), продолжая подстановку, находим х ( 2 ) и т.д. 32 Получим последовательность точек { х ( 0 ) , х (1) ,..., х ( к 1) } , которая приближается к искомому решению х. 3.1.1 Условия сходимости метода. Пусть '(x) – матрица Якоби, соответствующая сист еме (3.4) и в некоторой -окрестности решения х функции φi (х ) (i=1,2,…,m) дифференцируемы и выполнено неравенство вида: φ ( х) q , где (0 q < 1), q − постоянная. Тогда независимо от выбора х (0 ) из -окрестности корня итерационная последовательность {х k } не выходит за пределы данной окрестности, метод сходится со скоростью геометрической прогрессии и справедлива оценка погре шности x (n) x q n x (0) x . 3.1.2 Оценка погрешности. В данной окрестности решения системы, производные функции i (x) (i=1,…,m) должны быть достаточно малы по абсолютной величине. Таким образом, если нераве нство φ ( х ) q не выполнено, то исходную систему (3.1) следует преобразовать к виду (3.3). Пример. Рассмотрим предыдущий пример и приведем систему к удобному для итераций виду x1 x2 3 8 x1x2 x2 x23 , x2 ln x2 x1 . ln x1 Проверяем условие сходимости вблизи точки С. Вычи слим матрицу Якоби 33 8 x1 3 x 22 8 x2 3(8 x1 x 2 x 23 ) 2 1 1 φ ( x1 , x 2 ) 3 ln 2 x1 ln x1 3(8 x1 x 2 x 23 ) 2 3 . 1 1 1 ln x 2 ln 2 x 2 Так как x 1 3,8, x 2 2, то при этих значениях вычисляем норму матрицы φ ( x ) || φ ( x ) || || φ (3,8 ; 2) || 0,815. Запишем итерационную процед уру ( k 1) x1 3 ( k 1) x2 (k ) x2 (k ) (k ) 8 x1 x 2 (k ) x2 (k ) ln x 2 (k ) ( x2 )3 , (k ) x1 (k ) . ln x1 Следовательно, метод простых итераций будет сходит ься со скоростью геометрической прогрессии, знам енатель которой q 0,815. Вычисления поместим в таблице 1. Таблица 1 Решение системы нелинейных уравн ений k 1 … 8 9 3,8000 3,75155 …. 3,77440 x 1 =3,77418 x1(k ) 2,0000 x2(k ) 2,03895 … 2,07732 x 2 =2,07712 При k=9 критерий окончания счета выполняется при =10 -3 и можно положить x 1 =3,774 0,001 и x 2 =2,077 0,001. 3.2 Метод Ньютона Суть метода состоит в том, что система н елинейных уравнений сводится к решению систем линейных алгебра ических уравнений. Пусть дана система (3.1) и задано начальное приближение x (0 ) . Приближение к решению х строим в виде последовательности x ( 0) , x (1 ) , …, x (n ) . 34 В исходной системе (3.1) каждую фун кцию fi ( x1, x2,..., xn ), где i= 1, m , раскладывают в ряд Тейлора в точке х (n ) и заменяют линейной частью еѐ разложения fi ( x ) fi ( x (n) m fi ( x (n) ) j 1 xj ) ( x j x (jn ) ) . В результате получим систему линейных алгебраич еских уравнений f1 ( x (n) m f1 ( x ( n ) ) j 1 xj ) ( x j x (jn ) ) 0 . ........... fm ( x (n) m ) j 1 fm ( x (n) ) xj x (jn ) ) (x j (3.5) В матричной форме f(x (n) )+ f'(x (n ) )*(x−x ( n) )=0, ( 3.6) где f ' − матрица Якоби. Предположим, что матрица не вырожденная, то есть су1 ществует обратная матрица f ( x ( n ) ) . Тогда система (3.6) имеет единственное решение, которое и принимается за очередное приближение x (n+1) . Отсюда выражаем решение x ( n+1 ) по итерационной формуле: x (n 1) x (n) f ' ( x (n) ) 1 f ( x (n) ) . (3.7) Формула (3.7) и есть итерационная формула метода Ньютона для приближенного решения системы нелин ейных уравнений. Замечание. В таком виде формула (3.7) используется редко в виду того, что на каждой итерации нужно нах одить обратную матрицу. Поэтому поступают следующим обр азом: вместо системы (3.6) решают эквиваленту ей систему линейных алгебраических уравнений вида 35 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) функции f i (при i= 1, m ) дважды непрерывно дифференцируемы и матрица Якоби не вырождена. Тогда найдется такая малая окрестность вокруг решения х, что при выборе начального приближения x 0 из этой окрестности итерационный метод (3.7) не выйдет за пределы этой окрестности решения и справедлива оценка вида x (n 1) x 1 x (n) x 2 , где n − номер итерации. Метод Ньютона сходится с квадратичной скор остью. На практике используется следующий критерий остановки: x (n) x (n 1) . 36 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 , и таким образом они б удут иметь единичную длину), нахождение собственных 37 значений матрицы A и соответствующих им собственных векторов и есть полное решение проблемы собственных значений. А нахождение отдельных собственных знач ений и соответствующих им векторов − называется решением частной проблемы собственных значений. Эта проблема имеет самостоятельное значение на пра ктике. Например, в электрических и механических с истемах собственные значения отвечают собственным частотам к олебаний, а собственные вектора характеризуют соответс твующие формы колебаний. Эта задача легко решается для некоторых видов ма триц: диагональных, треугольных и трехдиагональных матриц. К примеру, определитель треугольной или диагонал ьной матрицы равен произведению диагональных элеме нтов, тогда и собственные числа равны диагональным эл ементам. а 0 0 Пример. Матрица А – диагональная А 0 а 0 . 0 0 а Тогда det(А− Е)= (а λ)3 , а характеристическое уравнение (а λ)3 0 имеет трехкратный корень =а. Собственными векторами для матрицы А будут единичные векторы e1 1 0 , e2 1 , e3 0 . 1 Пример. Найдем собственные числа матрицы А 2 9 5 1, 2 5,999 6 1 1 . 7 ,5 Составим характеристический многочлен 38 2 Р3 ( ) det ( A 3 10 ,8999 E ) det 2 1, 2 1 26 , 49945 9 5,3999 1 21,002 . 5 6 7 ,5 Используя метод Ньютона, определим один из корней уравнения Р 3 ( )=0, а именно -7,87279. Разделив многочлен P 3 (λ) на ( − 1 ) получим многочлен второй степени: P 2 (λ)= 2 +3,02711 +2,66765. Решив квадратное уравнение, находим оставшиеся два корня: (комплексные сопряженные 2 ,3 -1,51356 0,613841*i корни). Существуют прямые методы нахождения собственных значений и итерационные методы. Прямые методы неудо бны для нахождения собственных значений для матриц в ысокого порядка. В таких случаях с учетом возможностей компьютера более удобны итерационные методы. 4.1 Прямые методы нахождения собственных значений 4.1.1 Метод Леверрье Метод разделяется на две стадии: - раскрытие характеристического уравнения, - нахождение корней многочлена. Пусть det(A- E) – есть характеристический многочлен матрицы А={a ij } (i,j=1,2,…,m), m p1 m 1  p m , и 1 , 2 , …, m – есть полная т.е. det A E совокупность корней этого многочлена (полный спектр собственных значений). Рассмотрим суммы вида Sk λk1 λk2  λkm (k=1,2,…,m), т.е. 39 S1 λ1 λ 2 ... λ m SpA S 2 λ12 λ 22  λ 2m SpA 2 , ......... S m λ1m λ m2  λ mm SpA m где SpA m ( 4.3) a ii − след матрицы. i 1 В этом случае при k m справедливы формулы Ньют она для всех (1 k m) Sk p1Sk 1  pk 1S1 kpk , Откуда получаем при k=1 р 1 = -S 1 , при k=2 р 2 = -1/2∙ (S 2 + р 1 ∙S 1 ), . . . . . . . . . . . . . . при k=m р m = -1/n∙ (S m + р 1 ∙S m-1 + р 2 ∙S m-2 +. . .+ + р m- 1 ∙S 1 ). ( 4.4) ( 4.5) Следовательно, коэффициенты характеристического многочлена р i можно определить, если известны суммы S 1 ,S 2 , ... ,S m . Тогда схема алгоритма раскрытия характер истического определителя методом Леверрье будет следу ющей: 1) вычисляют степень матрицы: А k =А k- 1 ∙А для k=1,…,m; 2) определяют S k − суммы элементов, стоящих на главной диагонали матриц А k ; 3) по формулам (4.5) находят коэффициенты характер истического уравнения р i (i=1,2,…,m). 4.1.2 Усовершенствованный метод Фадеева Алгоритм метода: 1) вычисляют элементы матриц A 1 ,A 2 ,.., A m : 40 A1 A; SpA1 q1 ; B1 A1 q1 E ; SpA 2 A 2 A B1 ; q 2 ; B 2 A2 q 2 E ; 2 ................................. SpA m Am A B m 1 ; q m ; B m Am q m E , m (в конце подсчета B m – нулевая матрица для контроля); 2) определяют коэффициенты характеристического уравнения р i : q 1 = –р 1 , q 2 = –р 2 ,..., q m = –р 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 имеют одни и те же собственные числа. Это можно легко увидеть ра скрыв определитель det ( S 1 AS E ) det ( S 1 ( A det ( S 1 ) det ( A E)S) E ) det ( S ) det ( A E) . Следовательно, матрицы A и B − подобные, имеют одни и те же собственные значения. Н о собственные векторы х и 41 не совпадают, они связаны между собой простым соо тношением х = S∙х'. Такую матрицу A с помощью преобразования подобия или же последовательности таких преобразов аний можно привести к матрице Фробениуса вида: x f11 f12  f1,m 1  1        1 F 1 f1m . Детерминант матрицы F – det(F) можно разложить по элементам первой строки: det ( F E ) ( 1) m ( m p1 m 1  p m ) . Тогда коэффициенты характеристического многочл ена матрицы А будут р 1 = f 11 , p 2 = f 12 , …, p m = f 1m . . Второй случай. Матрицу А преобразованием подобия можно привести к матрице В верхнего треугольного вида B b11 b12  b1m b22  b2 m . . . .  bmm . Тогда собственными числами будут диагональные эл ементы матрицы B: det ( B E ) (b11 )( b22 )(bmm ). Третий случай. Матрицу A с помощью преобразования подобия можно привести к Жордановой форме S 1 AS Λ 42 Λ λ1 S1  λ2 S2         λm , где i – собственные числа матрицы A; S i – константы (0 или 1); если S i =1, то i = i+1 . К четвѐртому случаю относятся матрицы, которые с п омощью преобразования подобия можно привести к диагональному виду (матрица простой структуры):  2   ...  ...  n 1 S 1 AS D , у которой, как известно, собственными числами я вляются диагональные элементы. 43 4.1.4 Метод итераций определения первого собственн ого числа матрицы. Пусть дано характеристическое уравнение: det(A− ∙E) = 0, где − собственные значения матрицы А. Предположим, что | 1 |>| 2 | | 3 | … | m |, т.е. 1 – наибольшее по модулю собственное число. Тогда для нахождения приближенного значения 1 используется следующая схема: 1) выбирают произвольно начальный вектор у (0 ) ; 2) строят последовательность итераций вида: 1, 2 ,..., n y (1) Ay ( 0 ) , y ( 2 ) А Аy ( 0 ) A 2 y ( 0 ) , ........ y ( m ) А А m 1 y (0) A m y (0) , y ( m 1) A A m y ( 0 ) A m 1 y ( 0 ) . 3) выбирают y ( m ) A m y ( 0 ) и y ( m 1 lim n y i( m 1) 1) y i( m ) или A m 1 y ( 0 ) , тогда y i( m 1 1) y i( m ) , где y i( m ) и y i( m 1 ) – соответствующие координаты векторов y ( m) и y (m+1) . Возникает вопрос выбора начального вектора у ( 0 ) . При неудачном выборе можем не получить значения нужного корня, или же предела может не существ овать. Этот факт при вычислении можно заметить по прыгающим значениям этого отношения, следовательно, нужно изменить у ( 0 ) . В качестве первого собственного вектора можно взять вектор у (n +1 ) и пронормировать его. Пример. Найти наибольшее по модулю собственное значение и соответствующий ему собственн ый вектор матрицы А 44 A 3 1 1 2 2 . 1 1 1 1) Выбираем начальный вектор y (0) 1 . 1 2) Вычисляем последовательно векторы y ( 1) , y ( 2) , …, y (10 ) . Вычисления помещаем в таблицу 2. Таблица 2 – Вычисление векторов y (0 ) А∙y ( 0) А 2 ∙y (0 ) А 3 ∙y (0) …….. 1 4 17 69 1 5 18 67 1 2 7 25 y (n + 1 ) А 9 ∙y (0) 243569 210663 73845 А 1 0 ∙y (0 ) 941370 812585 284508 3) Вычисляем отношения координат векторов yi(10) и yi(9) (1) 1 y1(10 ) y1( 9 ) 4)Вычисляем 3,865 ; 1 y 2(10 ) (2) 1 y 2( 9 ) 3,857 ; (3) 1 y 3(10 ) y 3( 9 ) 3,853 . как среднее арифметическое λ(11), λ(12), λ(13) . λ1 λ(11) λ(12) 3 λ(13) 3,858 . 5) Определим соответствующий числу вектор: 1 собственный 941370 y (10 ) A (10 ) y ( 0 ) 812585 . 284508 Нормируем вектор y y (10 ) 3 941370 ( 10) 2 , разделив на его длину 812585 2 284508 2 1, 28 10 6 получим вектор 45 0,74 x1 0,64 . 0, 22 Далее можем определить второе с обственное число λ2 yi(n 1) yi(n 1) λ1 yi(n) λ1 yi(n 1) , где i=1,2,…,n. При вычислении собственных чисел подобным обр азом, будет накапливаться ошибка. Данная методика позволяет приближенно оценить собственные значения ма трицы. 46 5 Задача приближения функции Постановка задачи. Пусть на отрезке [a,b] функция у=f(x) задана таблицей своих значений y0 f ( x0), ..., yn f ( xn ) . Допустим, что вид функции f(x) неизвестен. На практике часто встречается задача вычисления значений фун кции у=f(x) в точках х, отличных от x 0 ,..., x n . Кроме того, в некоторых случаях, не смотря на то, что аналитическое в ыражение у=f(x) известно, оно может быть слишком гр омоздким и неудобным для математических прео бразований (например, специальные функции). Кроме эт ого значения y i могут содержать ошибки эксперимента. Определение . Точки x0, ..., xn называются узлами интерполяции. Требуется найти аналитическое выражение функции F(x), совпадающей в узлах интерполяции со значениями данной функции, т.е. F ( x0) y0, F ( x1) y1, ..., F ( xn ) yn. Определение. Процесс вычисления значений функции F(x) в точках отличных от узлов интерполирования называется интерполированием функции f(x). Если x x0, xn , то задача вычисления приближе нного значения функции в т. х называется интерполированием, иначе – экстраполированием. Геометрически задача интерполирования функции одной переменной означает построение кривой, проход ящей через заданные точки ( x0, y0), ( x1, y1),..., ( xn, yn ) (рисунок 5). То есть задача в такой постановке может иметь беск онечное число решений. 47 у x 0 x1 x2 x3 х Рисунок 5 − Геометрическая иллюстрация задачи интерполирования функции Задача становится однозначной, если в качестве F(x) выбрать многочлен степени не выше n, такой что: F n (x 0 )=y 0 , F n (x 1 )=y 1 , ..., F n (x n )=y n . Определение. Многочлен F n (x), отвечающий вышеназванным условиям, называется интерполяционным мног очленом. Знание свойств функции f позволяет осознанно выбирать класс G аппроксимирующих функций. Широко и спользуется класс функций вида Фm ( x) c0φ0( x) c1φ1( x) ... cmφm ( x), (5.1) являющихся линейными комбинациями некоторых бази сных функций 0 (x), ..., m (x). Будем искать приближающую функцию в виде многочлена степени m, с коэффициентами с 0 , ..., с m , которые находятся в зависимости от вида приближения. Функцию Ф m (х) называют обобщенным многочленом по системе функций 0 (х), 1 (х), …, m (х), а число m – его степенью. Назовем обобщенный многочлен Ф m (х) интерполяционным, если он удовлетворяет условию Ф m (х i )=y i , (i=0,1,…,n). (5.2) Покажем, что условие (5.2) позволяет найти прибл ижающую функцию единственным образом 48 c0 c0 ( x 0 ) c1 1 ( x 0 ) ... c m 0 ( x1 ) c1 1 ( x1 ) ... c m . . . . . . . . . . c 0 0 ( x n ) c1 1 ( x n ) ... c m ( x0 ) y 0 y1 m ( x1 ) m . yn, m ( xn ) (5.3) Система (5.3) есть система линейных алге браических уравнений относительно коэффициентов с 0 ,с 1 ,…,с m . Эта система n линейных уравнений имеет единстве нное решение, если выполняется условие m=n и определитель квадратной матрицы Р 0 ( x0 ) , 1( x 0 ) ,..., 0 ( x1 ) , 1( x1 ) n ( x0 ) ,..., n ( x1 ) det P 0. . . . 0 ( xn ) , . . . 1( x n ) ,..., . . n ( xn ) Определение. Система функций 0 (x),..., n (x) называется Чебышевской системой функций на [a,b],если определитель матрицы отличен от нуля detP≠0 при любом расположении узлов x i [a,b], i=0,1,…n, когда среди этих узлов нет совпадающих. Если мы имеем такую систему функций, то можно утверждать, что существует единственный для данной сист емы функций интерполяционный многочлен Ф m (х), коэффициенты которого определяются единственным обр азом из системы (5.3). Пример. При m n система функций 1, х, х 2 ,…, х m линейно независима в точках х 0 , х 1 , …, х n , если они попарно различны. 5.1 Интерполяционный многочлен Лагранжа Рассмотрим случай, когда узлы интерполирования не равноотстоят друг от друга на отрезке [a,b]. 49 Тогда шаг h=x i+1 −x i ≠const. Задача имеет единственное решение, если в качестве интерполирующей функции F(x) взять алгебраический многочлен L n (x)=a 0 +a 1 x+a 2 x 2 +…+a n x n , где а i неизвестные постоянные коэффициенты. Используя условие (5.2) можем записать (5.4) Ln ( x0 ) y 0 , Ln ( x1 ) y1 ,..., Ln ( x n ) y n . Запишем это в виде: a 0 a1 x0 a12 x02 ... a1n x 0n a 0 a1 x1 a12 x12 ... a1n x1n y0 y1 (5.5) . . . . . . . . . . a 0 a1 x n a12 x n2 ... a1n x nn y n . Эта система однозначно разрешима, так как система функций 1,х,х 2 ,…,х n линейно независима в точках х 0 ,х 1 ,…,х n . Однозначная разрешимость следует из того фа кта, что определитель этой системы (определитель Ванде рмонда) 1 x0 x02 ... x0n 1 x1 x12 ... x1n ... ... ... ... ... 1 xn x n2 ... x nn 0 j i n ( xi x j ) 0 . Без вывода приведем одну из форм записи интерпол яционного многочлена Лагранжа Ln ( x ) y 0 y1 ( x x1 ) ...( x x n ) ( x 0 x1 )...( x 0 x n ) ( x x 0 )( x x 2 )...( x x n ) ( x1 x 0 )( x1 x 2 )...( x1 x n ) ... y n ( x x 0 )...( x x n 1 ) ( x n x 0 )...( x n x n 1 ) (5.6) . 50 Определение. Этот многочлен называется интерпол яционным многочленом Лагранжа и сокращенно записыв ается в виде n L n (x)= yi i 0 ( x x0 )( x x1 )...( x xi 1 )( x xi 1 )...( x x n ) ( xi x0 )( xi x1 )...( xi xi 1 )( xi xi 1 )...( xi x n ) (5.7) . На практике часто пользуются линейной и квадрати чной интерполяцией. В этом случае формула Лагранжа им еет вид L 1 (x)= y 0 L 2 (x)= y 0 ( x x0 ) ( x0 x1 ) y1 ( x x0 ) − при линейной интерполяции; ( x1 x0 ) ( x x1 )( x x 2 ) ( x0 x1 )( x0 x 2 ) y1 ( x x0 )( x x 2 ) ( x1 x0 )( x1 x 2 ) y2 ( x x0 )( x x1 ) ( x 2 x0 )( x 2 x1 ) − при квадра- тичной интерполяции. Рассмотрим теперь случай с равноотстоящими узлами. Тогда интерполяционная формула Лагра нжа заметно упрощается. В этом случае шаг h=x i+1 -x i =const. Введем в рассмотрение многочлен вида Qi ( x ) ( x x0 )( x x1 )...( x xi 1 )...( x x n ) ( xi x0 )( xi x1 )...( xi xi 1 )...( xi x n ) x x0 Введем обозначение q= x x0 = q h , x x1 = q h h . h . , отсюда следует, что h (q 1) , . . . . . . . . . . x xi = q h i h h (q i), x x n = q h n h h (q n). Тогда многочлен Q i примет вид q ( q 1) [ q (i 1)] [ q (i 1)]...( q n ) h n Qi ( x ) . i h (i 1) h... h ( h )...[ ( n i ) h ] Произведя простейшие преобразования, получим выр ажение вида: 51 Cni q(q 1)...(q n) = ( 1) , q i n! n! где C ni – число сочетаний из n элементов по i C ni . i !( n i )! q (q 1)...( q n) ( 1)n Q i (q)= (q i) i!(n i)! i n i Тогда интерполяционный многочлен Лагранжа для ра вноотстоящих узлов имеет вид: Ln ( x) q(q 1)...( q n! n) n ( 1) i 0 n i Cni y. q i i 5.1.1 Оценка погрешности интерполяционного мног очлена Оценить погрешность интерполяционной формулы Л агранжа можно только тогда, когда известно аналитич еское выражение интерполируемой функции, а точнее, если и звестно максимальное значение ( n+1)-ой производной функции f(x) на отрезке [a,b]. Пусть |R n (x)| =| f(x) −L n (x)|, где R n (x) –погрешность; f(x) − точное значение функции в точке х; L n (x) − приближенное значение, полученное по полин ому Лагранжа. Если обозначить через M n 1 = f ( n 1) (ξ ) max f ( n 1) ( x ) , где x [ a ,b ] ξ [a,b], причем х 0 =а, х n = b, то Rn ( x ) Mn 1 ( n 1)! (ξ x0 )( ξ x1 ) ....( ξ xn ) . 5.2 Интерполяционные полиномы Ньютона 5.2.1 Интерполяционный многочлен Ньютона для ра вноотстоящих узлов Вычисление значений функции для значений аргуме нта, лежащих в начале таблицы удобно проводить, пользуясь первой интерполяционной формулой Ньютона. Для этого введем понятие конечной разности. 52 Определение. Конечной разностью перового порядка называется разность между значениями функции в сосе дних узлах интерполяции. Тогда конечные разности в то чках х 0 ,х 1 ,…,х n -1 Δ y0 y1 y0 f ( x1) f ( x0) Δ f ( x0) , Δ y1 . y2 . Δ yn y1 . 1 f ( x2) . yn yn Δ f ( x1) , f ( x1) . . . . 1 f ( xn ) . . f ( xn 1) . Δ f ( xn 1) . Конечная разность второго порядка имеет вид: 2 yi . yi . . yi , 1 . . . n yi ( n 1 y i ). Рассмотрим некоторые свойства конечных разностей. Вторая конечная разность в точке х i Δ 2 y i f ( xi 1 Δ x ) f ( xi Δ x ) f ( xi 2 ) 2 f ( xi 1 ) f ( xi ) y i f ( xi 1 ) f ( xi ) 2 yi 1 yi 2 . Аналогично третья конечная разность 3 yi yi 3 3 yi 2 3 yi 1 yi . Общее выражение для конечной разности n-го порядка имеет вид Δn yi  yn ( 1) m i Cn1 yn i 1 Cnm yn i m Cn2 yn i 2 ... n ...( 1) yi , а вообще, конечная разность порядка m от конечной разности порядка n Δ m (Δ n y ) Δ m n y . Конечные разности n-го порядка от многочлена степени n – есть величины постоянные, а конечные разности n+1-го порядка равны нулю. 53 Для вычисления значений функции в начале таблицы требуется построить интерполяционный многочлен степени n такой, что выполнены условия интерполяции Pn ( x0 ) y 0 ,..., Pn ( x n ) y n . В силу единственности многочлена степени n, построенного по n+1 значениям функции f(x) многочлен Pn ( x ) , в конечном счете, совпадает с многочленом Лагранжа. На йдем этот многочлен в виде: Pn ( x) a0 a1 ( x x0) a2( x x0)( x x1) ... an ( x x0)...( x xn 1) , где а i (i=0,1,…, n) – неизвестные коэффициенты. Для нах ождения а 0 положим x x0 . Тогда P ( x0 ) a 0 , отсюда а 0 =у 0 . Для вычисления a1 рассмотрим первую конечную разность для многочлена P n (x) в точке х. Pn ( x ) Pn ( x h ) Pn ( x ) [ a 0 a1 ( x x0 h ) ... a n ( x x0 h ) ... ...( x x n 1 h ) ] a 0 a1 ( x x0 ) ... a n ( x x0 )...( x x n 1 ) . В результате преобразований получим Pn ( x ) h a1 2 ha 2 ( x x0 ) ... n ha n ( x x0 )...( x x n 1 ). Вычислим первую конечную разность многочлена P n (x) в точке х 0 Pn ( x0 ) a1 h , но Pn ( x0 ) f ( x1 ) f ( x0 ) y1 y 0 y0 , откуда a1 Δ y0 . h Чтобы определить коэффициент а 2 , составим конечную Pn ( x h ) Pn ( x ). Отразность второго порядка 2 Pn ( x ) сюда после преобразования получим a 2 2 y0 2!h 2 . Вычисляя конечные разности более высоких порядков и п олагая х=х 0 , придем к общей формуле для определения коэффицие нi тов: ai y0 i !h i (i=0,1,2,…,n). 54 Подставим значения a i в многочлен, в результате пол учим первую интерполяционную формулу Ньют она: Pn ( x ) y 0 y0 1!h n ( x x0 ) ... y0 n !h n ( x x0 )...( x x n 1 ). Первую интерполяционную формулу можно зап исать в том виде, в котором ее удобнее использовать для интерп олирования в начале таблицы. Для этого введем пер еменную q=(x-x 0 )/h, где h– шаг интерполирования. Тогда первая формула примет вид Pn ( x ) y 0 q y 0 q ( q 1) 2 2! y 0 ... q ( q 1)...( q n 1) n! n y0 . 5.2.2 Вторая интерполяционная формула Ньютона Эта формула используется для интерполирования в конце таблицы. Построим интерполяционный многочлен вида Pn ( x ) a 0 a 1( x xn ) a 2 ( x xn )( x xn 1 ) ... a n ( x xn )...( x x1 ). Неизвестные коэффициенты а 0 ,а 1 ,…,а n подберем так, чтобы были выполнены равенства Pn ( x0 ) y 0 , Pn ( x1 ) y1 ,..., Pn ( x n ) y n . Для этого необходимо и достаточно, чтобы i Pn ( x n i ) i y n i (i=0,1,…,n). В случае, если положить x=x n , то сразу определяется коэффициент а 0 Pn ( x ) yn a 0 . Из выражения для первой конечной разности на йдем a1 : Δ Pn ( x) 1 ha 1 2 ha 2( x xn 1) ... n ha n( x Отсюда, полагая х=х n- 1 , получим a1 x n 1)( x Δ yn h 1 x n 2)...( x x1). . Из выражения 55 для второй конечной разности найдем а 2 : a2 Δ2 yn 2! h щая формула для коэффициента а i имеет вид ai 2 2 . Об- Δi yn i! h i i . Подставим эти коэффициенты в формулу многочлена и получим вторую интерполяционную формулу Ньютона: yn Pn ( x ) y n h n 1 ( x x n ) ... y0 n !h n ( x x n )...( x x1 ). На практике используют формулу Ньютона в другом виде. Положим q=(x-x n )/h. Тогда Pn ( x ) y n q y n q ( q 1) 1 2! 2 yn 2 ... q ( q 1)...( q n 1) n! n y0 . 5.3 Интерполирование сплайнами Многочлен Лагранжа или Ньютона на всем отре зке [a,b] с использованием большого числа узлов интерпол ирования часто приводит к плохому приближению, что объясняется накоплением погрешностей в ходе вычислений. Кроме т ого, из-за расходимости процесса интерполирования увел ичение числа узлов не обязательно приводит к повышению точности вычислений. Поэтому построим такой вид приближения, который: позволяет получить функцию, совпадающую с табли чной функцией в узлах; приближающая функция в узлах таблицы имеет непрерывную производную до нужного п орядка; В силу вышесказанного на практике весь отрезок [a,b] разбивается на частичные интервалы и на каждом из них приближающая функция f(x) заменяется многочленом невысокой степени. Такая интерполяция называется кусочнополиномиальной интерполяцией. Определение. Сплайн − функцией называют кусочнополиномиальную функцию, определенную на отрезке [a,b] и имеющую на этом отрезке некоторое число н епрерывных производных. 56 Слово сплайн означает гибкую линейку, которую и спользуют для проведения гладких кривых через определенное число точек на плоскости. Преимущество спла йнов – сходимость и устойчивость процесса вычисления. Рассмо трим частный случай (часто используемый на практике), к огда сплайн определяется многочленом третьей ст епени. 5.3.1 Построение кубического сплайна Пусть на отрезке [a,b] в узлах сетки заданы значения a x 0 x1 x 2 ... x n b , некоторой функции f(x), т.е. yi f ( xi ) (i= 0,1,…, n). Сплайном, соответствующим этим узлам функции f ( x ) называется функция S(х), которая: 1) на каждом частичном отрезке является многочл еном третьей степени; 2) функция S(x) и ее две первые производные S ( x), S ( x) непрерывны на [a,b]; 3) S ( xi ) f ( xi ) . На каждом частичном отрезке [ xi 1, xi ] будем искать сплайн S ( x ) Si ( x ) , где S i (x ) многочлен третьей степени Si ( x) ai bi ( x xi ) ci 2 (x xi ) 2 di 6 (x xi ) 3 . (5.8) То есть для x [ xi 1, xi ] нужно построить такую функцию Si (x) , где a i , bi , ci , d i подлежат определению. Для всего отрезка интерполирования [a,b], таким образом, необходимо определить 4 n неизвестных коэффициента. S ( x ) bi ci ( x xi ) di ( x xi ) 2 , 2 S ( x ) ci d i ( x xi ), S i ( x ) ai yi . Доопределим a 0 f ( x0 ) y0 . Требование непрерывности функции S(x) приводит к условия S i ( xi ) S i 1 ( xi ), (i=0,1,…,n-1). Отсюда из (5.8) получаем следующие уравнения: 57 di 1 ( xi xi 1)3 (i= 1,2,…,n-1). 1 6 Введем шаг интерполирования hi x i x i 1 . ai ai bi 1( xi xi 1) ci 1 (x 2 i xi 1)2 Тогда последнее равенство можно переписать в виде hi bi hi2 2 hi3 ci 6 di f i 1 (i= 1,2,…,n). fi Из непрерывности первой производной следует hi ci hi2 2 di bi bi 1 (i=2,3,…,n), а из непрерывности второй производной h i d i =c i −c i−1 (i=2,3,…,n). Объединив все три вида уравнений, получим си стему из 3n-2 уравнений относительно 3n неизвестных bi, ci, di . Два недостающих уравнения получим, задав граничные условия для функции S(x). Для этого воспользуемся граничными условиями для сплайн-функции в виде S (a) S (b) 0 (концы гибкой линейки свободны). Тогда получим систему уравнений hi d i hi ci hi bi ci ci 1 , c0 c n 0 , (i 1, 2 ,..., n ) hi2 d i bi bi 1 , (i 2 ,3,..., n ) 2 hi2 hi3 ci d i f i f i 1 , ( i 1, 2 ,..., n ). 2 6 (5.9) Решая систему методом подстановки (исключаем из (5.9) неизвестные b i ,d i ), получим систему: h i ci 1 2 ( h i h i 1) c i h i 1c i 1 6( yi yi 1 hi 1 y i yi hi 1 ), (i 1, 2,..., n ) (5.10) c0 c n 0 Система (5.10) имеет трехдиагональную матрицу. Эта система может быть решена методом прогонки или Гаусса. Метод прогонки рассматривается в пункте 6.7.2. 58 После решения системы коэффициенты сплайна di , bi определим через коэффициенты с i с помощью явных формул ci di ci 1 , hi hi bi hi2 ci 2 6 yi di yi (i= 1,2,…,n). 1 hi Существуют специальные виды записи сплайнов на к аждом из промежутков [x i ,x i+1 ] [9], которые позволяют уменьшить число неизвестных коэффициентов сплайна. Вводятся обозначения S ( xi ) hi mi , i xi и xi 1 0 ,1,..., n , t (x xi )/ hi . На отрезке [x i ,x i+1 ] кубический сплайн записывается в виде S ( x) yi (1 t ) 2 (1 2t ) yi 1t 2 m i h i t (1 t ) 2 (3 2t ) m i 1t 2 (1 t ) h i . Кубический сплайн, записанный в таком виде, на каждом из промежутков [x i ,x i+1 ] непрерывен вместе со своей пе рвой производной на [a,b]. Выберем m i таким образом, чтобы и вторая производная была непрерывна во всех внутренних узлах. Отсюда пол учим систему уравнений: λ imi где μ i 2 m i μ i mi 1 hi hi 1 1 hi , λi 1 1 3( μi μi yi hi 1 yi hi hi 1 hi λi yi , i yi hi 1 ), 1 1,2 ,..., n 1 . К этим уравнениям добавим уравнения, полученные из граничных условий 2 m0 m1 3 y1 y0 h0 , mn 1 2mn 3 yn yn hn 1 . 1 В результате получаем аналогичную систему с трехдиагональной матрицей. Решаем систему линейных уравнений относительно коэффициентов m i методом пргонки. 59 5.3.2 Сходимость процесса интерполирования кубич ескими сплайнами Доказывается, что при неограниченном увеличении числа узлов на одном и том же отрезке [a,b] S ( x) f ( x) . Оценка погрешности интерполяции R( x) f ( x) S ( x) зависит от выбора сетки и степени гладкости фун кции f(x). При равномерной сетке xi a i h (i=0,1,…,n) f(x) Sh( x ) где M 4 M 4 h4 8 max | f IV ( x) | . , [a, b ] Другие постановки задачи интерполирования фун кций. 1. Если функция периодическая, то используется тр игонометрическая интерполяция с периодом l, которая строится с помощью тригонометрического мн огочлена n Tn ( x ) a0 ( a k cos k 1 kx l bk sin kx l ), коэффициенты которого находятся из системы уравнений Tn ( xi ) f ( xi ) (i= 1,2,…, 2n+1). 2. Выделяют приближение функций рациональными, дробно – рациональными и другими функциями. В да нном пособии эти вопросы не рассматриваются. 5.4 Аппроксимация функций методом наименьших квадратов К такой задаче приходят при статистической обрабо тке экспериментальных данных с помощью регрессионного анализа. Пусть в результате исследования некоторой величины x значениям x1 ,x 2 ,x 3 ,...,x n поставлены в соответствие значения y1 ,y 2 ,y 3 ,...,y n некоторой величины у. Требуется подобрать вид аппроксимирующей зависим ости y=f(x), связывающей переменные х и у. Здесь могут иметь место следующие случаи. Во -первых: значения функции f(x) могут быть заданы в достаточно большом ко60 личестве узлов; во-вторых: значения таблично заданной функции отягощены погрешностями. Тогда проводить пр иближения функции с помощью интерполяционного многочлена нецелесообразно, т.к. - число узлов велико и пришлось бы строить нескол ько интерполяционных многочленов; - построив интерполяционные многочлены, мы повтор или бы те же самые ошибки, которые присущи та блице. Будем искать приближающую функцию из следующих соображений: 1) приближающая функция не проходит через узлы та блицы и не повторяет ошибки табличной функции; 2) чтобы сумма квадратов отклонений приближающей функции от таблично заданной в узлах таблицы была минимальной. у уn yn-1 y1 отклонения y0 х0 х1 … хn-1 хn х Рисунок 6 – Графическое изображение отклон ений приближающей функции от таблично заданной Рассмотрим линейную задачу наименьших квадр атов. Определение. Уровень погрешности, допускаемый при снятии характеристики измеряемой величины, наз ывается шумом таблицы. Пусть функция y=f(x) задана таблицей приближенных значений y i f(x i ), i=0,1,…,n, полученных с ошибками y i0 y i , где yi0 f ( xi ) . i Пусть даны функции 0 ( x ), 1( x ),..., m ( x ) , назовем их базисными функциями. Будем искать приближающую (аппроксимирующую) функцию в виде линейной комбинации базисных функций 61 Фm ( x ) y c0 0 ( x) c1 1 ( x) ... cm m ( x) . (5.11) Такая аппроксимация называется линейной, а Ф m (х) – обобщенным многочленом. Будем определять коэффициенты обобщенного многочлена c 0 ,…,c m используя критерий метода наименьших квадратов. Согласно этому критерию вычислим сумму квадратов отклонений таблично заданной функции от и скомого многочлена в узлах: n ( y i Ф m ( x i )) m 2 i 0 n ( y i c 0 0 ( x i ) ... c m m ( x i )) 2 . (5.12) i 0 Выражение для m можно рассматривать как функцию от неизвестных c 0 ,…, c m . Нас интересует, при каких значениях c 0 ,…, c m , значение m будет минимально. Для этого воспользуемся условием существования эк стремума функции, а именно, найдем частные прои зводные от m по всем переменным c 0 ,…, c m и приравняем их к нулю. Получим систему вида: ∂ m ∂ c0 n 2 ( yi c0 0( xi ) ... c m m ( x i )) 0 ( xi ) i 0 . . . . . . . . . . . . . . . . . . ∂ m ∂cm (5.13) n 2 ( yi c0 0 ( xi ) ... c m m ( x i )) m ( xi ) i 0 Система (5.13) – система линейных уравнений относ ительно c 0 ,…, c m . Чтобы систему (5.13) записать компактно, ведем определение. Определение. Скалярным произведением функций f на g на множестве точек x 0 ,..., x n называется выражение n f ( xi ) g ( xi ) . ( f ,g) i 0 Тогда систему (5.13) можно записать в виде: 62 c0 ( 0, 0) c1 ( 0 , 1) c0 ( 1, 0) c1 ( 1 , 1) . c0 ( . . m , . 0) . c1 ( m ... c m ( ... c m ( . . , 1) . . ... c m ( 0, 1, m) m) . . m , ( ( 0 , y) 1 , y) . m) . ( . m . . (5.13а) , y) Системы (5.13) или (5.13а) будем называть нормальной системой уравнений. Решив ее, мы найдем коэффициенты c 0 ,…,c m и следовательно, найдем вид аппроксимирующего многочлена. Н апомним, что это возможно, если базисные функции линейно независимы, а все узлы различны. Осталось определить степень многочлена m. Прямому вычислению поддаются только значения среднеквадрати чного отклонения m , анализируя которые будем выб ирать степень многочлена. Алгоритм выбора степени многочлена m. В случае, когда m=n мы получим интерполяционный многочлен, поэтому выберем m<0 и 2 >0 должны быть такими, чтобы m находилось между ними; 2) первоначально m выбирают произвольно, но учит ывая условие, что m< 1 , то степень необходимо уменьшить х отя бы на единицу; б) если m < 2 , то степень необходимо увеличить хотя бы на единицу. 5) затем строим приближающую функцию . 63 Очень часто для приближения по методу наименьших квадратов используются алгебраические многочлены степени m n, т.е. k ( x ) x k . Тогда нормальная система (5.13) принимает следующий вид: m n ( n x ij k ) c j j 0 i 0 y i x ik , (k= 0,1,…,m). (5.14) i 0 Запишем систему (5.14) в развернутом виде в двух на иболее простых случаях m =1 и m =2. В случае многочлена первой степени P 1 (x)=c 0 +c 1 x, нормальная система имеет вид n ( n 1) c 0 n ( x i ) c1 i 0 n n ( x i ) c0 i 0 2 i ( yi i 0 (5.15) n x ) c1 i 0 yi x i . i 0 P 2 (x)=c 0 +c 1 x+c 2 x 2 , Для многочлена второй степени нормальная система имеет вид n (n 1)c0 i 0 n ( i 0 xi )c1 i 0 n n ( ( xi )c0 xi2 )c0 n ( xi2 )c1 i 0 n ( i 0 xi3)c1 ( xi2 )c2 i 0 n ( xi3)c2 i 0 n ( i 0 xi4 )c2 n yi i 0 n yi xi . i 0 n (5.16) yi xi2 i 0 64 6 Численные методы решения задачи Коши для обыкновенных дифференциальных уравнений и систем дифференциальных уравнений Будем рассматривать задачу Коши для системы обы кновенных дифференциальных уравнений(ОДУ).Запишем систему в векторной форме du dt f (t, u ) , (6.1) где: u – искомая вектор-функция; t– независимая переменная; u (t ) (u1 (t ) ,..., u m (t )) ; f (t, u ) ( f 1 ,..., f m ) , m– порядок системы; u1 (t ) ,..., u m (t ) координаты; t 0; u (0) u 0 . Запишем систему (6.1) в развернутом виде du i dt f i (t , u1 ,..., u m ) , ( 6.2) где: i=1,...,m; ui (0) ui0 . В случае i=1 – это будет ОДУ 1-го порядка, а при i=2 – система из двух уравнений первого порядка. В случае i=1 решение задачи Коши предполагает н ахождение интегральной кривой, проходящей через зада нную точку и удовлетворяющую заданному начальном у условию. Задача состоит в том, чтобы найти искомую ве кторфункцию u, удовлетворяющую (6.1) и заданным начал ьным условиям. Известны условия, гарантирующие существование и единственность решения (6.1) или (6.2). Предположим, что функции f i (i=1, …, m) непрерывны по всем аргументам в некоторой замкнутой области D={t a , u i u i0 b }, где a,b– известные константы. Из непрерывности функций следует их ограниче нность, т.е. функции fi сверху ограничены некоторой константой М: | f i |0. По оси t введем равномерную сетку с шагом , т.е. t n n , n 0 , 1, 2 , .. . . Оборассмотрим систему точек ω τ значим через u(t) точное решение (6.4) , а через yn y(tn ) приближенные значения функций u в заданной системе точек. Приближенное решение является сеточной функцией, т.е. определено только в точках сетки ω τ . 6.1 Семейство одношаговых методов решения задачи Коши 6.1.1 Метод Эйлера (частный случай метода Рунге-Кутта) Уравнение (6.4) заменяется разностным уравнением yn 1 τ yn f ( t n , y n ) , n=0,1,2,…,. y0 u0 66 В окончательной форме значения y n 1 можно определить по явной формуле yn yn 1 τ f (t n , y n ) . (6.5) Вследствие систематического накопления ошибок м етод используется редко или используется только для оценки вида интегральной кривой. Определение 1. Метод сходится к точному решению в некоторой точке t , если yn u(tn ) 0, при , tn t . Метод сходится на интервале (0, t], если он сходится в каждой точке этого интервала. Определение 2. Метод имеет р-й порядок точности, если существует такое число р>0, для которого y n u (t n ) O ( p ) , при , где: – шаг интегрирования; O – малая величина порядка τ p . Так как u n 1 τ un u (t n ) O ( τ ) , то метод Эйлера имеет первый порядок точности. Порядок точности разностного метода совпадает с порядком аппроксимации исходного дифференциального уравнения. 6.1.2 Методы Рунге-Кутта Метод Рунге-Кутта второго порядка точности Отличительная особенность методов Рунге -Кутта от метода (6.5) заключается в том, что значение правой части уравнения вычисляется не только в точках сетки, но и та кже в середине отрезков(промежуточных точках). Предположим, что приближенное значение y n решения задачи в точке t=t n уже известно. Для нахождения y n+1 поступают следующим образом: 1) используют схему Эйлера в таком виде yn 1 2 0 ,5 τ yn f (t n , y n ) ( 6.6) 67 и отсюда вычисляют промежуточное значение y n+ 1/2 ; 2) воспользуемся разностным уравнением вида yn 1 yn f (t n τ 0 ,5 τ , y n 12 ), ( 6.7) откуда найдем значение y n +1 . Далее подставим значение y n 1 2 y n 0 , 5 τ f n в уравнение (6.7). Тогда получим разностное уравнение yn 1 yn τ 0 ,5 τ, y n f (t n 0 ,5 τ f n ) , (6.8) где fn f (tn, yn ) . Можно показать, что метод (6.8) имеет второй пор ядок точности, т.е. y n u (t ) O(τ 2 ) . Реализация метода (6.8) в виде двух этапов (6.6),(6.7) называется методом предиктор–корректор (прогноза и коррекции) в том смысле, что на первом этапе (6.6) приближенное значение предсказывается с невысокой точностью O (τ ) , а на втором этапе (6.7) это предсказанное значение исправляется и результирующая погрешность им еет второй порядок O ( τ 2 ) . Будем рассматривать явные методы. Задаем числовые коэффициенты a i , b ij , i=2,...,m; j=1,2,...,m-1 и σ i , i=1,2,...,m. Последовательно вычисляем функции k1 f (t n ,y n ) ; k2 f (t n a 2 τ, y n b21 τ k1 ) ; k 3 f (t n a3 τ, y n b31 τ k1 b32 τ k 2 ) ; …………………………………………….. kn f (t n a n τ, y n bm1 τ k1 ... bmm 1 τ k m 1 ) . Затем из формулы yn yn m i k i находим новое знаτ i 1 чение y n 1 =y(t n +1 ). Здесь ai, bij , σ i –числовые параметры, ко1 торые определяются или выбираются из соображений то чности вычислений. Чтобы это уравнение аппроксимировало уравнение (6.4), необходимо потребовать выполнения ус68 m ловия 1 . При m=1 получается метод Эйлера, при m=2 i i 1 получаем семейство методов yn yn 1 τ( 1k1 2k2 ) , (6.9) где: k1 f (t n , y n ) ; k 2 f (t n a 2 τ , y n b21 τ k1 ) ; y 0 =u 0 . Семейство определяет явные методы Рунге -Кутта. Подставив нужные 1 и 2 , получаем окончательную формулу. Точность этих методов совпадает с точностью аппроксим ирующего метода и равна O ( τ 2 ) . Невязкой, или погрешностью аппроксимации метода (6.9) называется величина un n un 1 f (t n , u n ) 1 τ 2 a 2 τ ,u n f (t n b21 τ f (t n , u n )) , полученная заменой в (6.9) приближенного решения то чным решением. При 1 + 2 =1 получим первый порядок точности. Если же потребовать дополнительно σ 2b21 σ 2a2 0,5 , то получим методы второго порядка точности вида yn 1 yn (1 ) f (t n , y n ) f (t n a , yn a f (t n , y n )) при σ a 0,5 . Приведем один из методов Рунге -Кутта третьего порядка точности yn yn 1 1 τ где: k 1 f ( t n , y n ) ; k3 f (t n 6 k2 f (t n , yn k1 ( k1 k3 ) , 4k 2 τ , yn k1 ) ; 2 2 2 k2 ) . Метод Рунге-Кутта 4-го порядка точности yn 1 yn τ где: k 1 k3 f ( t n , y n ); f (t n τ 2 , yn 6 f (tn k2 τ 2 1 k 2); ( k 1 2 k 2 2 k 3 k 4) , τ 2 , yn k4 τ 2 k 1) ; f (tn τ, y n τ k 3) . 69 Методы Рунге-Кутта сходятся и порядок их точности совпадает с порядком аппроксимации разностным отнош ением. Теорема. Пусть правая часть уравнения (6.4) f (t, u ) удовлетворяет условию Липшица по аргументу u с константой L, и пусть n – невязка метода Рунге-Кутта. Тогда для погрешности метода при n τ T справедлива оценка y n u ( t n ) Te T max δn , где: σ L m (1 Lb τ ) m 1 ; max ; b i max bij . i, j i На практике обычно пользуются правилом Рунге. Для этого сначала проводят вычисления с шагом τ , затем – τ/2. τ2 Если y n – решение при шаге τ, а y2n – при шаге τ 2 , p – порядок точности метода, то справедлива оценка. y 2 n2 y 2 n2 u (t n ) 2p yn 1 Тогда за оценку погрешности для метода четвертого порядка точности при шаге τ 2 принимают величину max y nτ y 2τ n2 15 i . 6.2 Многошаговые разностные методы решения з адачи Коши для обыкновенных дифференциальных уравнений Рассмотрим задачу Коши для обыкновенного дифф еренциального уравнения du dt f (t, u) , (6.10) где u (0) u 0 . Для решения задачи Коши для уравнения ( 6.10) при t>0 введем равномерную сетку с постоянным шагом ωτ t n n τ, n 0 ,1,... . 70 Введем понятие линейного m шагового разностного метода для решения задачи (6.10). Линейным m-шаговым разностным методом называется система разностных уравн ений a0 y n a1 y n 1 ... a m y n m b0 f n τ b1 f n ... bm f n 1 m , (6.11) где: n=m,m+1...; a k , b k – числовые коэффициенты, не зависящие от n; k=0,1,…,m, причем a 0 ≠0. Систему (6.11) будем рассматривать как рекуррент ные соотношения, выражающие новые значения yn y(tn ) через ранее найденные значения yn 1, yn 2,..., yn m , причем расчет начинают с индекса n=m, т.е. с уравнения a0 y m a1 y m 1 ... a m y 0 b0 f m τ b1 f m 1 ... bm f 0 . Отсюда следует, что для начала расчета по формулам (6.11) надо знать m предыдущих значений функции y, причем y 0 =u 0 определяется исходной задачей (6.10). Эти предыдущие m значений могут быть найдены одним из одн ошаговых методов Рунге-Кутта. Отличие от одношаговых методов состоит в том, что по формулам (6.11) расчет ведется только в то чках сетки. Определение. Метод (6.11) называется явным, если коэффициент b 0 =0. Тогда значение y n легко выражается через yn 1, yn 2,..., yn m . В противном случае метод назыв ается неявным, и для нахождения y n придется решать нелинейное уравнение вида a0 τ ak m yn b0 f (t n ,y n ) ( bk t n k 1 k τ yn k ). (6.12) Обычно это уравнение решают методом Н ьютона при (0) y n 1 . Коэффициенты уравнения начальном значении y n (6.11) определены с точностью до множителя, тогда, чт обы устранить этот произвол, вводят условие m bk 1 , с тем ус- k 0 ловием, что правая часть (6.11) аппроксимирует п равую часть дифференциального уравнения (6.10). 71 На практике используют частный случай методов (6.11), т.н. методы Адамса, когда производная u (t ) аппроксимируется разностным отношением, включающим две с оa1 1 ; a k =0, k=2,...,m и седние точки t n и t n 1 . Тогда a 0 yn yn m 1 bk f n τ k . (6.13) k 0 Это и есть методы Адамса. При b 0 =0 метод будет явным, в противном случае – неявным. 6.2.1 Задача подбора числовых коэффици ентов a к , b к Выясним, как влияют коэффициенты a k , b k на погрешность аппроксимации уравнения (6.11), на устойчивость и сходимость. Определение. Невязкой, или погрешностью аппроксимации методов (6.11) называется функция m m ak rn un bk f ( t n k k 0 k ,un k ), (6.14) k 0 которая получается в результате подстановки точного решения u(t) дифференциального уравнения (6.10) в разностное уравнение (6.11) . Если разложить функции u n k u (t n k ) в ряд Тейлора в точках t t n равномерной сетки, окончательно получим функцию m rn ( k 0 ak τ p )u (t n ) m ( ( kτ ) l 1 l 1 k 0 (ak k l bk ) ) u ( l ) (t n ) (l 1)! O ( τ p ) . (6.15) Из вида функции rn следует, что порядок аппроксим ации будет равен p, если выполнены условия m ak 0; k 0 m k l 1 ( ka k lb k ) k 0 0; (6.16) m bk 1, k 0 где l=1,...,p. 72 Условия (6.16) представляют собой систему линейных уравнений из p+2 уравнений относительно неизвестных a 0 ,…,a m , b 0 ,…,b m . Их количество равно 2(m+1). Решив систему (6.16), получаем неизвестные числовые коэффицие нты. Для неявных m-шаговых методов наивысшим достижимым порядком аппроксимации будет p=2m, а для явных – p=2m-1. Запишем систему (6.16) для методов Адамса m l k l 1b k 1; k 1 (6.17) m b0 1 bk , k 1 где l=2,...,p . Отсюда видно, что наивысший порядок аппроксимации для неявного m-шагового метода Адамса p=m+1, для явного (b 0 =0) – p=m. 6.2.2 Устойчивость и сходимость многошаговых разн остных методов Наряду с системами уравнений (6.11) будем рассмат ривать т.н. однородные разностные уравнения вида a 0V n a1V n 1 ... a m V n m 0 , (6.18) где n=m,m+1,... . Будем искать его решение в виде функ ции Vn q n , где q – число, подлежащее определению. Подставив Vn в (6.18) получаем уравнение для нахождения q a0 q m a1 q m 1 ... a m 1q am 0. (6.19) Уравнение (6.19) принято называть характеристич еским уравнением для разностных методов (6.11). Говорят, что разностный метод (6.11) удовлетворяет условию корней, если все корни уравнения (6.19) q 1 , …, q m лежат внутри или на границе единичного круга комплексной плоскости, пр ичем на границе нет кратных корней. Разностный метод (6.11), удовлетворяющий условию корней, называется устойчивым методом. 73 Теорема. Пусть разностный метод (6.11) удовлетворяf (t , u ) ет условию корней и выполнено условие t T . Тогда при m малых T , n m и достаточно tn n будет выполнена оценка yn u (t n ) M ( max 0 j m 1 yj L при u u (t j ) max 0 k n m rk ) , (6.20) где: r k − погрешность аппроксимации; y j u (t j ) − погрешность в задании начальных условий; M –константа , зависящая от L,T и не зависящая от n. Методы Адамса удовлетворяют условию корней, т.к. для них a 0 =-a 1 =1, следовательно, q=q 1 =1. 6.2.3 Примеры m-шаговых разностных методов Ада мса Явные методы. При m=1 порядок точности p=1. Тогда метод описывается формулой yn yn 1 fn 1. В этом случае получаем метод Эйл ера. При m=2 порядок точности p=2. Тогда метод описывается формулой yn yn 3 1 2 fn 1 1 2 fn 2. При m=3 порядок точности p=3. Тогда метод описывается формулой yn yn 1 τ 1 12 ( 23 f n 1 16 f n 2 5 fn 3 ) . При m=4 порядок точности p=4. Метод описывается формулой yn yn 1 1 24 (55 f n 1 59 f n 2 37 f n 3 9 fn 4) . Неявные m-шаговые методы Адамса: m=1, p=2, m=2, p=3, yn yn 1 yn yn 1 1 ( fn f n 1 ) – метод трапеций; 2 1 (5 f n 8 f n 1 f n 2 ) ; 12 74 yn m=3, p=4, yn 1 1 24 (9 f n 19 f n 1 5 fn 2 fn 3). Неявные методы содержат искомое значение y n нелинейным образом, поэтому для его нахождения п рименяют итерационные методы решения нелинейных уравн ений. 6.3 Численное интегрирование жестких обыкновенных дифференциальных уравнений систем Жесткие системы можно сравнить с плохо обусловле нными системами алгебраических уравнений. Рассмотрим систему дифференциальных уравнений(ДУ) (6 d u f (t, u ) , .21) dt u (0) u0 . Для решения (6.21) рассмотрим разн остные где методы вида 1 τ m m ak yn k 0 k b k f( t n k k 0 , yn k ), (6.22) где n= m, m+1, m+2,…. Устойчивость и сходимость методов (6.22) определяе тся расположением корней характеристического уравн ения, т.е. |q| 1 – корни принадлежат единичному кругу. Среди методов (6.22) выделим те, которые позволяют получить асимптотически устойчивые решения. Пример. В качестве частного случая (6.21) рассмотрим уравнение вида du u, dt <0; u(t ) u0e λt – решение ДУ. При (6.23) где: u(0) u0 ; <0 решение есть монотонно убывающая функция при t . Для этого решения можно записать при любом шаге >0 u (t τ) u (t ) , (6 .24) что означает устойчивость решения u(t). 75 Рассмотрим для задачи (6.23) метод Эйлера yn 1 yn yn , τ где: n=0,1,2,…, yn 1 q yn , q – промежуточный параметр, равный 1+ . Оценка (6.24) для метода Эйлера будет выполнена т огда и только тогда, когда |q| . Шаг лежит в интервале 0 < ≤ . Метод Эйлера для задачи (6.23) устойчив при в ыполнении этого условия. Определение 1. Разностный метод (6.22) называется абсолютно устойчивым, если он устойчив при любом >0. Определение 2. Разностный метод называется условно устойчивым, если он устойчив при некоторых ограничен иях на шаг . Например, метод Эйлера для (6.23) условно устойчив, при 0 < ≤ . Примером абсолютно устойчивого мет ода является неявный метод Эйлера yn 1 yn yn 1 , τ т.к. в этом случае q (1 τ) 1 1 , при любых >0. Замечание. Условная устойчивость является недоста тком явных методов в связи с тем, что приходится выбирать мелкий шаг интегрирования. Пример для задачи (6.23). Если =−200, тогда 0.01. Если мы рассмотрим интервал (0,1], то необходимо б удет 100 шагов. Неявные методы со своей стороны прив одят к решению на каждом шаге нелинейного уравнения, но это уже недостаток неявного метода. 6.3.1 Понятие жесткой системы ОДУ Замечание. Все вышерассмотренные методы легко реализуются на примере одного уравнения и легко перенося тся на системы ДУ, но при решении систем возникают д ополнительные трудности, связанные с разнома сштабностью описанных процессов. Рассмотрим пример системы двух уравнений: 76 d u1 dt du 2 a1 u 1 , a2 u2 dt где: t >0; a 1 ,a 2 >0. Эта система однородных независимых ДУ имеет решение u1 ( t ) u1 ( 0 ) e u 2 (t ) u 2 (0) e a1 t a2 t . Это решение монотонно убывает с ростом t. Пусть коэффициент а 2 на порядок больше а 1 , т.е. а 2 >>a 1 . В этом случае компонента u 2 затухает гораздо быстрее u 1 , и тогда, начиная с некоторого момента времени t , решение задачи u(t) почти полностью будет определяться поведением ко мпоненты u 1 . Однако при численном решении данной задачи шаг интегрирования будет определяться, как правило, компонентой u 2 , не существенной с точки зрения поведения решения системы. Рассмотрим метод Эйлера для р ешения данной системы u1( n 1) u1( n ) τ u 2( n 1) u 2( n ) τ a1 u1( n ) . a 2 u 2( n ) Он будет устойчив, если на шаг наложены ограниче- ния τ a1 2 τ a2 2 . Учитывая, что a 2 >> a1 >0, получаем окончательное ограничение на τ 2 a2 . Такие трудности могут возникнуть при решении л юбых систем ОДУ. 77 Рассмотрим в качестве примера систему du A u, dt (6.25) где А– квадратная матрица m*m. Если матрица А имеет большой разброс собственных чисел, то возникают проблемы с разномасштабностью оп исываемых системой процессов. Допустим, что матрица А постоянна (т.е. не зависит от t). Тогда система (6.21) будет называться жесткой, е сли: 1) вещественные части собственны х чисел k 0 для всех k, где k=1,…,m; 2) число S max | Re 1 k m min | Re 1 k m k | k | велико (десятки и сотни), и число S называется числом жесткости системы. Если же матрица А зависит от t, то и собственные числа зависят от t и k зависят от t. Решение жесткой системы (6.25) содержит как быстро убывающие, так и медленно убывающие составляющие. Начиная с некоторого t >0 решение системы определяется медленно убывающей составляющей. При использов ании явных разностных методов быстро убывающая с оставляющая отрицательно влияет на устойчивость, поэтому приходится брать шаг интегрирования слишком мелким. 6.3.2 Некоторые сведения о других методах решения жестких систем На практике разностные методы (6.22) для решения ж естких систем используются в виде методов Гира (неявный разностный метод) и метода матричной экспоне нты (метод Ракитского). 6.3.2.1 Методы Гира b0 Это частный случай методов (6.22), когда коэффиц иент 1 , b1 b2 ... bm 0 . Запишем числовые коэффициенты, 78 которые определяются из условия p-го порядка точности аппроксимации системы разностными методами m m ak ; a0 m 1; k ak k 1 k 1 k l ak 0, k 1 (6.26) где l=2,...,p. Решив систему линейных уравнений (6.26) с учетом предыдущих условий, получаем все нужные коэффицие нты. Трехшаговый метод Гира (частный случай методов (6.22) с учетом условий (6.26)) имеет вид 11 6 yn 3 3 yn 1 2 1 yn 2 3 yn τ f (t n , y n ). 3 (6.27) Метод имеет третий порядок точности. При m=4, получаем четырехшаговый метод Гира 25 y n 48 y n 1 36 y n 2 16 y n 3 3 yn 4 12 τ (6.28) f (t n , y n ) . Запишем систему (6.26) в виде a1 a1 a1 2 a2 .... m am 1 2 2 2 a 2 .... m a m 0 . ...................................... 2 m a2 .... m m am (6.29) Решив (6.29) для каждого случая можем найти коэфф ициенты a k , k=1,2,…,т. Чисто неявные разностные методы обладают хорош ими свойствами устойчивости, поэтому используются для р ешения жестких систем уравнений. 79 6.3.2.2 Метод Ракитского (матричной экспоненты) решения систем ОДУ du (6 .30) A u, dt где: u (u1 ,..., u n ) ; u (0) u 0 ; А– матрица размерности n*n. Допустим, что матрица А – постоянная, т.е. ее элементы не зависят от времени. Система (6.30) – однородная, с постоянными коэффициентами. Запишем аналитическое р ешение (6.30) (6.31) u e A t u0 , где e A t – матричная экспонента и e At E A t ( A t)2 ... 2! ( A t)n n! +…. (6.32) Проинтегрируем уравнение (6.30) при значениях t= , 2 , 3 ,…. Если точно знать матрицу e A τ , то точное решение в указанных точках можно получить по формуле (6.31), т.е. решение можно записать u|t u|t 2τ eA τ u0; τ eA τ u x τ ; Таким образом, задача сводится к тому, чтобы достаточно точно знать матрицу e A τ . На практике поступают следующим образом: при больших нельзя воспользоваться рядом Тейлора в связи с его бесконечностью, т.е. для удовлетворительной точности пришлось бы взять мно го членов ряда, что трудно. Поэтому поступают так: отрезок [0, ] разбивают на k частей, чтобы длина h= /k удовлетворяла условию ||A∙h||<0.1. Тогда запишем по схеме Горнера e Ah E A h (E A h 2 (E A h 3 (E A h 4 ))) . 80 Каждый столбец матрицы e A h − w j вычисляют по формуле A h k w j (e ) w j , где w 0j – вектор столбец, в i-ой строке которого 1, а в остальных – нули. Если эта матрица найдена, то решение находится по (6.31). Для исследования разностных методов при решении ж естких систем рассматривают модельное ура внение du dt λ u, (6.33) где – произвольное комплексное число. Для того, чтобы уравнение (6.33) моделировало исхо дную систему (6.30) его нужно рассматривать при таких значениях , которые являются собственными числами матрицы А. Многошаговые разностные методы (6.31) имеют вид m μ bk ) y n (a k k (6.34) , k 0 где: n=m, m+1…; ∙ . n Если решение уравнения (6.34) искать в виде y n q , то для нахождения числа q получим характеристическое уравнение вида m (a k μ bk ) q m k 0. k 0 Для устойчивости метода достаточно выполнения усл овия корней | q k | 1 . В случае жестких систем испол ьзуются более узкие определения устойчивости. Предварительные сведения. Областью устойчивости разностных методов называется множество всех точек комплексной плоскости ∙ , для которых разностный метод применительно к уравнению (6.33) устойчив. 81 Определение 1. Разностный метод называется Аустойчивым, если область его устойчивости содержит л евую полуплоскость Re <0. Замечание. Решение модельного уравнения (6.33) асимптотически устойчиво при значениях Re <0, поэтому сущность А-устойчивого метода заключается в том, что Аустойчивый разностный метод является абсолютно усто йчивым, если устойчиво решение исходного дифференц иального уравнения. Так как класс А-устойчивых методов узок, то пользую тся А( )-устойчивым методом. Определение 2. Разностный метод (6.31) называется А( )-устойчивым, если область его устойчивости соде ржит угол меньший , т.е. |arg(- )|< , где ∙ . Исходя из этого, определяется, что при А( ) устойчивость совпадает с определением Аустойчивого метода. 6.4 Краевые задачи для обыкновенных дифференциальных уравнений Постановка краевой задачи. Рассматриваем дифференциальное уравнение порядка n (n 2) (6.35) y ( n ) f ( x , y , y ,..., y ( n 1) ) . Если сделаем замену переменных вида y0 y1; y1 y 2 ; .......... .......... .......... .. y n 1 f ( x , y1 ,..., y n 1 ); y k ( x0 ) y0( k ) , (6 .36) где: yk ( x0) y0(k ) ; k=0,1,…,(n-1), то задача (6.35) сводится к задаче Коши для нормальной системы ОДУ порядка n. Типовые примеры краевых задач. Рассмотрим дифференциальное уравнение 82 F ( x, y, y ,..., y (n)) 0. (6.37) Для уравнения (6.37) краевая задача формулируется следующим образом: найти решение y=y(x), удовлетворяющее уравнению (6.37), для которой значения ее прои зводных в заданной системе точек x xi удовлетворяют n независимым краевым условиям, в общем виде нелине йным. Эти краевые условия связывают значения искомой функции y и ее производных до (n-1) порядка на границах заданного отрезка. y B A a b x Рисунок 7 – Краевые условия для случая 1 ние y y 1. Рассмотрим уравнение второго порядка y f ( x , y , y ) . Необходимо найти решение уравнения, удовлетворяющее заданным краевым усл овиям: y(a)=A, y(b)=B, т.е. необходимо найти интегральную кривую, проходящую через две заданные точки (рисунок 7). 2. Рассмотрим уравнеf ( x , y , y ) с краевыми условиями y ( a ) A1 , y (b ) B1 . Из графика на рисунке 8 видно, что tg( )=A 1 , tg( )=B 1 . Здесь интегральная кривая пересекает прямые x=a и x=b под α β заданными углами и a b x соответственно. Рисунок 8 – Краевые условия для случая 2 83 3. Смешанная краевая задача включает оба случая. f ( x, y , y ) , Рассматривается то же самое уравнение y но с краевыми условиями y(a)=A 1 , y (b)=B 1 . Геометрическую иллюстрацию этих краевых условий легко представить, используя рисунки 7 и 8. Замечание. Краевая задача для уравнения (6.37) в общем случае может не иметь решений, иметь единс твенное решение, иметь несколько решений или бесконечное мн ожество решений. 4. Поражение заданной цели баллистическим снар ядом. Дифференциальные уравнения движения сн аряда с учетом сопротивления воздуха имеют вид E cos Θ x E sin Θ y g , где: x – вторая производная по времени; E=E(y,v) – известная функция высоты и скорости; v x 2 y 2 ; g=g(y) – ускорение силы тяжести; – угол наклона к горизонту касательной к траектории движения снаряда; Θ M(x1, y1) y1 x1 y x . Предполагая, что при t=t 0 снаряд выпущен из точки, совпадающей с началом координат с начальной скоростью v 0 под углом 0 , а в момент t=t 1 он поразит неподвижную мишень в точке M ( x1 , y1 ) получаем краевые условия y v0 arctg x Рисунок 9 – Траектория снаряда x 0, y x 0, x v 0 cosΘ0 , x1 , y y1 , y v 0 sinΘ0 , при t при t t0 , t1 . 84 Здесь неизвестны значения и t 1 . Решив данную крае- вую задачу, можем найти начальный угол Θ0 x 0 x (t 0 ); y 0 y (t 0 ) ; arctg y 0 x 0 , где: – угол, при котором поражается цель в точке M. Аналогично ставятся краевые задачи для систем дифф еренциальных уравнений. 6.5 Решение линейной краевой задачи Рассмотрим важный частный случай решения краевой задачи, когда дифференциальное уравнение и краевые у словия линейны. Для этого рассмотрим уравнение p 0 ( x) y (n) ( x) p1 ( x ) y ( n 1) ( x ) ... f ( x) , p n ( x) y ( x) (6.38) где: pi (x) и f(x) известные непрерывные функции на отре зке [a, b]. Предположим, что в краевые условия входят две аб сциссы x=a, x=b. Это двухточечные краевые задачи. Кра евые условия называются линейными, если они имеют вид n 1 R ν ( y) (ν) k ( y (k ) (ν) k (a ) k 0 y (k ) (b )) = , (6.39) где: − заданные константы. Причем они одновр еменно не равны нулю, т.е. n 1 [| k 0 ( ) k | ( ) k | 0 , при v=1,2,…,n. |] Например, краевые условия во всех трех рассмотре нных ранее задачах линейны, т.к. их можно записать в в иде причем 1, 1 0, 1 y (a ) y (b ) A, 1, 1 1 y (a ) γ1 1 y (b ) γ2 0, 2 , B – для первой задачи. 85 6.6 Решение двухточечной краевой задачи для линейного уравнения второго порядка сведением к задаче Коши Запишем линейное уравнение второго порядка в виде y p( x) y q( x) y f ( x) , (6.40) где: p, q, f − известные непрерывные функции на некот ором отрезке [a,b]. Требуется найти решение уравнения (6.40), удовлетв оряющее заданным краевым условиям y (a) 1 y (a) A, y (b ) 1 y (b ) B. Причем константы и (6.41) одновременно не равны нулю 1 0, 1 0. Решение задачи (6.40), (6.41) будем искать в виде л инейной комбинации y=C u+V, где С – константа, u – общее решение соответствующего однородного уравнения u p( x) u q( x) u 0, (6.42) а V– некоторое частное решение неоднородного уравн ения V p( x) V q( x) V f ( x) . (6.43) Потребуем ,чтобы первое краевое условие было выпо лнено при любом C, C ( 0 u (a) ( 0 V (a) A, 1 u ( a )) 1 V ( a )) 0, откуда следует, что 0 u ( a ) 1 u (a) A. 0 V (a) 1 V (a) Тогда u( a ) u (a) 1k 0k , (6.44) где k– некоторая константа, не равная нулю. Значение функции V и ее производная в точке а могут быть, например,выбраны равными 86 V (a ) если коэффициент α 0 A V (a ) V (a ) , (6.45) , (6.46) и V (a ) A 1 если коэффициент α 1 . Из этих рассуждений следует, что функция u – есть решение задачи Коши для однородного уравнения (6.42) с начальными условиями (6.44), а функция V – есть решение задачи Коши для неоднородного уравнения (6.43) с начал ьными условиями (6.45) или (6.46) в зависимости от усл овий. Константу C надо подобрать так, чтобы выполнялись условия (6.41) (вторая строчка) в точке x=b C ( 0 u (b ) B. 1 u (b )) 0 V (b ) 1 V (b ) Отсюда следует, что B C ( V (b ) u (b ) 1 1 V (b ) u (b ) , где знаменатель не должен быть равен нулю, т. е. u (b ) u (b ) 0. (6.47) Если условие (6.47) выполнено, то краевая задача (6.35), (6.36) имеет единственное решение. Если же (6.47) не выполняется, то краевая задача (6.35), (6.36) либо не имеет решения, либо имеет бесконечное множество реш ений. 1 6.7 Методы численного решения двухточечной краевой задачи для линейного уравнения второго порядка 6.7.1 Метод конечных разностей Рассмотрим линейное дифференциальное уравн ение y p( x) y q( x) y f ( x) (6.48) 87 с двухточечными краевыми условиями y (a ) 1y 0 y (b ) ( (a ) A 1 y (b ) B 1 0, 1 ; (6.49) 0) , где: p, q, f – известные непрерывные функции на некот ором отрезке [a,b]. Одним из наиболее простых методов решения этой краевой задачи является сведение ее к системе конечноразностных уравнений. Основной отрезок [a,b] делим на n равных частей с шагом h=(b-a)/n, то есть рассматриваем равноме рную сетку xi x0 i h , i=0,1,…,n. В каждом внутреннем узле дифференциальное уравн ение (6.48) аппроксимируем, используя формулы численн ого дифференцирования второго порядка точности yi yi 1 ; 2 h yi 2 2 yi yi y 1 i h yi (6.50) 1 2 , где i=1,..., n-1. Для граничных точек x0 a и xn b , чтобы не выходить за границы отрезка, используем формулы численного дифференцирования первого порядка точн ости y0 y1 y0 h , y yn yn 1 n h . (6.51) Используя отношение (6.50) исходное дифференциальное уравнение (6.48) аппроксимируем конечно разностными уравнениями yi 1 2 yi h 2 yi 1 pi yi yi 2 h 1 1 qi y i fi , (6.52) 88 где i=1,...,n-1. Учитывая краевые условия, получим еще два уравнения y1 y0 1 h yn yn y0 1 yn 1 h A . (6.53) B Таким образом, получена линейная система n+1 уравнений с n+1 неизвестными y0, y1,..., yn , представляющими собой значения искомой функции y y(x) в точках x0, x1,..., xn . Разностная схема (6.52)-(6.53) аппроксимируеткраевую задачу (6.48)- (6.49) с порядком 1 по h за счет краевых условий. Решив эту систему, получим таблицу значений искомой функции y. 6.7.2 Метод прогонки (одна из модификаций метода Гаусса) При применении метода конечных разностей к кра евым задачам для дифференциальных уравнений второго порядка получается система линейных алгебраических уравнений с трехдиагональной матрицей, т.е. каждое ура внение системы содержит три соседних неизвестных. Для решения т аких систем разработан специальный метод – «метод прогонки». Для этого систему (6.52) перепишем в виде yi 1 m i yi n i yi 1 ~ 2 fi h (6.54) для внутренних точек (i=1,…,n-1), где: 2 mi 1 h2 qi ; ni pi h 2 1 1 pi h 2 ; ~ fi pi h 2 1 fi . pi h 2 На концах отрезка x 0 =a и x n =b производные заменяем разностными отношениями 89 y y1 y0 и h y yn yn 1 n h . Учитывая эту замену, получим еще два уравнения y0 yn y1 1 y0 A; h yn 1 1 (6.55) yn B. h Обратим внимание на внешний вид записи системы (6.54), (6.55). В каждом уравнении системы пр исутствует три ненулевых элемента. В первом и последнем – по два ненулевых коэффициента. Разрешая уравнение (6.54) относительно y i , получим ~ fi yi mi h2 1 mi yi ni 1 mi yi 1 . (6.56) Предположим, что с помощью полной системы (6.54), (6.55) из уравнения (6.56) исключена неизвестная y i-1 . Тогда это уравнение примет вид y i ci (d i y i 1 ) , (6.57) где: c i , d i – некоторые коэффициенты; i=1,2,…,n-1. Отсюда y i 1 c i 1 (d i 1 y i ) . Подставляя это выражение в уравнение (6.54), пол учим ~ yi 1 mi yi nici 1(di 1 yi ) fi h 2 , а отсюда ~ ( fi h 2 yi nici 1di 1) mi nici 1 yi 1. (6.58) Сравнивая (6.57) и (6.58), получим для определ ения c i и d i рекуррентные формулы ci mi 1 nici 1 ; di ~ 2 fi h nici 1di 1; i=1,…,n-1. (6.59) Из первого краевого условия (6.55) и из формулы (6.57) при i=0 находим 90 c0 1 0h ; d0 Ah . (6.60) 1 1 На основании формул (6.59), (6.60) последовательно определяются коэффициенты c i , d i (i=1,…,n-1) до c n-1 и d n 1 включительно (прямой ход). Обратный ход начинается с определения y n . Для этого из второго краевого условия (6.55) и из формулы (6.57) при i=n-1 найдем yn 0h 1c n 1 d n 1 0h 1 (c n 1 1) . (6.61) Далее по формуле (6.57) последовательно находим y n 1 , y n 2 , ..., y 0 . Заметим, что метод прогонки обладает устойчивым в ычислительным алгоритмом. 7 Приближенное решение дифференциальных уравнений в частных производных В реальных физических процессах искомая функция з ависит от нескольких переменных, а это приводит к уравн ениям в частных производных от искомой функции. Как и для обыкновенных дифференциальных уравн ений(ОДУ), в этом случае для выбора одного конкретного решения, удовлетворяющего уравнению в частных производных , кроме начальных условий, необходимо задавать дополнительные условия (т.е. краевые условия ). Чаще всего такие задачи на практике не имеют аналитического решения и приходится использовать численные методы их решения, в том числе метод сеток, метод конечных разностей и так далее. Мы будем рассматривать класс линейных уравнений в частных производных второго порядка. В общем виде в случае двух переменных эти уравнения записываются в виде 2 A( x,y ) u x 2 2 B ( x,y ) u x y 2 C ( x,y ) u y 2 a ( x,y ) u x (7.1) 91 u b ( x,y ) y F ( x,y ) , c ( x,y )u где: A, B, C, a, b, c − заданные непрерывные функции двух переменных, имеющие непрерывные частные прои зводные, u – искомая функция. Для сокращения записи введем об означения u xx u ; u xy x2 2 u ; u yy x y 2 u; u x y2 u; x uy u. y Будем рассматривать упрощенную форму записи (7.1) вида A( x, y)u xx b( x, y)u y B( x, y)u xy c( x, y)u C ( x, y)u yy a( x, y)u x F ( x, y) ( 7.2) и рассмотрим частный случай (7.2), когда a=b=c=F 0, т.е. A( x, y)u xx B( x, y)u xy C ( x, y)u yy 0. (7.3) Путем преобразований уравнение (7.3) может быть приведено к каноническому виду (к одному из трех ста ндартных канонических форм) эллиптическому типу, гипе рболическому типу, параболическому типу. Причем тип уравнения будет определяться коэффициентами А, В, С, а именно – знаком дискриминанта D=B 2 −4 A C. Если D <0, то имеем уравнение эллиптического типа в точке с координатами x, y; если D=0, то (7.3) − параболического типа; если D>0, то (7.3) − гиперболического типа; если D не сохраняет постоянного знака, то (7.3) − смешанного типа. Замечание. Если А, В, С − константы, тогда каноническое уравнение (7.3) называется полностью эллиптическ ого, параболического, гиперболического типа. Введем понятие оператора Лапласа для сокращенной записи канонических уравнений вида Δu 2 u 2 x 2 u . 2 y 92 Используя это определение, запишем сокращенные канонические уравнения всех трех типов 1. u=0. Это уравнение эллиптического типа, так наз ываемое уравнение Лапласа. В механике это уравнение оп исывает стационарные тепловые поля, установившееся т ечение жидкости и т.д. 2. u=−f , где f − заданная непрерывная функция. Это уравнение Пуассона имеет эллиптический тип и опис ывает процесс теплопередачи с внутренним источником тепла. 3. a 2Δ u u t , где a − константа. Не во всех уравнениях в качестве переменных будут выступать станда ртные переменные x, y. Может быть также переменная времени. Это уравнение диффузии описывает процесс теплопрово дности и является уравнением параболич еского типа. 4. 2 u 2 t a 2Δ u , а– константа. Это уравнение гипербол и- ческого типа − так называемое волновое уравнение и оно описывает процесс распространения волн. 7.1 Метод сеток для решения смешанной задачи для уравнения параболического типа (уравнения теплопроводности) Смешанная задача означает, что следует найти иск омую функцию, удовлетворяющую заданному уравнению в частных производных, краевым, а так же начальным условиям. Различить эти условия можно в том случае, если одна из независимых переменных − время, а другая − пространственная координата. При этом условия, относящиеся к н ачальному моменту времени, называются начальными, а условия, относящиеся к фиксированным зн ачениям координат − краевыми. Рассмотрим смешанную задачу для однородного ура внения теплопроводности u t 2 k u x2 , k=const>0. (7.4) Задано начальное условие (7.5) 93 u( x,0) f ( x) и заданы краевые условия первого рода u (0 ,t ) μ1 (t ) ; u ( a ,t ) μ 2 (t ) . (7.6) Требуется найти функцию u(x,t), удовлетворяющую в области D (00. для любых значений параметра Есть и явная схема (рисунок 11), но она устойчива только при λ 1 2 , т.е. при τ h 2 2 . Вычисления по этой схеме придется вести с малым шагом по , что приводит к большим затратам машинного времени. Рисунок 11 – Четырехточечный шаблон явной схемы 95 7.2 Решение задачи Дирихле для уравнения Лапласа методом сеток Рассмотрим уравнение Лапласа 2 2 u 2 x u 2 y 0. (7.8) Будем рассматривать уравнение Лапласа в прямоугол ьной области Ω u(0, y) ( x, y ) , 0 x f1( y) ; u(a, y) a, 0 y b с краевыми условиями f 2 ( y) ; u( x,0) f3( x) ; u( x, b) f4( x) , где f1, f 2, f3, f 4 − заданные функции. Заметим, что чаще вс его область бывает не прямоугольной. Введем обозначения u ij =u(x i ,y j ). Накладываем на прямоугольную область сетку xi h i ; i=0,1,…,n; y j l j ; j=0,1,…,m. Тогда xn h n , ym l m =b. Частные производные аппроксимируем по формулам 2 u x 2 2 u i, j 1 2 u y ui ui 1, j O ( h 2); 2 u i, j 2 1 h 2 u i, j l u i, j 2 1 O (l 2 ) , и заменим уравнение Лапласа конечно-разностным уравнением Рисунок 12 – Схема узлов «крест» u i 1, j 2 u i, j u i 1, j u i, j h 2 1 2 u i, j u i, j l 2 1 0, (7.9) 96 где: i=1,…,n-1, j=1,...,m-1 (т.е. для внутренних узлов). Погрешность замены дифференциального уравнения разностным составляет величину О( h 2 +l 2 ). Уравнения (7.9) и значения u i,j в граничных узлах образуют систему л инейных алгебраических уравнений относительно прибл иженных значений функции u(x,y) в узлах сетки. Выразим u i,j при h=l, и заменим систему ( ui u ij ui 1, j u i0 f 3 ( x i ); u im f 4 ( x i ); u0 j f 1 ( y j ); u nj f 2 ( y j ). ui, j 1, j 1 u i , j 1) / 4 ; (7.10) Систему (7.10) линейных алгебраических уравнений можно решить любым итерационным методом (Зейделя , простых итераций и т.д.). При построении системы использовалась схема типа «крест» (рисунок 12). Строим последовательность итер аций по методу Зейделя ( s 1) u i,j 1 4 (s) (s (u i( s 1,j1) u i( s )1,j u i,j 1 u i,j 1) 1 ), где s– текущая итерация. Условие окончания итерационн ого процесса max u ij( s 1) u ij( s ) . (7.11) i, j Условие (7.11) ненадежно и на практике используют другой критерий max u ij( s 1) u ij( s ) (1 v ) , i,j где v max u ij( s max u (s) ij u ij( s ) 1) u ( s 1) ij . 97 Схема «крест»– явная устойчивая схема ( малое изменение входных данных ведет к малому изменению выходных данных). 7.3 Решение смешанной задачи для уравнения ги перболического типа методом сеток Рассмотрим уравнение колебания однородной и огран иченной струны. Задача состоит в отыскании функции u(x,t) при t>0, удовлетворяющей уравнению гиперболического типа 2 2 u t c 2 u x 2 , (7.12) где: 0