Alg i metodi vichisl / Teorija / lek_num_met

xk = xk −1 − F(xk −1) ,

F '(xk −1)

F (x) = a0 + a1x +...+ an xn , F '(x) = a1 + 2a2 x +...+ nan xn−1.

Для вычисления значений многочленов F(x) и F'(x) в точке х=xk-1 может быть использована схема Горнера. Естественно, при использовании метода Ньютона должны выполняться условия сходимости. При их соблюдении в результате численного решения получается значение того корня, который находится вблизи заданного начального приближения x0.

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

Комплексные корни.

При использовании компьютера имеется возможность работать с комплексными числами; поэтому изложенный метод Ньютона может быть использован (с необходимым обобщением) и для нахождения

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

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

(3.17)

Линейный остаток	равен нулю, если р, q выражаются с помощью найденных корней:
	р = - 2с, q =	х = с ± id.

Представление (3.17) может быть также использовано для нахождения р, q, а значит, и для определения корней. Эта процедура лежит в основе метода Лина. Суть этого метода состоит в следующем. Предположим,

что коэффициенты	равны нулю. Тогда, сравнивая коэффициенты при одинаковых	степенях х
многочлена F(x) в выражениях (3.16) и (3.17), можно получить (для упрощения выкладок		)

(3.18)

…………

(3.19)

В соотношения (3.19) входят коэффициенты и которые являются функциями p и q. Действительно,

задав значения р и q, из соотношений (3.18) можно последовательно найти и. Поэтому соотношения (3.19) представляют собой систему двух нелинейных уравнений относительно р и q:

Такая система в методе Лина решается методом простой итерации: задаются начальные приближения для р, q

которые используются для вычисления коэффициентов затем из уравнений (3.19) уточняются значения р, q. Итерационный процесс вычисления этих величин продолжается до тех пор, пока их изменения в двух последовательных итерациях не станут малыми.

Широко распространен также другой метод, основанный на выделении квадратичного множителя х2 + рх+q,

— метод Бэрстоу. Он использует метод Ньютона для решения системы двух уравнений.

§3 Системы уравнений. Метод простой итерацию. Метод Ньютона.

3.1 Системы уравнений.

Вводные замечания.

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

Пусть для вычисления неизвестных	требуется решить систему п нелинейных уравнений
F(x ,x ,...,x )=0,
1 1 2	n
F2(x1,x2,...,xn)=0,		(15)
. . . . . . . . . . .		(3.20)
Fn(x1,x2,...,xn)=0.
В векторной форме эту систему можно записать как

F(x) = О,

где
F={	}, x={	}.

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

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

3.2 Метод простой итерации.
Систему уравнений (3.20) представим в виде
x1 = f1(x1,x2,...,xn),
x2 = f2(x1,x2,...,xn),	(16)
. . . . . . . . . . .	(3.21)

xn = fn(x1,x2,...,xn)

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

значений на предыдущей итерации	как
xi(k )	= fi (x1(k −1) , x2(k −1) ,…, xn(k −1) ),i =1,2,…,n	(3.22)

Систему (3.21) можно решать и методом Зейденя, напоминающим метод Гаусса-Зейделя решения систем линейных уравнений. Значение находится из i-го уравнения системы (3.21) с использованием уже

вычисленных на текущей итерации значений неизвестных. Таким образом, значения неизвестных на k-й. итерации будут находиться не с помощью (3.22), а с помощью соотношения

xi(k ) = f i (x1(k ) ,…, xi(−k1) , xi(k −1) ,…, xn(k −1) ), i = 1,2,…, n

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

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

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

Этот метод обладает гораздо более быстрой сходимостью, чем метод простой итерации и метод Зейделя. В случае одного уравнения F(x)= 0 алгоритм метода Ньютона был легко получен путем записи уравнения касательной к кривой у =F(x). По сути для нахождения нового приближения функция F(x) заменялась линейной функцией, т. е. раскладывалась в ряд Тейлора, при этом член, содержащий вторую производную, отбрасывался (как и все последующие члены). Та же идея лежит в основе метода Ньютона для системы

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

Рисунок 3.6 !!!!!!!

Пусть приближенные значения неизвестных системы (3.20), полученные на предыдущей итерации, равны

соответственно	Задача состоит в нахождении приращений (поправок) к этим
значениям	благодаря которым следующее приближение к решению системы (3.20)
запишется в виде
	(3.23)

Проведем разложение левых частей уравнений (3.20) в ряд Тейлора, ограничиваясь лишь линейными членами относительно приращений:

F1 (x1(k ) , x2(k ) ,…, xn(k ) ) ≈ F1 +	∂F1	x1 +…+	∂F1		xn ,
	∂x		∂x	n
	1

F2 (x1(k ) , x2(k ) ,…, xn(k ) ) ≈ F2 +	∂F2
	∂x1
. . . . . . . . . . .
Fn (x1(k ) , x2(k ) ,…, xn(k ) ) ≈ Fn	+ ∂Fn
	∂x
	1

x1 +	∂F2	xn	,
	… + ∂xn			(3.24)
. . . . . . . . . . . . .
x1	+…+ ∂Fn		xn .
	∂xn

В правых частях этих соотношений значения и их производных вычисляются в точке

x(k−1) = (x1(k−1) , x2(k−1) ,…, xn(k−1) )

Поскольку в соответствии с (3.20) левые части (3.24) должны обращаться в нуль, приравняем нулю и правые

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

∂F1 ∂x1 ∂F2 ∂x1

∂Fn ∂x1

x1	+	∂F1	x2	+…+	∂F1	xn = −F1 ,
		∂x2			∂xn
x1	+	∂F2	x2	+…+	∂F2	xn = −F2 ,
		∂x2			∂xn
		∂Fn	. . . . . . . . . . . .
x1	+		x2	+…+	∂Fn	xn = −Fn ,	(3.25)
		∂x2			∂xn

Определителем системы (3.25) является якобиан.

Для существования единственного решения системы (3.25) он должен быть отличным от нуля на каждой итерации.

Таким образом, итерационный процесс решения системы уравнений (3.20) методом Ньютона состоит в

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

,

(3.26)

Пусть приближенные значения неизвестных равны а, b. Предположим, что якобиан системы (3.26) при х=а, у=b отличен от нуля, т. е.

	∂F1		∂F1
J =	∂x		∂y	≠ 0
	∂F		∂F
	2	2
	∂x		∂y

Тогда следующие приближения неизвестных можно записать в виде

x =a −	1		(F		∂F2		−F		∂F1	),
	J				∂y
		1				2			∂y
y =b +	1		(F		∂F2		−F		∂F1		).
	J				∂x
		1				2			∂x

Величины, стоящие в правых частях, вычисляются при х= а, у= b.

Алгоритм метода Ньютона для решения системы двух уравнений изображен на рис. 3.6. В качестве исходных

данных задаются начальные приближения неизвестных х, у, погрешность и допустимое число итераций М. В условной конструкции, проверяющей выполнение критерия завершения итерации, используется логическая операция «и», имеющаяся в современных языках программирования. Если итерации сойдутся, то выводятся значения х, у; в противном случае — текущие значения х, у и соответствующее сообщение.

Глава 4 ПРЕДСТАВЛЕНИЕ МАТРИЦ И МНОГОМЕРНЫХ МАССИВОВ

НА ЯЗЫКАХ ВЫСОКОГО УРОВНЯ

§ 1. Представление матриц и многомерных массивов на языках C, C++.

Специального типа данных матрица или многомерный массив в Си (Си++) нет, однако, можно использовать массив элементов типа массив. Например, переменная a представляет матрицу размера 3×3 с вещественными элементами:

double a[3][3];

Элементы матрицы располагаются в памяти последовательно по строкам: сначала идут элементы строки с индексом 0, затем строки с индексом 1, в конце строки с индексом 2 (в программировании отсчет индексов всегда начинается с нуля, а не с единицы!). При этом выражение

a[i],

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

a[i][j],

представляет собой элемент матрицы a в строке с индексом i и столбце с индексом j. Элемент матрицы можно использовать в любом выражении как обычную переменную (например, можно читать его значение или присваивать новое).

Такая реализация матрицы удобна и максимально эффективна с точки зрения времени доступа к элементам. У нее только один существенный недостаток: так можно реализовать только матрицу, размер которой известен заранее. Язык Си не позволяет описывать массивы переменного размера, размер массива должен быть задан до начала стадии компиляции.

Пусть нужна матрица, размер которой определяется во время работы программы. Тогда пространство под нее надо захватывать в динамической памяти с помощью функции malloc языка Си или оператора new языка C++. При этом в динамической памяти захватывается линейный массив и возвращается указатель на него. Рассмотрим вещественную матрицу размером n строк m на столбцов. Захват памяти выполняется с помощью функции malloc языка Си

double *a;

. . .

a = (double *) malloc(n * m * sizeof(double));

или с помощью оператора new языка C++:

double *a;

int m, n;

. . .

a = new double[n * m];

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

синдексом m − 1. Каждая строка состоит из n элементов, следовательно, индекс элемента строки i и столбца j

влинейном массиве равен

i * n + j

(действительно, поскольку индексы начинаются с нуля, то i равно количеству строк, которые нужно пропустить, i * n - суммарное количество элементов в пропускаемых строках; число j равно смещению внутри последней строки). Таким образом, элементу матрицы в строке i и столбце j соответствует выражение

*(a + i * n + j)

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

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

double **a; // Адрес массива указателей

int m, n; // Размеры матрицы: m строк, n столбцов int i;

. . .

// Захватывается память под массив указателей a = (double **) malloc(m * sizeof(double *));

for (i = 0; i < m; ++i) {

// Захватывается память под строку с индесом i *(a+i) = (double *) malloc(n * sizeof(double));

}

После этого к элементу a ij можно обращаться с помощью выражения

*(*(a+i)+j)

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

Многомерные массивы реализуются аналогично матрицам. Например, вещественный трехмерный массив размера 4 × 4 × 2 описывается как

double a[4][4][2];

обращение к его элементу с индексами x, y, z осуществляется с помощью выражения

a[x][y][z]

Многомерные массивы переменного размера с числом индексов большим двух встречаются в программах довольно редко, но никаких проблем с их реализацией нет: они реализуются аналогично матрицам. Например, пусть надо реализовать трехмерный вещественный массив размера m × n × k. Захватывается линейный массив вещественных чисел размером m * n * k:

double *a;

. . .

a = (double *) malloc(m * n * k * sizeof(double));

Доступ к элементу с индексами x, y, z осуществляется с помощью выражения

a[(x * n + y) * k + z]

§ 2. Представление матриц и многомерных массивов на языке Pascal

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

Type

Vector = array [1..5] of <тип_элементов>;

Matrix= array [1..10] of vector;

Var m: matrix;

Type

Matrix= array [1..5] of array [1..10] of < тип элементов >;

или еще проще:

type

matrix = array [1..5, 1..10] of <тип элементов>;

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

m [i, j]

Только для инициализации двумерного массива используется вложенный цикл for. Например:

. . . .

For i:= 1 to 10 do For j:= 1 to 20 do

m[i, j] := 0;

Многомерные массивы реализуются аналогично матрицам. N-мерный массив характеризуется N индексами. Формат описания такого типа данных:

Type
<Имя	типа>	=	Array[	<диапазон	индекса1>,	<диапазон	индекса2>,...
<диапазон индекса N> ] Of <тип компонент>;

Отдельный элемент именуется так:

<Имя массива>[<Индекс 1>,<Индекс 2>,...,<Индекс N>]

§ 3. Пример приведения матрицы к ступенчатому виду методом Гаусса на языке С.

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

Напомним, что матрица A с элементами aij называется ступенчатой, если она обладает следующими двумя свойствами:

1.если в матрице есть нулевая строка, то все строки ниже нее также нулевые;

2.пусть aij не равное 0 -- первый ненулевой элемент в строке с индексом i, т.е. элементы ail = 0 при l < j. Тогда все элементы в j-м столбце ниже элемента aij равны нулю, и все элементы левее и ниже aij также равны

нулю: akl = 0 при k > i и l =< j.

Ступенчатая матрица выглядит примерно так: !!!!СДЕЛАТЬ СТУПЕНЧАТУЮ матрицу формулой

Здесь тёмными квадратиками отмечены первые ненулевые элементы строк матрицы. Белым цветом изображаются нулевые элементы, серым цветом - произвольные элементы.

Алгоритм Гаусса использует элементарные преобразования матрицы двух типов.

•Преобразование первого рода: две строки матрицы меняются местами, и при этом знаки всех элементов одной из строк изменяются на противоположные.

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

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

Метод Гаусса в математическом варианте состоит в следующем:

1. ищем сначала ненулевой элемент в первом столбце. Если все элементы первого столбца нулевые, то переходим ко второму столбцу, и так далее. Если нашли ненулевой элемент в k-й строке, то при помощи элементарного преобразования первого рода меняем местами первую и k-ю строки, добиваясь того, чтобы первый элемент первой строки был отличен от нуля; 2. используя элементарные преобразования второго рода, обнуляем все элементы первого столбца, начиная

со второго элемента. Для этого от строки с номером k вычитаем первую строку, умноженную на коэффициент

ak1/a11 .

3. переходим ко второму столбцу (или j-му, если все элементы первого столбца были нулевыми), и в дальнейшем рассматриваем только часть матрицы, начиная со второй строки и ниже. Снова повторяем пункты 1) и 2) до тех пор, пока не приведем матрицу к ступенчатому виду.

Программистский вариант метода Гаусса имеет три отличия от математического:

1.индексы строк и столбцов матрицы начинаются с нуля, а не с единицы;

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

операцию проверки равенства вещественных чисел. Поэтому вместо проверки на равенство нулю числа aij следует сравнивать его абсолютную величину aij с очень маленьким числом ε (например, ε = 0.00000001). Если aij =< ε, то следует считать элемент aij нулевым;

3.при обнулении элементов j-го столбца, начиная со строки i + 1, мы к k-й строке, где k > i, прибавляем i-ю

строку, умноженную на коэффициент r = -akj/aij : r = -akj/aij.

ak = ak + r * ai

Такая схема работает нормально только тогда, когда коэффициент r по абсолютной величине не превосходит

единицы. В противном случае, ошибки округления умножаются на большой коэффициент и, таким образом, экспоненциально растут. Математики называют это явление неустойчивостью вычислительной схемы. Если вычислительная схема неустойчива, то полученные с ее помощью результаты не имеют никакого отношения к исходной задаче. В нашем случае схема устойчива, когда коэффициент r = -akj/aij не превосходит по модулю единицы. Для этого должно выполняться неравенство

aij >= akj при k > i

Отсюда следует, что при поиске разрешающего элемента в j-м столбце необходимо найти не первый попавшийся ненулевой элемент, а максимальный по абсолютной величине. Если он по модулю не превосходит ε, то считаем, что все элементы столбца нулевые; иначе меняем местами строки, ставя его на вершину столбца, и затем обнуляем столбец элементарными преобразованиями второго рода.

Ниже дан полный текст программы на Си, приводящей вещественную матрицу к ступенчатому виду. Функция, реализующая метод Гаусса, одновременно подсчитывает и ранг матрицы. Программа вводит размеры матрицы и ее элементы с клавиатуры и вызывает функцию приведения к ступенчатому виду. Затем программа печатает ступенчатый вид матрицы и ее ранг. В случае квадратной матрицы также вычисляется и печатается определитель матрицы, равный произведению диагональных элементов ступенчатой матрицы.

При реализации метода Гаусса используется схема построения цикла с помощью инварианта,. В цикле меняются две переменные -- индекс строки i, 0 =< i < m - 1, и индекс столбца j, 0 =< j < n - 1. Инвариантом цикла является утверждение о том, что часть матрицы (математики говорят минор) в столбцах 0,1,...j - 1 приведена к ступенчатому виду и что первый ненулевой элемент в строке i - 1 стоит в столбце с индексом меньшим j. В теле цикла рассматривается только минор матрицы в строках i,...,m - 1 и столбцах j,...,n - 1. Сначала ищется максимальный по модулю элемент в j-м столбце. Если он по абсолютной величине не превосходит ε, то j увеличивается на единицу (считается, что столбец нулевой). Иначе перестановкой строк разрешающий элемент ставится на вершину j-го столбца минора, и затем столбец обнуляется элементарными преобразованиями второго рода. После этого оба индекса i и j увеличиваются на единицу. Алгоритм завершается, когда либо i = m, либо j = n. По окончании алгоритма значение переменной i равно числу ненулевых строк ступенчатой матрицы, т.е. рангу исходной матрицы.

Для вычисления абсолютной величины вещественного числа x типа double мы пользуемся стандарной математической функцией fabs(x), описанной в стандартном заголовочном файле "math.h.

#include <stdio.h> // Описания функций ввода-вывода #include <math.h> // Описания математических функций

#include <stdlib.h> // Описания функций malloc и free

//Прототип функции приведения матрицы

//к ступенчатому виду.

//Функция возвращает ранг матрицы int gaussMethod(

int m,	// Число строк матрицы
int n,	// Число столбцов матрицы
double *a,	//	Адрес массива элементов матрицы
double eps	//	Точность вычислений
);

int main() {

int m, n, i, j, rank; double *a;

double eps, det;

printf("Введите размеры матрицы m, n: ");

scanf("%d%d", &m, &n);

// Захватываем память под элементы матрицы

a = (double *) malloc(m * n * sizeof(double));

printf("Введите элементы матрицы:\n"); for (i = 0; i < m; ++i) {

for (j = 0; j < n; ++j) {

// Вводим элемент с индексами i, j scanf("%lf", (a+i*n + j));

}

}

printf("Введите точность вычислений eps: "); scanf("%lf", &eps);

// Вызываем метод Гаусса

rank = gaussMethod(m, n, a, eps);

// Печатаем ступенчатую матрицу printf("Ступенчатый вид матрицы:\n"); for (i = 0; i < m; ++i) {

// Печатаем i-ю строку матрицы for (j = 0; j < n; ++j) {

printf(	// Формат %10.3lf означает 10
"%10.3lf ", // позиций на печать числа,
*(a+i*n + j) // 3 знака после точки
);
}	// Перевести строку
printf("\n");

}

// Печатаем ранг матрицы printf("Ранг матрицы = %d\n", rank);

if (m == n) {

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

//ее определитель

det = 1.0;

for (i = 0; i < m; ++i) { det *= *(a+i*n + i);

}

printf("Определитель матрицы = %.3lf\n", det);

}
free(a);	// Освобождаем память
return 0;	// Успешное завершение программы

}

//Приведение вещественной матрицы

//к ступенчатому виду методом Гаусса с выбором

//максимального разрешающего элемента в столбце.

//Функция возвращает ранг матрицы

int gaussMethod(	// Число строк матрицы
int m,
int n,	// Число столбцов матрицы
double *a,	// Адрес массива элементов матрицы
double eps	// Точность вычислений
) {
int i, j, k, l;
double r;

i = 0; j = 0;

while (i < m && j < n) {

//Инвариант: минор матрицы в столбцах 0..j-1

//уже приведен к ступенчатому виду, и строка

//с индексом i-1 содержит ненулевой эл-т

//в столбце с номером, меньшим чем j

//Ищем максимальный элемент в j-м столбце,

//начиная с i-й строки

r =	0.0;
for (k = i; k < m; ++k) {
	if (fabs(a[k*n + j]) > r) {
	l = k;	// Запомним номер строки

r = fabs(*(a+k*n + j)); // и макс. эл-т

}

}

if (r <= eps) {

//Все элементы j-го столбца по абсолютной

//величине не превосходят eps.

//Обнулим столбец, начиная с i-й строки for (k = i; k < m; ++k) {

*(a+k*n + j) = 0.0;

}

++j;	//	Увеличим индекс столбца
continue; //		Переходим к следующей итерации

}

if (l != i) {

// Меняем местами i-ю и l-ю строки for (k = j; k < n; ++k) {

r = a[i*n + k];

*(a+i*n + k) = *(a+l*n + k);

*(a+l*n + k) = (-r); // Меняем знак строки

}

}

//Утверждение: fabs(*(a+i*n + k)) > eps

//Обнуляем j-й столбец, начиная со строки i+1,

//применяя элем. преобразования второго рода for (k = i+1; k < m; ++k) {

r = (-(*(a+k*n + j)) / *(a+i*n + j));

//К k-й строке прибавляем i-ю, умноженную на r

*(a+k*n + j) = 0.0;

for (l = j+1; l < n; ++l) {

*(a+k*n + l) += r * (*(a+i*n + l));

}

}

++i; ++j; // Переходим к следующему минору

}

return i; // Возвращаем число ненулевых строк

}

Приведем два примера работы этой программы. В первом случае вводится вырожденная матрица размера 4 × 4:

Введите размеры матрицы n, m: 4 4 Введите элементы матрицы:

1 2 3 4

4 3 2 1

5 6 7 8

8 7 6 5