Итерационные методы решений систем алгебраических уравнений

Точные методы

 

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

Метод Гаусса

Прямой ход метода

1-й шаг. Предположим, что а₁₁¹0. Поделим первое уравнение на этот элемент:

 

. (1.3)

Остальные уравнений системы (1.1) запишем в виде

 

, (1.4)

 

где i= .

Уравнение (1.3) умножаем на a_i1 и вычитаем из i-го уравнения системы (1.4).

Получим систему вида:

 

. (1.5)

 

.

 

Система (1.5) имеет матрицу вида:

 

.

 

Работаем с укороченной системой, т.к. х₁ входит только в 1-ое уравнение

 

.

2-й шаг. Если , то из укороченной системы аналогично исключаем неизвестное x₂ и получаем матрицу коэффициентов такого вида:

 

.

 

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

х_m-1 и приходим к системе с матрицей вида:

 

. (1.6)

 

Эта матрица является верхней треугольной:

 

.

 

Обратный ход метода. Из последнего уравнения системы (1.6) находим х_m, из предпоследнего х_m-1,..., из первого уравнения – х₁.

Общая формула:

 

x_m=y_m/c_mm,

, (i=m-1,…,1).

 

Для реализации метода Гаусса требуется примерно (2/3)m³ арифметических операций, причем большинство из них приходится на прямой ход.

Ограничение метода единственного деления заключается в том, что угловые элементы не равны нулю, т.е. .

Ведущие элементы – – элементы на k-ом шаге исключения. Но если ведущий элемент близок к нулю, то в процессе вычисления может накапливаться погрешность. В этом случае на каждом шаге исключают не х_k, a х_j (при j¹k). Такой подход называется методом выбора главного элемента. Для этого выбирают неизвестные x_j с наибольшим по абсолютной величине коэффициентом либо в строке, либо в столбце, либо во всей матрице. Для его реализации требуется - арифметических действий.

 

1.1.2 Связь метода Гаусса с разложением матрицы на множители. Теорема об LU разложении.

 

Пусть дана система A х =B (1.1), которая при прямом ходе преобразуется в эквивалентную систему (1.6) и которую запишем в виде

C х = y, (1.7)

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

 

Как связаны в системе (1.1) элементы В и элементы y из (1.7)?

Если внимательно посмотреть на прямой ход метода Гаусса, то можно увидеть, что

 

.

 

Для произвольного j имеем

 

, (1.7)

 

где j= , d_ji – числовые коэффициенты:

 

. (1.8)

 

Это можно записать в виде:

 

В=D y,

где D -нижняя треугольная матрица с элементами на главной диагонали (j= , ).

 

В связи с тем, что в методе Гаусса угловые коэффициенты не равны нулю , то на главной диагонали матрицы D стоят не нулевые элементы. Следовательно, эта матрица имеет обратную, тогда y=D^-1В, С x = D^-1В. Тогда

D´C x =B. (1.9)

 

В результате использования метода Гаусса, получили разложение матрицы А на произведение двух матриц

A = D´C,

где D - нижняя треугольная матрица, у которой элементы на главной диагонали не равны нулю, а C - верхняя треугольная матрица с единичной диагональю.

Таким образом, если задана матрица A и вектор B, то в методе Гаусса сначала производится разложение этой матрицы А на произведение D и C, а затем последовательно решаются две системы:

 

D y =B,

C x = y. (1.10)

 

Из последней системы находят искомый вектор x. При этом разложение матрицы А на произведение СD – есть прямой ход метода Гаусса, а решение систем (1.10) обратный ход. Обозначим нижнюю треугольную матрицу через L, верхнюю треугольную матрицу - U.

 

Теорема об LU разложении

 

Введем обозначения: - угловой минор порядка j матрицы А, т.е.

 

 

Теорема. Пусть все угловые миноры матрицы А не равны нулю (Δj¹0 для j= ). Тогда матрицу А можно представить единственным образом в виде произведения А=L*U.

 

Идея доказательства. Рассмотрим матрицу А второго порядка и будем искать разложение этой матрицы в виде L и U.

 

.

 

Сопоставляя эти два равенства, определяем элементы матриц L и U (перемножим и приравняем неизвестные). Система имеет единственное решение. Методом математической индукции сказанное можно обобщить для матрицы размерности m´m.

Следствие. Метод Гаусса (схема единственного деления) можно применять только в том случае, когда угловые миноры матрицы А не равны нулю.

 

1.1.3 Метод Гаусса с выбором главного элемента

 

Может оказаться так, что система (1.1) имеет единственное решение, хотя какой либо из миноров матрицы А равен нулю. Заранее неизвестно, что все угловые миноры матрицы А не равны нулю. Избежать этого можно с выбором главного элемента.

1. Выбор главного элемента по строке, т.е. производится перенумерация неизвестных системы.

 

ПРИМЕР. Пусть дана система второго порядка

 

 

при ½ а₁₂ ½>½ а₁₁ ½, тогда на первом шаге вместо неизвестного х₁ исключают х₂:

 

.

 

К этой системе применяем первый шаг прямого хода метода Гаусса.

2. Выбор главного элемента по столбцу.

Предположим, что ½ а₂₁ ½>½ а₁₁ ½и переставим уравнения

 

и применяем первый шаг прямого хода метода Гаусса. Здесь имеет место перенумерация строк.

3. Поиск главного элемента по всей матрице заключается в совместном применении методов 1 и 2. Всё это приводит к уменьшению вычислительной погрешности.

 

 

1.1.4 Метод Холецкого (метод квадратных корней)

 

Пусть дана система

 

А х =В, (1.11)

 

где матрица А - симметричная матрица. Тогда решение системы (1.11) проводится в два этапа:

1. Симметричная матрица А представляется как произведение двух матриц

 

А = L * L^Т.

 

Рассмотрим метод квадратных корней на примере системы 4-го порядка:

 

.

 

Перемножаем матрицы в левой части разложения и сравниваем с элементами в левой части:

 

, , , , ,

 

, , .

 

2. Решаем последовательно две системы

 

L y =B,

L^T x = у.

 

Особенности.

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

2) Возможно переполнение - если угловые элементы близки к нулю.

 

 

Рассмотрим систему

 

A x= B,

 

где А - невырожденная действительная матрица.

Метод простых итераций

 

Систему (3.1) преобразуем к следующему эквивалентному виду:

 

.

(3.3)

Или в векторной форме

(3.4)

 

Пусть задано начальное приближение . Подставляем его в правую часть системы (3.4) и получаем x⁽¹⁾=j(x⁽⁰⁾), продолжая подстановку, находим х⁽²⁾ и т.д. Получим последовательность точек , которая приближается к исходному решению х.

 

3.1.1 Условия сходимости метода.

 

Пусть j '(x) - матрица Якоби (якобиан), соответствующая системе (3.4) и в некоторой D-окрестности решения х функции (i=1,2,…,m) дифференцируемы и выполнено неравенство вида:

 

,

где (0 £ q < 1), q - постоянная.

 

Тогда независимо от выбора х⁽⁰⁾ из D-окрестности корня итерационная последовательность {х^k} не выходит за пределы данной окрестности, метод сходится со скоростью геометрической прогрессии и справедлива оценка погрешности

.

 

3.1.2 Оценка погрешности.

 

В данной окрестности решения системы, производные функции j_i(x) (i=1,…,m) должны быть очень малы по абсолютной величине, т.е. сами функции должны быть почти постоянными. Тогда исходную систему (3.1) следует преобразовать к виду (3.3) с учетом условий сходимости.

Пример. Рассмотрим предыдущий пример и приведем систему к удобному для итераций виду

 

Проверяем условие сходимости вблизи точки С. Вычислим матрицу Якоби

 

.

 

Так как x₁»3.8, x₂»2, то при этих значениях вычисляем норму матрицы

 

|| ||»|| ||»0.815.

 

Запишем итерационную процедуру

 

Следовательно, метод простых итераций будет сходиться со скоростью геометрической прогрессии, знаменатель которой q»0.815. Вычисления поместим в таблице 1.

 

Таблица 1 Решение системы нелинейных уравнений

К			…
	3.80000	3,75155	….	3,77440	x1=3,77418
	2.00000	2,03895	…	2,07732	x2=2,07712

 

При К=9 критерий окончания счета выполняется при e=10^-3 и можно положить x₁ =3.774 0.001

x₂ =2.077 0.001.

Метод Ньютона

 

Суть метода состоит в том, что система нелинейных уравнений сводится к решению систем линейных алгебраических уравнений. Пусть дана система (3.1) и задано начальное приближение x ⁽⁰⁾, приближение к решению х строим в виде последовательности .

В исходной системе (3.1) каждую функцию где i= , раскладывают в ряд Тейлора в точке х⁽ⁿ⁾ и заменяют линейной частью её разложения

 

.

 

Для каждого уравнения получаем

 

.

(3.5)

 

В матричной форме

 

(3.6)

где f ' - матрица Якоби.

 

Предположим, что матрица не вырождена, то есть существует обратная матрица .

Тогда система (3.6) имеет единственное решение, которое и принимается за очередное приближение x ⁽ⁿ⁺¹⁾. Отсюда выражаем решение x ⁽ⁿ⁺¹⁾ по итерационной формуле:

 

.

(3.7)

Сходимость метода

 

Теорема. Пусть в некоторой окрестности решения х системы (3.1) функции f_i (при i= ) дважды непрерывно дифференцируемы и матрица Якоби не вырождена. Тогда найдется такая малая d окрестность вокруг решения х, что при выборе начального приближения x ⁰ из этой окрестности итерационный метод (3.7) не выйдет за пределы этой окрестности решения и справедлива оценка вида

,

 

где n - номер итерации.

Метод Ньютона сходится с квадратичной скоростью. На практике используется следующий критерий остановки:

 

.

Прямые методы

4.1.1 Метод Леверрье

 

Метод разделяется на две стадии:

- раскрытие характеристического уравнения,

- нахождение корней многочлена.

 

Пусть det(A-lE) - есть характеристический многочлен матрицы А ={a_ij} (i,j=1,2,…,m), т.е. , и l₁,l_2,…,l_m - есть полная совокупность корней этого многочлена (полный спектр собственных значений).

Рассмотрим суммы вида (k=1,2,…,m), т.е.

 

,

(4.3)

где - след матрицы.

В этом случае при k£m справедливы формулы Ньютона для всех (1£k£ m)

,

(4.4)

 

Откуда получаем

при k=1 р1 = -S1 , при k=2 р2 = -1/2 * (S2 + р1*S1), .............. при k=m рm = -1/n * (Sm + р1*Sm-1 + р2*Sm-2 +... + рm-1*S1).

(4.5)

 

Следовательно, коэффициенты характеристического многочлена р _i можно определить, если известны суммы S₁,S₂,...,S_m. Тогда схема алгоритма раскрытия характеристического определителя методом Леверрье будет следующей:

1) вычисляем степень матрицы: А^к=А^к-1*А для k=1,…,m;

2) определяют S_k - суммы элементов стоящих на главной диагонали матриц А^к;

3) по формулам (4.5) находят коэффициенты характеристического уравнения р_i (i=1,2,…,m).

4.1.2 Усовершенствованный метод Фадеева

 

Алгоритм метода:

1) вычисляют элементы матриц A₁,A₂,..,A_m:

 

 

(в конце подсчета B_m нулевая матрица для контроля);

 

2) определяют коэффициенты характеристического уравнения р_i

q₁ = -р₁, q₂ = -р₂,..., 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 имеют одни и те же собственные числа. Это можно легко увидеть раскрыв определитель

 

.

Следовательно, матрицы A и B - подобные, имеют одни и те же собственные значения. Но собственные векторы х и х’ – не совпадают, они связаны между собой простым соотношением

 

х = S* х '.

 

Такую матрицу A с помощью преобразования подобия или же последовательности таких преобразований можно привести к матрице Фробениуса вида:

.

 

Детерминант матрицы F det (F) можно разложить по элементам первой строки:

 

.

 

Тогда коэффициенты характеристического многочлена матрицы А будут р₁ = f₁₁ , p₂ = f_12,…, p_n = f_1m.

 

Второй случай. Матрицу А преобразованием подобия можно привести к матрице В верхнего треугольного вида

.

 

Тогда собственными числами будут диагональные элементы матрицы B:

 

.

 

Третий случай. Матрицу A с помощью преобразования подобия можно привести к Жордановой форме

 

,

 

где l_i - собственные числа матрицы A; S_i - константы (0 или 1); если S_i=1, то l_i=l_i+1.

 

К четвёртому случаю относятся матрицы, которые с помощью преобразования подобия можно привести к диагональному виду (матрица простой структуры):

 

,

 

Задача приближения функции

 

Постановка задачи. Пусть на отрезке [a,b] задана функция у= f(x) своими n+1 значениями , в точках .

Допустим, что вид функции f(x) неизвестен. На практике часто встречается задача вычисления значений функции у= f(x) в точках х, отличных от . Кроме того, в некоторых случаях, не смотря на то, что аналитическое выражение у=f(x) известно, оно может быть слишком громоздким и неудобным для математических преобразований (например, специальные функции). Кроме этого значения y_i могут содержать ошибки эксперимента.

Определение. Точки называются узлами интерполяции.

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

 

 

Определение. Процесс вычисления значений функции F(x) в точках отличных от узлов интерполирования называется интерполированием функции f(x). Если , то задача вычисления приближенного значения функции в т. х называется интерполированием, иначе - экстраполированием.

Геометрически задача интерполирования функции одной переменной означает построение кривой, проходящей через заданные точки (рисунок 5). То есть задача в такой постановке может иметь бесконечное число решений.

у

x x x x х

 

Рисунок 5 - Геометрическая иллюстрация задачи интерполирования функции

 

Задача становится однозначной, если в качестве F(x) выбрать многочлен степени не выше n, такой что:

 

F (x )=y , F (x )=y ,..., F (x )=y .

 

Определение. Многочлен F (x), отвечающий вышеназванным условиям, называется интерполяционным многочленом.

Определение. Когда многочлен F(x) выбирается в классе степенных функций, то интерполяция называется параболической.

Знание свойств функции f позволяет осознанно выбирать класс G аппроксимирующих функций. Широко используется класс функций вида

 

(5.1)

 

являющихся линейными комбинациями некоторых базисных функций (x),..., _m(x). Будем искать приближающую функцию в виде многочлена степени m, с коэффициентами с ,...,с , которые находятся в зависимости от вида приближения. Функцию Ф_m(х) называют обобщенным многочленом по системе функций j₀(х),j₁(х),…,j_m(х), а число m – его степенью. Назовем обобщенный многочлен Ф_m(х) интерполяционным, если он удовлетворяет условию

 

Ф_m(х_i)=y_i, (i=0,1,…,n). (5.2)

 

Покажем, что условие (5.2) позволяет найти приближающую функцию единственным образом

 

(5.3)

 

Система (5.3) есть система линейных алгебраических уравнений относительно коэффициентов с₀,с₁,…,с_m.

Эта система n линейных уравнений имеет единственное решение, если выполняется условие m=n и определитель квадратной матрицы Р

 

 

Определение. Система функций (x),..., _m(x), линейно независимая в точках х₀,х₁,…,х_n, которые попарно различны и выписанный выше определитель не равен нулю, называется Чебышевской системой функций. Если мы имеем такую систему, то можно утверждать, что существует единственный для данной системы функций интерполяционный многочлен Ф_m(х), коэффициенты которого определяются единственным образом из системы (5.3).

Пример. При m£n система функций 1 ,х,х²,…,х^m линейно независима в точках х₀,х₁,…,х_n, если они попарно различны.

Отсюда следует, что

 

,

 

где знаменатель не должен быть равен нулю, т. е.

 

.

(6.47)

 

Если условие (6.47) выполнено, то краевая задача (6.35), (6.36) имеет единственное решение. Если же (6.47) не выполняется, то краевая задача (6.35), (6.36) либо не имеет решения, либо имеет бесконечное множество решений.

 

 

Задано начальное условие

 

(7.5)

 

и заданы краевые условия первого рода

 

(7.6)

 

Требуется найти функцию u(x,t), удовлетворяющую в области D (0< x a, 0< t T) условиям (7.5) и (7.6). Физически это можно представить как стержень, на концах которого поддерживается требуемый температурный режим, заданный условиями (7.6).

Рисунок 10 – Неявная схема

При проведении замены получим , т.е. k =1. Задача решается методом сеток: строим в области D равномерную сетку с шагом h по оси x и шагом t по t (см. рисунок 10).

Приближенное значение искомой функции в точке - обозначим через . Тогда ; ; i =0,1,..., n; ;

j =0,1,..., m; .
Заменим производные разностными отношениями

 

;

.

 

В результате получим неявную двухслойную схему с погрешностью O(t+h²)

 

.

 

Используя подстановку , выразим из этой схемы u_i,j-1

 

,

(7.7)

 

где: u_0,j=m₁ (t_j); u_n,j=m₂ (t_j).

Получаем разностную схему, которой аппроксимируем уравнение (7.4). Эта схема (7.7) неявная, и выглядит так, как показано на рисунке 10. При построении схемы (7.7) получается система линейных уравнений с трехдиагональной матрицой. Решив ее любым способом (в частности, методом прогонки), получаем значения функции на определенных временных слоях. Так, на нулевом временном слое используем начальное условие U_i,0=f (x_i), т.к. j =0. Эта неявная схема более устойчива для любых значений параметра l>0.

Есть и явная схема (рисунок 11), но она устойчива только при , т.е. при .

Рисунок 11 - Явная схема

Лабораторная работа № 1

Варианты заданий

 

Таблица 2

Но­мер вари­анта	Матрица системы	Правая часть
	.4000 -.0029 -.0055 -.0082	.0003 -.5000 -.0050 -.0076	.0008 -.0018 -1.4000 -.0070	.0014 -.0012 -.0039 -2.3000	.1220 -.2532 -.9876 -2.0812
	1.7000 .0000 -.0003 -.0005	.0003 .8000 -.0002 -.0004	.0004 .0001 -.1000 -.0003	.0005 .0002 .0000 -1.0000	.6810 .4803 -.0802 -1.0007
	3.0000 .0011 -.0005 -.0022	.0038 2.1000 .0005 -.0011	.0049 .0032 1.200C -.0001	.0059 .0043 .0026 .3000	1.5136 1.4782 1.0830 .3280
	4.3000 .0100 .0037 -.0027	.0217 3.4000 .0090 .0027	.0270 .0207 2.5000 .0080	.0324 .0260 .0197 1.6000	2.6632 2.7779 2.5330 1.9285
	5.6000 .0147 .0087 .0028	.0268 4.7000 .0150 .0090	.0331 .0271 3.8000 .0153	.0393 .0334 .0274 2.9000	4.0316 4.3135 4.2353 3.7969
	6.9000 .0191 .0134 .0077	.0319 6.0000 .0205 .0149	.0390 .0333 5.1000 .0220	.0461 .0405 .0348 4.2000	5.6632 6.1119 6.2000 5.9275
	8.2000 .0234 .0179 .0124	.0370 7.3000 .0260 .0205	.0451 .0396 6.4000 .0286	.0532 .0477 .0422 5.5000	7.5591 8.1741 8.4281 8.3210
Продолжение таблицы 2
	9.5000 .0278 .0224 .0170	.0422 8.6000 .0315 .0261	.0513 .0459 7.7000 .0351	.0604 .0550 .0496 6.8000	9.7191 10.5000 10.9195 10.9775
	10.8000 .0321 .0268 .0215	.0475 9.9000 .0369 .0316	.0576 .0523 9.0000 .0416	.0676 .0623 .0570 8.1000	12.1430 13.0897 13.6744 13.8972
	12.1000 .0365 .0312 .0260	.0528 11.2000 .0423 .0370	.0639 .0586 10.3000 .0481	.0749 .0697 .0644 9.4000	14.8310 15.9430 16.6926 17.0800
	13.4000 .0408 .0356 .0304	.0581 12.5000 .0477 .0425	.0702 .0650 11.6000 .0546	.0822 .0770 .0718 10.7000	17.7828 19.0599 19.9744 20.5261
	14.7000 .0452 .0400 .0349	.0635 13.8000 .0531 .0479	.0765 .0714 12.9000 .0610	.0896 .0844 .0793 12.0000	20.9985 22.4406 23.5195 24.2353
	16.0000 .0496 .0444 .0393	.0688 15.1000 .0585 .0534	.0829 .0777 14.2000 .0674	.0970 .0918 .0867 13.3000	24.4781 26.0849 27.3281 28.2078
	17.3000 .0539 .0488 .0437	.0741 16.4000 .0639 .0588	.0892 .0841 15.5000 .0739	.1043 .0992 .0941 14.6000	28.2215 29.9928 31.4001 32.4435
	23.8000 .0757 .0707 .0656	.1010 22.9000 .0909 .0858	.1212 .1161 22.0000 .1060	.1414 .1363 .1313 1.1000	50.8968 53.4873 55.7118 57.5703
	19.9000 .0626 .0576 .0525	.0849 19.0000 .0747 .0696	.1020 .0969 18.1000 .0867	.1191 .1140 .1090 17.200	36.5001 38.5997 40.3345 41.7045
	21.2000 .0670 .0619 .0569	.0902 20.3000 .0801 .0750	.1084 .1033 19.4000 .0932	.1265 .1215 .1164 18.5000	41.0351 41.2986 45.1968 46.7299
Продолжение таблицы 2
	22.5000 .0714 .0663 .0612	.0956 21.6000 .0855 .0804	.1148 .1097 20.7000 .0996	.1339 .1289 .1238 19.8000	45.8340 48.2611 50.3226 52.0184
	23.8000 .0757 .0707 .0656	.1010 22.9000 .0909 .0858	.1212 .1161 22.0000 .1060	.1414 .1363 .1313 21.1000	50.8968 53.4873 55.7118 57.5703
	25.1000 .0801 .0750 .0700	.1063 24.2000 .0963 .0912	.1276 .1225 23.3000 .1124	.1488 .1437 .1387 22.4000	56.2234 58.9772 61.3645 63.3853
	26.4000 .0844 .0794 .0744	.1117 25.5000 .1017 .0966	.1339 .1289 24.6000 .1189	.1562 .1512 .1461 23.7000	61.8139 64.7307 67.2806 69.4636
	27.7000 .0888 .0838 .0788	.1171 26.8000 .1070 .1020	.1403 .1353 25.9000 .1253	.1636 .1586 .1536 25.0000	67.6682 70.7478 73.4601 75.8051
	29.0000 .0932 .0882 .0831	.1225 28.1000 .1124 .1074	.1467 .1417 27.2000 .1317	.1710 .1660 .1610 26.3000	73.7864 77.0286 79.9030 82.4098
	30.3000 .0975 .0925 .0875	.1278 29.4000 .1178 .1128	.1531 .1481 28.5000 .1381	.1784 .1734 .1684 27.6000	80.1684 83.5730 86.6095 89.2778
	31.6000 .1019 .0969 .0919	.1332 30.7000 .1232 .1182	.1595 .1545 29.8000 .1445	.1859 .1809 .1759 28.9000	86.8143 90.3811 93.5793 96.4090
	32.9000 .1062 .1013 .0963	.1386 32.0000 .1286 ..1236	.1659 .1610 31.1000 .1510	.1933 .1883 .1833 30.2000	93.7240 97.4528 100.8126 103.8034
	34.2000 .1106 .1056 .1006	.1400 33.3000 .1340 .1290	.1724 .1674 32.4000 .1574	.2007 .1957 .1907 31.5000	100.8.976 104.7881 108.3093 111.4610
Продолжение таблицы 2
	35.5000 .1150 .1100 .1050	.1494 34.6000 .1394 .1344	.1788 .1738 33.7000 .1638	.2082 .2032 .1982 32.8000	108.3351 112.3871 116.0694 119.3819
	36.8000 .1193 .1143 .1094	.1547 35.9000 .1448 .1398	.1852 .1802 35.0000 .1702	.2156 .2106 .2056 31.1000	116.0363 120.2497 124.0930 127.5660
	38.1000 .1237 .1187 .1137	.1601 37.2000 .1502 .1452	.1916 .1866 36.3000 .1766	.2230 .2180 .2131 35.4000	124.0015 128.3760 132.3800 136.0134

 

 

Лабораторная работа № 3

Лабораторная работа № 4

Лабораторная работа № 5

Точные методы

Метод Леверрье

 

Входные параметры: n—целое положительное число, равное порядку n системы;

Выходные параметры: b—массив из n действительных чисел, при выходе из программы содержит решение системы.

Пример. Вычислить собственные значения и собственные вектора матрицы

 

|2 2 –2|

|2 5 –4|

|-2 –4 5|

 

Текст программы:

 

PROCEDURE Liver (Var VeS:Tvector;Nn:Integer);;

Label 35,270,220,100,80,190,500,250,350,120,10;

Var

I,J,K,T:Integer; Y,Yy:tvector; {A,Aa:Mas;} S,C,Q,P,H,D,U,M,V,F,W,Yyy,Z,L,R,X,Aaa,B:Real;

Sq,Sq1: Real; F1: TEXT;

Begin

Assign(F1,'last.out');

Rewrite(F1);

For I:=Nn Downto 1 Do Y[i+1]:=VeS[i];

Y[1]:=1;

K:=1;

35: T:=1;

C:=Y[2]/Y[1];

If (Nn=1) Then Begin

P:=-C; Q:=0; Goto 270;

End;

If (Nn=2) Then Begin

H:=C*C/4-Y[3]/Y[1];

Goto 220;

End;

M:=10;C:=4;D:=8;U:=4;

V:=8;F:=1;W:=2;T:=0;

80: If(M<>10) Then Goto 100;

P:=C;M:=0;Q:=D;C:=U;D:=V;

U:=P;V:=Q;Yyy:=C;Z:=D;F:=-F;

100:M:=M+1;H:=0;Q:=Y[1];P:=Y[2]-C*Q;L:=Q;

120:For J:=3 To Nn Do Begin

R:=P;P:=Y[j]-C*R-D*Q;Q:=R;

R:=L;L:=Q-C*R-H*D;H:=R;

End;

Q:=Y[Nn+1]-D*Q;

S:=L+C*R;

If (T=0) Then Begin

X:=D*R;H:=R*X+S*L;

If H=0 Then Goto 80;

End;

C:=C+(P*S-Q*R)/H;

D:=D+(P*X+Q*L)/H;

If ((C-Yyy+D-Z)<>0) Then Goto 190;

If (F=-W) Then Begin writeln('Љ®аҐ­м ­Ґ ­ ©¤Ґ­');

Goto 500;

End;

W:=-F;

190:H:=C*C/4-D;

Sq:=Sqrt((Q-P*C/2)*(Q-P*C/2)+P*P*Abs(H));

If (Sq>0.0000001) Then Goto 80;

T:=0;

Y[2]:=Y[2]-C*Y[1];

For J:=3 To Nn-1 Do Begin

Y[j]:=Y[j]-C*Y[j-1]-D*Y[j-2];End;

220:P:=-C/2;

Q:=Sqrt(Abs(H));

If(h>=0) Then Begin

M:=P+Q;P:=P-Q;Q:=0;

Goto 250; End;

M:=P;

250:write(F1,'p',k,'=',m:9:4,' + i',q:7:4);K:=K+1;{Ves[k-1]:=M;}WriteLn(F1,'');

270:write(F1,'p',k,'=',p:9:4,' + i',-q:7:4);K:=K+1;{Ves[k-1]:=P;}WriteLn(F1,'');

If(t<>0) Then Goto 350;

Nn:=Nn-2;

Aaa:=0;

Sq1:=Sqrt((S-R*C/2)*(S-R*C/2)+R*R*Abs(H));

If (Sq1<=0.0000001) Then Aaa:=1;

B:=0;

If (Nn>=2) Then B:=1;

If((Aaa+B)=2) Then Begin

T:=1;

Goto 120; End; Goto 35;350:500:

Close(F1);

End;

 

Вычисления по программе привели к следующим результатам:

Собственные числа Собственные вектора

.10000Е+01 -.94281Е+00.23570Е+00 -.23570Е+00

.10000Е+02 -.33333Е+00 -.66667Е+00.66667Е+00

.10000Е+01.00000Е+00 -.70711Е+00 -.70711Е+00

 

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

Метод Фадеева

 

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

Выходные параметры: q—массив из n действительных чисел при выходе из программы содержит решение системы.

Метод Крылова

 

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

Выходные параметры: q—массив из n действительных чисел при выходе из программы содержит решение системы.

Схема алгоритма приведена на рисунке 24.

Пример. Вычислить собственные значения и собственные вектора матрицы

 

|2 2 –2|

|2 5 –4|

|-2 –4 5|

 

Текст программы:

 

PROCEDURE Krylov(Nn:integer; A: TMatr;Var Y:TVector);

Label 2;

Var I,J,K,T,Nn1:Integer;

yy:TVector;

aa:TMatr;

s,c,q,p,h,d,u,

m,v,f,w,yyy,z,

l,r,x,aaa,b: Real;

Begin

form3.Caption:='Крылов';

Nn1:=Nn;

For I:=1 To Nn Do Yy[i]:=I*Trunc(10*ranDom);

2:For I:=Nn Downto 1 Do Begin

Y:=Yy;

For J:=1 To Nn Do Aa[j,i]:=Yy[j];

For J:=1 To Nn Do Begin

S:=0;

For K:=1 To Nn Do S:=S+A[j,k]*Y[k];

Yy[j]:=S; End;

End;

For I:=1 To Nn Do Y[i]:=-Yy[i];

Simq(Nn1,Aa,Y,K