Конспект лекций Новосибирск 2003



страница1/4
Дата12.09.2017
Размер1.58 Mb.
Просмотров63
Скачиваний0
ТипКонспект
  1   2   3   4


Министерство образования Российской Федерации
НОВОСИБИРСКИЙ ГОСУДАРСТВЕННЫЙ ТЕХНИЧЕСКИЙ УНИВЕРСИТЕТ


ЧИСЛЕННЫЕ МЕТОДЫ

Конспект лекций

Новосибирск

2003


Главное отличие численных методов от аналитических в том, что действия выполняются не с выражениями, а со значениями, когда соотношение или функция имеют определенный вид и поиск решения заключается в применении алгоритма. Это любые действия с матрицами, решение систем линейных уравнений, нахождение корней, интегралов, производных и экстремумов функций заданного вида (чаще всего — полинома). Если вид функции заранее не известен, доступны лишь значения функции, выполняется либо перебор значений, либо подмена исследуемой функции функцией известного вида с коэффициентами, найденными по значениям исследуемой функции.
Численные методыАлгоритмическиеПеребораИтерационныеИнтерполяционныеРешение уравненийПрямого перебора

Половинного деленияПростых итерацийСекущих

Касательных

НьютонаРешение систем уравненийГаусса

LU-разложения

ПрогонкиЯкоби

ЗейделяНьютонаИнтегрированиеПрямоугольников

Трапеций


Симпсона

Сплайн-интерполяцииДифференцированиеЧисленного дифференцированияНахождение экстремумовПрямого перебора

ТрихотомииПараболической интерполяцииДругоеДействия над матрицами

Нахождение значения полинома большой степени



Некоторые алгоритмы работы с матрицами
Ниже приводится словесное и формульное описания некоторых алгоритмов работы с матрицами. Данная информация используется для разработки соответствующих программ и подпрограмм зада­ния.
Вычисление нормы матрицы

Норма матрицы ||A|| является скалярной характеристикой абсолютной величины ее элементов. Норма может быть построена разными способами, но наиболее распространенной является вычисляемая по формуле, аналогичной формуле для расчета модуля вектора:



В этой и в последующих формулах m – число строк, n – число столбцов матрицы. Можно использовать приведенную норму, которая рассчитывается на один элемент:



В матричной алгебре также используются m-норма, представляющая собой максимальную из сумм модулей элементов по строкам:



и l-норма, представляющая собой максимальную из сумм модулей элементов по столбцам:



.
Транспонирование матрицы

Операция транспонирования часто используется в матричной алгебре. Транспонированная матрица A обозначается при записи AT.

Транспонирование заключается в обмене местами строк и столбцов матрицы. Элемент, расположенный в i-той строке и j-том столбце помещается в j-тую строку и i-тый столбец.
Умножение матрицы на вектор

В результате умножения матрицы A из m строк и n столбцов на вектор , имеющий n компонентов, получается вектор , компоненты которого вычисляются по формуле:



, (i = 1, 2, …, m)
Умножение двух матриц

В результате умножения двух матриц и получает­ся матрица , элементы которой находятся по следующей прос­той формуле:



(i = 1, 2, …, m; j = 1, 2, …, p)

Для осуществления процедуры умножения необходимо, чтобы чис­ло столбцов первой матрицы A было обязательно равно числу строк второй матрицы B.


Возведение матрицы в целую степень

Возведение матрицы в целую степень производится путем ее повторного умножения на основе вышеприведенной формулы. Эту операцию можно выполнять только с квадратной матрицей.


Вычисление определителя квадратной матрицы

Для вычисления определителя квадратной матрицы произвольного порядка наиболее удобным является метод, основанный на преобразовании матрицы к треугольному виду и последующего вычисления произведе­ния элементов главной диагонали.

Алгоритм преобразования матрицы к треугольному виду, не из­меняющий значение определителя, заключается в следующем.

Преоб­разование выполняется за (n – 1) шагов (n - число строк и столбцов матрицы). На k-том шаге все элементы матрицы, лежащие ниже k-той строки и правее k-того столбца пересчитываются по формуле:



, где ,

(i = k+1, k+2, …, n; j = k+1, k+2, …, n)

В результате данного преобразования получается матрица, со­держащая ниже главной диагонали только нулевые элементы. Поэтому ее определитель равен произведению элементов главной диагонали:

Для работы алгоритма при возможных нулевых элементах в глав­ной диагонали (в исходной матрице или полученных при пересчете) рекомендуется использовать процедуру выбора главного элемента. Она заключается в перестановке строк матрицы на каждом шагу ее приведения к треугольному виду перед осуществлением пересчета. Отыскивается строка, содержащая в k-том столбце наибольший по аб­солютной величине элемент, а затем производится перестановка мес­тами k-той строки и строки с этим наибольшим элементом. Только после этого осуществляется деление на . Каждая перестановка меняет знак определителя на противоположный. С учетом перестановок



,

где m - число перестановок строк при выборе главного элемента.


Вычисление обратной матрицы

Один из возможных методов вычисления обратной матрицы осно­ван на решении матричного уравнения , где A - исходная, B - обратная и E - единичная матрицы:



Это уравнение при умножении матрицы A на каждый столбец мат­рицы B распадается на n систем линейных алгебраических уравнений (СЛАУ):



, где

(k = 1,2, …, n)

Решение каждой системы при k =1, 2 , ..., n позволяет определить k-тый столбец искомой обратной матрицы, т.е. процедура на­хождения обратной матрицы сводится к n-кратному решению СЛАУ. Этот метод называется методом Жордана.

Вычисление матричных функций

После определения операции возведения квадратной матрицы в целую степень можно ввести понятие матричной функции матричного аргумента, используя представление функций в виде степенного ряда. Например, по аналогии с экспоненциальной функцией скалярного аргумента x



вводится экспоненциальная функция матричного аргумента A:



Аналогично можно представить sin(A), cos(A) и другие функции:



Таким образом, матричная функция матричного аргумента представляет собой квадратную матрицу того же порядка, что и аргумент A. Она получается бесконечным суммированием матриц-слагаемых. При расчете суммирование членов ряда следует выполнять до тех пор, пока норма очередного слагаемого не станет меньше заданной погрешности .



Решение систем линейных уравнений прямыми методами
1. Введение

Система линейных алгебраических уравнений (СЛАУ), записанная в векторно-матричной форме, имеет вид:

(1)

где даны - матрица коэффициентов СЛАУ и - вектор свободных членов, а нужно определить - вектор неизвестных.

Для решения таких систем применяются как прямые (регулярные), так и ите­рационные методы /1, 2, 3, 5/. В вычислительной практике наиболее распространены прямые методы Гаусса и LU-разложения. В лабораторной работе изучаются эти методы, и с помощью реализующих их подпрограмм осуществляется наглядное решение СЛАУ, которое помогает понять их суть.

В ряде задач вычислительной математики, связанных с модели­рованием физических процессов, приходится решать системы линейных алгебраических уравнений с трехдиагональной матрицей. В этом слу­чае вместо методов Гаусса или LU-разложения более рационально ис­пользовать специальный метод - прогонки, который также изучается в данной работе.



2. Метод Гаусса

Метод Гаусса основан на замене исходной системы (1) на экви­валентную (т.е. имеющую то же решение ), но с матрицей верхней треугольной формы:

(2)

или в развернутом виде:

(2а)

Матрица получается путем преобразования исходной матрицы , а вектор - из вектора . Так как матрица - треугольная, все компоненты вектора , начиная с , легко находятся последо­вательно.

Алгоритм метода Гаусса состоит из двух этапов:


  1. Преобразование и . - прямой ход;

  2. Последовательное вычисление – обратный ход.

Далее приводится словесно-формульное описание алгоритмов прямого и обратного хода метода Гаусса. Обоснование проводимых операций можно найти в соответствующей литературе / 1, 2, 3 /.

1. Прямой ход

Процедура прямого хода осуществляется за шагов ( - поря­док СЛАУ). На k-том шаге (k = 1, 2, 3, … , n) формируется k-тая строка матрицы и k-тый компонент вектора . Каждый шаг состоит из двух частей.

1.1. На первой части шага вычисляются ис­комые элементы матрицы и компоненты вектора :

(3)

Примечание: Индекс вверху (в скобках) означает индекс шага, от которого следует взять соответствующие значения элементов матрицы и компонентов вектора . Он необходим, так как на каждом шаге эти величины преобразуются. Индекс соответствует исходным элемен­там и компонентам. В дальнейшем будет видно, что этот индекс яв­ляется фиктивным, т.е. при программировании не порождает соот­ветствующий индекс массивов. Но при алгебраической записи формул в литературе он обычно приводится.

1.2. На второй части шага оставшаяся часть матрицы (ниже k-той строки и правее k-того столбца) и вектора (ниже k-того компо­нента) пересчитываются по формулам:



(4)

При вторая часть шага не выполняется, т.к. пересчитывать больше нечего.

2. Обратный ход

Обратный ход начинается с определения , которое следует из (2а):

(5)

Затем последовательно находятся по фор­муле:

(6)

Эта формула следует из обобщения уравнений, которые получа­ются при умножении матрицы на вектор из (2а) и рассмотрении их в последовательности снизу вверх.

Действия с компонентами векторов и в (3) и (4) совпадают с аналогичными действиями над элементами матриц и , поэтому в алгоритме Гаусса часто присоединяют эти вектора в качестве допол­нительного (n+1)-го столбца соответствующих матриц. При этом ал­горитм и реализующая его программа становятся более компактными и красивыми.

3. Метод LU-разложения

Метод LU-разложения основан на том, что матрица может быть представлена в виде произведения двух матриц:

(7)

где - нижняя треугольная матрица,



- верхняя треугольная матрица с единичной главной диагональю.

При этом разложении исходная система (1) разлагается на две, но с треугольными матрицами:

(8а)

(8б)

которые легко решаются последовательно, т.е. сначала из решения системы (8a) находится вектор промежуточных неизвестных , а затем, используя найденный вектор в качестве вектора свободных членов, решается система (8б) и находится искомый вектор неизвестных .

1. Прямой ход

Алгоритм разложения (факторизации) матрицы на две треугольных также осуществляется за шагов. На k-том шаге формирует­ся k-тая строка матрицы и k-тый столбец матрицы . Каждый шаг состоит из двух частей.

1.1. На первой части вычисляются



(9а)

(9б)

1.2. На второй части шага пересчитывается оставшаяся часть матри­цы , лежащая ниже k-той строки и правее k-того столбца:



(10)
2. Обратный ход

2.1. Сначала решается система . Нижний треугольный вид мат­рицы позволяет последовательно найти :



(11)
2.2. Затем решается система . Так как она совпадает с сис­темой, получаемой в методе Гаусса, то алгоритм ее решения есть алгоритм обратного хода метода Гаусса.

4. Связь методов Гаусса и LU-разложения

Методы Гаусса и LU-разложения имеют глубокую связь. Подробно об этом можно прочитать, например, в / 2 /. Отметим только, что матрицы и векторы в обоих методах одинаковы. Из (8a) следует, что .

Методы Гаусса и LU-разложения имеют одинаковую трудоемкость по количеству вычислительных операций и примерно одинаково широко используются на практике. LU-разложение матрицы не затрагивает вектор свободных членов , поэтому этот метод удобно использовать при решении СЛАУ с одной и той же матрицей, но с разными вектора­ми свободных членов.

5. Выбор главного элемента

Методы Гаусса и LU-разложения не будут работать, если в главной диагонали окажется нулевой элемент, стоящий там изначаль­но или полученный при пересчете. Кроме того, если диагональный элемент не равен точно нулю, но мал по величине по сравнению с другими элементами матрицы, то при делении на него возникает большая погрешность результата.

Поэтому метод Гаусса и LU-разложения можно дополнять выбором главного элемента. Суть этой процедуры заключается в том, что на каждом шаге преобразования матриц сначала осуществляется поиск в оставшейся части k-го столбца (начиная с k-той строки) элемента, максимального по абсолютной величине. После его нахождения строка с максимальным элементом и текущая k-тая строка матрицы меняются местами.

В методе Гаусса производится и соответствующий обмен компо­нентов вектора свободных членов. Такой обмен не изменяет корни СЛАУ и их порядок в векторе .

Процедура выбора главного элемента усложняет алгоритм, но зато исключает деление на ноль или малую величину. Она обычно используется и при вычислении определителя матрицы путем приведения ее к треугольному виду.
6. Метод прогонки

Метод прогонки применяется для СЛАУ с трехдиагональной мат­рицей, в которой все элементы, не входящие в главную и прилегаю­щие к ней диагонали, равны нулю:





При этом нет необходимости хранить всю матрицу и коэффициен­ты СЛАУ хранятся и обрабатываются в виде трех одномерных масси­вов: - массива элементов главной диагонали, - массива элемен­тов ниже главной диагонали и - массива элементов выше главной диагонали. Индекс элемента массивов , или соответствует ин­дексу строки матрицы.



(12)

Идея метода прогонки основана на представлении каждой компо­ненты вектора неизвестных через следующую:



(13)

где и - прогоночные коэффициенты, которые необходимо опре­делить. Умножение первой строки СЛАУ (12) на вектор и приведение полученного уравнения к виду (13) позволяет найти выражения для первых прогоночных коэффициентов:



(14)

Формулы для последующих прогоночных коэффициентов можно по­лучить, умножая последовательно строки матрицы СЛАУ (12) на век­тор и приводя записи к виду (13):



(15)

где ( i =2, 3, ..., n)

После определения прогоночных коэффициентов совершается об­ратный ход - последовательное определение компонент вектора неиз­вестных, начиная с :

(i = n-1, n-2, …, 1) (16)
7. Погрешности решения СЛАУ

Казалось бы прямые методы решения СЛАУ являются достаточно точными; погрешность результата должна быть порядка погрешности вычисления - выполнения арифметических операций, которая определяется типом используемых данных. Для типа double константа машинного нуля DBL_EPSILON, характеризующая эту погрешность, составляет . Для методов Гаусса и LU-разложения выполняется порядка операций, поэтому ожидаемая погрешность даже при n = 100 должна быть порядка . Но на самом деле, величина погрешности сильно зависит и от конкретных значений элементов матрицы.

Каждая матрица, например , имеет так называемое число обусловленности . В вычислительной математике его определяют через нормы прямой и обратной матриц соотношением

Эта величина может изменяться от до . Чем она больше, тем больше будет погрешность вычисления вектора неизвестных

(17)

Число обусловленности определяет также связь погрешности вычисления вектора неизвестных с погрешностями задания матрицы коэффициентов и вектора свободных членов, которые в реальных задачах часто берутся из физических измерений и могут иметь значительную погрешность

; (18)

(19)

Реальная погрешность будет определяться большей из трех величин, получаемых по соотношениям (17), (18) и (19). Системы с большим значением числа обусловленности называются плохо обусловленными. Для них получение достаточно точного решения невозможно или представляет сложную проблему.

Для определения числа обусловленности матрицы используются весьма сложные алгоритмы /5, 6/ . В данной лабораторной работе для этих целей применяется функция Decomp( ), которая взята из источника / 6 /.


8. Расчет определителя матриц

Преобразование матрицы в прямом ходе алгоритмов Гаусса и LU-разложения позволяет легко вычислить ее определитель. В методе LU-разложения определитель матрицы есть произведение определителей сомножителей:



(20)

Так как матрицы и треугольные, их определители равны произведению элементов главной диагонали. При этом (главная диагональ матрицы содержит единицы) и



(21)

В методе Гаусса главная диагональ матрицы после ее окончательного пересчета совпадает с главной диагональю матрицы , поэтому



(22)

Алгоритм расчета определителя трехдиагональной матрицы более сложен. Он основывается на рекурентной формуле, которая может быть получена разложением определителя по его минорам / 8 /:



; (i = n, n-1, n-2, … , 3) (23)

В этих формулах , и - элементы соответствующих массивов-диагоналей.

Так как язык С/С++ допускает рекурсию, то вычисление определителя трехдиагональной матрицы удобно оформить в виде следующей рекурсивной функции:

double det3d(Vector R, Vector P, Vector Q, int n)

{

if (n==1) return P[1];

if (n==2) return P[1]*P[2]-Q[1]*R[2];
return P[n]*det3d(R,P,Q,n-1)-Q[n-1]*R[n]*det3d(R,P,Q,n-2);

}

Решение нелинейных алгебраических уравнений
1. Введение

В научно-технических расчетах часто встречаются нелинейные уравнения вида



, (1)

где F(x) - некоторая нелинейная (трансцендентная) функция. Вид ее не позволяет выразить x и, тем самым, найти корень аналитическими методами. Численные методы позволяют найти приближенное значение корня уравнения с заданной точностью. Некоторые из них требуют задания интервала, на котором предполагается корень, другие - начального (нулевого) приближения корня. Значение корня ищется путем после­довательных приближений (итераций), которые циклически совершают­ся по определенной итерационному алгоритму. Число итераций заранее определить нельзя – они совершаются до выполнения некоторого условия. В лабораторной работе изучаются методы прос­тых итераций, Ньютона, половинного деления и хорд.


Методы,_требующие_задания_интервала_поиска_корня'>2. Основные методы решения нелинейных уравнений
2.1. Методы, требующие задания интервала поиска корня

Проблема выбора достаточно точного нулевого приближения, обеспечивающего сходимость итераций описанных выше методов, для некоторых уравнений становится трудноразрешимой. Зато часто можно легко указать интервал, в котором ожидается корень уравнения. Поэтому широко используются методы, основанные на за­дании интервала поиска, к которым относятся метод половинного де­ления и метод хорд.

Суть метода половинного деления (бисекции) заключается в том, что интервал, на котором находится корень, последовательно делится пополам. Знак левой части уравнения в середине интервала сравнивается со знаком на его границах. Далее, в среднюю точку переносится та граница, где знак левой части уравнения совпадает со знаком на середине интервала; при этом интервал поиска умень­шается в два раза. Затем все повторяется. Итерации продолжаются до тех пор, пока ширина интервала не станет меньше заданной погрешности , а значение левой части уравнения на середине интервала не бу­дет по абсолютной величине меньше погрешности . Тогда за иско­мый корень принимается середина интервала.

Метод обладает высокой надежностью и позволяет найти корни уравнений, не решаемых методами простых итераций и Ньютона.

Блок схема алгоритма метода половинного деления представлена на рис. 5.3.

Метод хорд в чем-то похож на метод половинного деления. Суть метода заключается в следующем. На каждой итерации по двум точкам - границам интервала – и левая часть уравнения линеаризиру­ется и находится корень полученного линейного уравнения

(9)

Затем та граница интервала ( или ), где знак левой части уравнения F совпадает со знаком ее в точке x, переносится в эту точку и все повторяется до тех пор, пока не будут выполняться ус­ловия (2). Блок-схема алгоритма метода хорд представлена на рис. 5.4.

Метод хорд часто путают с методом секущих, который ближе к методу Ньютона и использует для вычисления производной в формуле (7) значения функции на двух последних итерациях.
2.2. Методы, требующие задания нулевого приближения

Методы, требующие задания нулевого приближения, работают по следующей схеме. Задается нулевое приближение корня , затем, по некоторой итерационной формуле находится первое приближение, по нему - второе и т.д.



Итерации выполняются до тех пор, пока не будет выполняться условия:



(2a)

(2b)

где - заданная абсолютная погрешность определения корня;



- заданная погрешность невязки левой части уравнения с нулем.

Из методов этого типа наиболее распространены методы простых итераций и Ньютона.


Метод простых итераций применим для уравнений, которые можно привести к виду

, (3)

из которого вытекает итерационная формула:



(4)

Итерации выполняются по ней до тех пор, пока не будет выпол­няться условие (2a). Условие (2b) в этом методе не может быть использовано.

Итерационный процесс сходится к значению корня, если на ин­тервале итераций выполняется условие:

(5)
Метод Ньютона является наиболее быстро сходящимся и поэтому одним из наиболее используемых.

Суть метода основана на том, что левая часть уравнения в окрестности линеаризируется по вычисленным значениям самой левой части и ее первой производной



(6)

и находится корень полученного линейного уравнения



,

который и принимается за следующее приближение .

Это приводит к итерационной формуле:

(7)

Итерации выполняются из заданного нулевого приближения до тех пор, пока не будут выполняться условия (2). Блок-схема алгоритма метода Ньютона приведена на рис. 5.2.

Метод Ньютона в своем классическом варианте не всегда сходится. Итерации сходятся к значению корня лишь тогда, когда на интервале итераций знаки первой и второй производных от левой части уравнения не изменяются. Для других случаев существуют спе­циальные разновидности метода Ньютона с более сложными алгоритмами.

Простым и эфффективным способом добиться сходимости во многих случаях является ограничение ньютоновского приращения. Если величина приращения на k-той итерации



получается больше некоторого задаваемого предельного значения , то оно искусственно ограничивается этой величиной:



(8)

При ограничении знак приращения должен быть выбран равным знаку приращения .



Поделитесь с Вашими друзьями:
  1   2   3   4


База данных защищена авторским правом ©zodorov.ru 2017
обратиться к администрации

войти | регистрация
    Главная страница


загрузить материал