Размещено на http://www.allbest.ru/

ГОУ ВПО «Магнитогорский государственный университет»

Кафедра прикладной математики и вычислительной техники

КУРСОВАЯ РАБОТА

по дисциплине Уравнения математической физики

Метод вращений (Якоби) для нахождения собственных значений и собственных векторов матриц

Выполнила студентка 37 гр., ФМФ

Харина М.Н.

Научный руководитель: Торшина О. А.

Магнитогорск, 2010

Содержание

якоби вращение вектор матрица

1. Введение

2. Основные определения

2.1 Определения собственного значения и собственного вектора матрицы. Пример нахождения собственных значений и собственных векторов матрицы

2.2. Некоторые свойства собственных значений векторов

2.3 Системы линейных алгебраических уравнений

3. Методы решения систем линейных алгебраических уравнений

3.1 Метод Зейделя

3.2 Метод верхних релаксаций

3.3 Метод Якоби

4. Программа на C++ для решения СЛАУ методом Якоби

5. Заключение

6. Список использованной литературы

1. Введение

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

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

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

Одним из самых распространенных методов решения систем линейных уравнений является метод вращений Якоби. Этот метод (который также называют методом простых итераций) известен в различных вариантах уже более 200 лет.

2. Основные определения

2.1 Определения собственного значения и собственного и вектора матрицы

Рассмотрим квадратную матрицу n-ого порядка:

Собственные значения i квадратной матрицы A есть действительные или комплексные числа, удовлетворяющие условию:

,

E - единичная матрица,

- собственный вектор матрицы A, соответствующий некоторому собственному значению .

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

Определитель этой матрицы называется характеристическим определителем и равен:



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

и называется характеристическим многочленом. Корни этого многочлена - собственные значения или характеристические числа матрицы A. Числа называются коэффициентами характеристического многочлена.

Ненулевой вектор называется собственным вектором матрицы A, если эта матрица переводит вектор X в вектор

,

т.е. произведение матрицы A на вектор X и произведение характеристического числа на вектор X есть один и тот же вектор. Каждому собственному значению матрицы соответствует свой собственный вектор .

Для определения координат собственного вектора составляется характеристическое уравнение: . Переписав его в векторном виде и выполнив умножение, получим систему линейных однородных уравнений:

Определитель этой системы равен нулю, т.к. из этого условия были определены собственные значения матрицы A. Следовательно, система имеет бесконечное множество решений. Ее можно решить с точностью до постоянного множителя (как систему однородных уравнений). Решив эту систему, мы найдем все координаты собственного вектора X. Подставляя в систему однородных уравнений поочередно , получаем n собственных векторов.

При определении собственных значений и принадлежащих им собственных векторов решается одна и двух задач:

1. Определение все собственных значений и принадлежащих им собственных векторов матриц;

2. Определение одного или нескольких собственных значений и принадлежащих им собственных векторов.

Первая задача состоит в развертывании характеристического определителя в многочлен n-й степени (т.е. в определении коэффициентов ) с последующим вычислением собственных значений и, наконец, в определении координат собственного вектора .

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

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

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

Пример нахождения собственных значений и собственных векторов матрицы

Найти собственные значения и собственные векторы матрицы

Решение. Составляем характеристическую матрицу Размещено на http://www.allbest.ru/

:

Находим характеристический многочлен

Решим характеристическое уравнение

Подбором находим, что один корень уравнения равен -1. Есть теорема, которая говорит, что если число c является корнем многочлена , то многочлен делится на разность , то есть , где  -- многочлен. В соответствии с этой теоремой многочлен должен делиться на . Выделим в характеристическом многочлене этот множитель :

Находим корни трехчлена . Они равны -1 и 3. Таким образом,

 -- корень кратности 2 ,  -- простой корень. Итак, собственные числа матрицы A равны , . Найдем соответствующие им собственные векторы.

Пусть , тогда для собственного вектора получаем матричное уравнение

что соответствует системе уравнений

Решаем ее методом Гаусса. Выписываем расширенную матрицу системы

Первую строку, умноженную на числа -2 и -3 прибавляем соответственно ко второй и третьей строкам

Меняем местами вторую и третью строки

Возвращаемся к системе уравнений

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

Полагаем , находим , . Итак, собственному числу соответствует собственный вектор Размещено на http://www.allbest.ru/

.

Пусть , тогда для собственного вектора Размещено на http://www.allbest.ru/

получаем матричное уравнение

что соответствует системе уравнений

Решаем ее методом Гаусса. Выписываем расширенную матрицу

Первую строку умножаем на числа 2 и 3 и прибавляем соответственно ко второй и третьей строкам

Вторую строку умножаем на -1 и прибавляем к третьей

Возвращаемся к системе уравнений

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

Полагаем , находим , . Итак, собственному числу соответствует собственный вектор . Чтобы избавиться от дроби, умножим собственный вектор на 2, получим собственный вектор с тем же самым собственным числом. В итоге собственному числу соответствует собственный вектор .

Ответ: Собственные числа: , , соответствующие собственные векторы: , .

2.2 Некоторые свойства собственных значений векторов

Все n собственных значений любой симметричной матрицы (aij=aji; i,j = 1,2,…,n) вещественны.

Собственные векторы, отвечающие различным собственным значениям симметричной матрицы, ортогональны:

, при ,

, при .

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

Подобные матрицы

, где P - неособая матрица, имеют одинаковые собственные значения, их собственные вектора связаны соотношением:

2.3 Системы линейных алгебраических уравнений

Система m линейных уравнений с n неизвестными (или, линейная система) в линейной алгебре -- это система уравнений вида

(2.3.1)

Здесь x1, x2, …, xn -- неизвестные, которые надо определить. a11, a12, …, amn -- коэффициенты системы -- и b1, b2, … bm -- свободные члены -- предполагаются известными. Индексы коэффициентов (aij) системы обозначают номера уравнения (i) и неизвестного (j), при котором стоит этот коэффициент, соответственно.

Система (2.3.1) называется однородной, если все её свободные члены равны нулю (b1 = b2 = … = bm = 0), иначе -- неоднородной.

Система (2.3.1) называется квадратной, если число m уравнений равно числу n неизвестных.

Решение системы (2.3.1) -- совокупность n чисел c1, c2, …, cn, таких что подстановка каждого ci вместо xi в систему (1) обращает все её уравнения в тождества.

Система (2.3.1) называется совместной, если она имеет хотя бы одно решение, и несовместной, если у неё нет ни одного решения.

Совместная система вида (1) может иметь одно или более решений.

Решения c1(1), c2(1), …, cn(1) и c1(2), c2(2), …, cn(2) совместной системы вида (2.3.1) называются различными, если нарушается хотя бы одно из равенств:

c1(1) = c1(2), c2(1) = c2(2), …, cn(1) = cn(2).

Совместная система вида (2.3.1) называется определённой, если она имеет единственное решение; если же у неё есть хотя бы два различных решения, то она называется неопределённой. Если уравнений больше, чем неизвестных, она называется переопределённой.

3. Методы решения систем линейных алгебраических уравнений

3.1 Метод Зейделя

Метод Зейделя является классическим итерационным методом решения системы линейных уравнений.

Возьмём систему:

, где

Или

Перепишем задачу в виде:

Здесь в j-м уравнении мы перенесли в правую часть все члены, содержащие xi , для i > j. Эта запись может быть представлена:

где в принятых обозначениях D означает матрицу, у которой на главной диагонали стоят соответствующие элементы матрицы A, а все остальные нули; тогда как матрицы U и L содержат верхнюю и нижнюю треугольные части A, на главной диагонали которых нули.

Итерационный процесс в методе Зейделя строится по формуле после выбора соответствующего начального приближения .

Новые значения используются здесь сразу же по мере получения

где

Таким образом i-тая компонента (k + 1)-го приближения вычисляется по формуле:

Условие окончания итерационного процесса Зейделя при достижении точности в упрощённой форме имеет вид:

Более точное условие окончания итерационного процесса имеет вид

и требует больше вычислений. Хорошо подходит для разреженных матриц.

3.2 Метод верхних релаксаций

Метод верхней релаксации - это есть метод Зейделя с заданным числовым параметром w.

Метод верхних релаксаций является одним из наиболее распространенных одношаговых методов, который имеет следующий вид (3.2.3), где w заданный числовой параметр (0<w<2). Изменяя w можно получать различную скорость сходимости итерационного процесса. Этот параметр выбирается таким образом, чтобы на каждом шаге итерационного процесса уменьшалась величина, характеризующая близость полученного решения к искомому решению системы.

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

Для получения расчетных формул (3.2.3) перепишем в виде: или в компонентной записи получим (3.2.4) - это есть основная вычислительная формула.

В выражение (3.2.4) и входят одинаковым образом => при вычислениях они могут быть записаны в один и тот же массив. При реализации метода верхних релаксаций используется следующая форма записи алгоритма вычислений .

Действительно, при последовательном нахождении элемента (i+10 итерации) на каждом шаге будут использоваться найденные ранее значения, которые при k<j соответствуют i+1 итерации, а при k<j-i итерации.

3.3 Метод Якоби

Метод Якоби состоит из цепочки ортогональных преобразований подобия матрицы. Каждое преобразование (ротация Якоби) -- это плоский поворот с целью обнуления одного из внедиагональных элементов матрицы. Последовательные преобразования не сохраняют уже установленные нулевые элементы, но вместе с тем внедиагональные элементы становятся меньше и меньше до тех пор, пока матрица не станет диагональной с точностью до машинного нуля. Накопление в процессе преобразований произведения трансформационных матриц дает матрицу собственных векторов, в то время как диагональные элементы являются собственными значениями.

Метод Якоби всегда является робустным для действительных симметричных матриц. Для матриц размера, скажем больше 10, этот алгоритм более медленный, чем метод QR (см. ниже). Однако алгоритм Якоби значительно проще других более эффективных методов, таким образом он рекомендуется в случаях, когда время выполнения не является критическим.

Базовая матрица ротации Якоби имеет вид:

Ppq = (1 )

( … )

( c…s )

( …1… )

( -s…c )

( … )

( 1)

Все диагональные элементы матрицы ротации Якоби равны 1, за исключением двух элементов в строках (и столбцах) p и q. Все внедиагональные элементы равны нулю, за исключением двух элементов, равных s и (-s). Числа c и s являются косинусом и синусом угла поворота f, так что c2+s2=1.

Ротация Якоби преобразует матрицу A согласно формуле:

A' = PpqTAPpq.

Произведение PpqTA изменяет только строки p и q матрицы A, в то время как APpq меняет только столбцы p и q. Отметим, что нижние индексы p и q в обозначении Ppq не означают элемент этой матрицы, а скорее индексируют тип ротации, т.е. плоскость, в которой она происходит. Таким образом, измененные элементы в матрице A' расположены только в строках и столбцах p и q:

( ... a'1p ... a'1q ... )

(… ... ... ... )

(a'p1 ... a'pp ... a'pq ... a'pn)

A' =( ... ... ... ... )

(a'q1 ... a'qp ... a'qq ... a'qn)

(… ... ... ... )

( ... a'np … a'nq ... )

Для измененных элементов матрицы (с учетом приведенного выше и ее симметрии) получаются явные формулы:

a'rp = carp - sarq (r!=p, r!=q)

a'rq = carq + sarp (r!=p, r!=q)

a'pp = c2app + s2aqq - 2csapq

a'qq = s2app + c2aqq + 2csapq

a'pq = (c2 - s2)apq + cs(app- aqq)

Идея метода Якоби в том, чтобы постараться обнулить внедиагональные элементы серией ротаций. Для того, чтобы обнулилось a'pq, угол ротации, согласно формулам, приведенным выше, должен быть следующим:

q = cot(2f) = (c2 - s2)/(2cs) = (aqq- app)/(2apq).

Полагая t=s/c, для определения t прийдем к следующему уравнению:

t2 + 2tq - 1 = 0.

Меньший по модулю корень этого уравнения соответствует вращению на угол, меньший p/4; выбор этого угла на каждой из стадий дает наиболее устойчивый алгоритм. Используя формулу для решения квадратного уравнения с дискриминантом в знаменателе, меньший корень определяется как:

t = sign(q)/(|q|+sqrt(q2+1)).

Если величина q такова, что q2 даст переполнение, полагаем t = 1/(2q). По значению t определяются c и s:

c = 1/sqrt(t2+1), s = tc.

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

a'pp = app - tapq

a'rp = arp - s(arq + garp)

a'rq = arq + s(arp - garq) ,

где t (тангенс половины угла поворота) определяется, как t = s/(1+c).

Сходимость метода Якоби можно оценить, рассматривая сумму квадратов внедиагональных элементов S; уравнения преобразований дают для этой суммы после трансформации следующее соотношение:

S = Sr!=s(|ars|2); S' = S - 2|apq|2.

Поскольку преобразования ортогональные, сумма диагональных элементов при этом возрастает соответственно на |apq|2. Сумма же квадратов внедиагональных элементов монотонно убывает. Поскольку она ограничена снизу нулем и поскольку мы можем каждый раз выбирать какой хотим элемент apq для обнуления, сумма будет стремиться к нулю.

После цепочки преобразований получается матрица D, диагональная с точностью до машинного нуля. Ее диагональные элементы образуют собственные значения первоначальной матрицы A, поскольку выполняется

D = VTAV, где V = P1P2P3...,

а Pi -- матрицы ротации Якоби. Столбцы матрицы V являются собственными векторами. Они вычисляются рекуррентно, на каждой стадии по формулам:

V' = VPi ,

где первоначально V -- единичная матрица. Для преобразования элементов V используются формулы:

v'rs = vrs (s!=p, s!=q)

v'rp = cvrp - svrq

v'rq = svrp + cvrq

Для уменьшения ошибки округления вышеприведенные формулы следует переписать в терминах g (как выше для элементов матрицы A).

Единственным остающимся вопросом теперь является стратегия, по которой нужно перебирать уничтожаемые элементы. Оригинальный алгоритм Якоби от 1846 года на каждой стадии искал в верхней треугольной области наибольший по модулю элемент и обнулял его. Это весьма рациональная стратегия для ручного счета, но она неприемлема для компьютера, поскольку делает число операций на каждой стадии ротации Якоби порядка N2 вместо N.

Для наших целей лучшей стратегией является циклический метод Якоби, где элементы удаляются в строгом порядке. Например, можно просто проходить по строкам: P12 , P13 , ... P1n ; P23 ,...,P2n ;... Можно показать, что скорость сходимости будет в общем случае квадратичной, как для оригинального, так и для циклического метода Якоби. Назовем одно такое множество последовательных n(n-1)/2 преобразований Якоби проходом.

Программа, помещенная ниже, содержит две тонкости:

· На первых трех проходах мы проводим ротацию относительно индекса (pq) только тогда, когда |apq|>e для некоторого порогового значения e = S0/(5n2), где S0 -- сумма модулей внедиагональных элементов верхней треугольной матрицы.

· После четырех проходов в том случае, когда |apq|<<|app| и |apq|<<|aqq|, полагаем |apq|=0 и пропускаем ротацию. Критерием сравнения является превышение диагональных элементов в 10D+2 раз, где D -- число значащих десятичных цифр.

В приведенной ниже программе симметричная матрица a размером (n x n) на входе загружена в массив a[1...n][1...n]. На выходе наддиагональные элементы a разрушаются, но диагональные и поддиагональные остаются на месте с тем, чтобы сохранить информацию об исходной симметричной матрице. Вектор d[1...n] возвращает собственные значения a. Во время вычисления он содержит текущую диагональ. Матрица v[1...n][1...n] возвращает набор нормализованных собственных векторов, относящихся к d[k] для ее k-го столбца. Параметр nrot -- число ротаций Якоби, необходимых для достижения сходимости.

Типичные матрицы требуют от 6 до 10 проходов для достижения сходимости, или от 3n2 до 5n2 ротаций. Каждая ротация требует порядка 4n операций умножения и сложения, так что общие затраты имеют порядок от 12n3 до 20n3 операций. При вычислении собственных векторов наряду с собственными значениями, число операций на ротации возрастает с 4n до 6n, что означает превышение только на 50 процентов.

4. Программа на C++ для решения СЛАУ методом Якоби

#include <math.h>

#include "nrutil.h" /* Здесь определяются некоторые утилиты типа выделения памяти */

/* Преобразование элементов при ротации */

#define ROTATE(a,i,j,k,l) g=a[i][j];h=a[k][l];a[i][j]=g-s*(h+g*tau);a[k][l]=h+s*(g-h*tau)

/* максимальное число проходов */

#define MAXSWEEP 50

/* Программа jacobi вычисляет все собственные значения и собственные векторы

действительной симметричной матрицы a[1...n][1...n]. На выходе разрушаются

все наддиагональные элементы a. d[1...n] возвращает собственные значения

матрицы, v[1...n][1...n] -- в столбцах нормализованные собственные векторы,

nrot -- число ротаций Якоби, потребовавшихся для данной программы.

*/

void jacobi(float **a, int n, float *d, float **v, int *nrot) {

int j, iq, ip, i;

float tresh, theta, tau, t, sm, s, h, g, c, *b, *z;

/* выделить память для временно используемых векторов (декларации в nrutil.h) */

b=vector(1,n); z=vector(1,n);

/* инициализировать v как единичную матрицу */

for(ip=1;ip<=n;ip++) {

for(iq=1;iq<=n;iq++) v[ip][iq]=0.;

v[ip][ip]=1.;

}

/* инициализировать b диагональю a, z -- нулем */

for(ip=1;ip<=n;ip++) {b[ip]=a[ip][ip]; z[ip]=0.;}

/* вначале число ротаций нулевое */

*nrot=0;

/* делаем не более MAXSWEEP проходов */

for(i=1;i<=MAXSWEEP;i++) {

/* вычисляем сумму модулей внедиагональных элементов */

for(sm=0.,ip=1;ip<=n;ip++) for(iq=ip+1;iq<=n;iq++) sm += fabs(a[ip][iq]);

/* диагональная матрица -> нормальный выход */

if(sm==0.) {

free_vector(z,1,n); free_vector(b,1,n); return;

}

/* пороговое значение элемента для ротации */

tresh=(i<4 ? 0.2*sm/(n*n) : 0.);

/* проход осуществляется по строкам, в каждой строке по столбцам */

for(ip=1;ip<=n-1;ip++) for(iq=ip+1;iq<=n;iq++) {

/* отследить случай малого элемента после 4 проходов */

g=100.*fabs(a[ip][iq]);

if(i>4 && (float)fabs(d[ip]+g)==(float)fabs(d[ip])

&& (float)fabs(d[iq]+g)==(float)fabs(d[iq])) a[ip][iq]=0.;

/* и случай малого элемента на первых 3 проходах

(обработать только превысившие порог) */

else if(fabs(a[ip][iq])>tresh) {

h=d[ip]-d[iq];

/* вычислить значение t=s/c по формуле корня квадратного уравнения */

if((float)(fabs(h)+g)==(float)fabs(h)) t=a[ip][iq]/h;

else {

theta=0.5*h/a[ip][iq];

t=1./(fabs(theta)+sqrt(1.+theta*theta));

if(theta<0.) t = -t;

}

/* вычислить c, s, tau, и др. Изменить диагональ. Обнулить (ip,iq)-элемент */

c=1./sqrt(1+t*t); s=t*c; tau=s/(1.+c); h=t*a[ip][iq];

z[ip] -= h; z[iq] += h; d[ip] -= h; d[iq] += h;

a[ip][iq]=0.;

/* поворот при 1<=j<ip */

for(j=1;j<=ip-1;j++) {ROTATE(a,j,ip,j,iq);}

/* поворот при ip<j<iq */

for(j=ip+1;j<=iq-1;j++) {ROTATE(a,ip,j,iq,j);}

/* поворот при iq<j<=n */

for(j=iq+1;j<=n;j++) {ROTATE(a,ip,j,j,iq);}

/* добавка для матрицы собственных векторов */

for(j=1;j<=n;j++) {ROTATE(v,j,ip,j,iq);}

/* приращение счетчика ротаций */

++(*nrot);

}

}

/* добавить до диагонали и реинициализировать z */

for(ip=1;ip<=n;ip++) {

b[ip] += z[ip]; d[ip]=b[ip]; z[ip]=0.;

}

}

/* если мы здесь, то число ротаций превысило лимит. Функция nerror (выход с диагностикой ошибки) описана декларацией в nrutil.h. */

nerror("Too many iterations in the routine jacobi");

}

5. Заключение

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

Но такие системы уравнений могут быть, во-первых, очень большого размера, например, NxN=10000х10000, и даже более; во-вторых, система уравнений может оказаться недоопределенной; в-третьих, она может оказаться с линейно зависимыми уравнениями; в-четвертых, она может оказаться переопределённой и несовместной. Кроме того, в-пятых, вычислительная техника может иметь далеко не рекордное быстродействие и объём оперативной памяти, и заведомо конечную разрядность двоичного представления чисел и связанные с этим ненулевые вычислительные погрешности. Поэтому итерационные методы, и метод вращений в частности, получили большое применение в решении СЛАУ. Современная вычислительная техника позволяет проводить исследование устойчивости и сходимости итерационного метода в зависимости от параметров задачи.

6. Список использованной литературы

1 Вержбицкий В.М. Численные методы (математический анализ и обыкновенные дифференциальные уравнения). - М.: Высшая школа, 2001.

2.Демидович Б.П., Марон И.А. Численные методы анализа. - М.: Наука, 1970.

3.Мысковских И.П. Лекции по методам вычислений, СПб. 1998

4. Самарский А.А. «Введение в численные методы». Москва «Наука», 1987.

5. Фаддеев Д.К., Фаддеева В.Н. «Вычислительные методы линейной алгебры». Москва «Физматгиз», 1963

Размещено на Allbest.ru