Алгоритмы C++

for (int i=0; i<n; ++i)

printf ("%.15lf\n", a[i][n]);

Применения

Помимо непосредственно решения систем, метод Гаусса находит применение и в других задачах:

●Нахождение определителя матрицы

●Нахождение ранга матрицы

Нахождение ранга матрицы

Ранг матрицы - это наибольшее число линейно независимых строк/столбцов матрицы. Ранг определён не только для квадратных матриц; пусть матрица прямоугольна и имеет размер NxM.

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

Заметим, что если матрица квадратная и её определитель отличен от нуля, то ранг равен N(=M), иначе он будет меньше. В общем случае, ранг матрицы не превосходит min(N,M).

Алгоритм

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

операции, что и при решении системы или нахождении её определителя, но если на каком-либо шаге в i-ом столбце среди невыбранных до этого строк нет ненулевых, то мы этот шаг пропускаем, а ранг уменьшаем на единицу (изначально ранг полагаем равным max(N,M)). Иначе, если мы нашли на i-ом шаге строку с ненулевым элементом в i-

ом столбце, то помечаем эту строку как выбранную, и выполняем обычные операции отнимания этой строки от остальных.

Реализация

const double EPS = 1E-9;

int rank = max(n,m); vector<char> line_used (n); for (int i=0; i<m; ++i) {

int j;

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

if (!line_used[j] && abs(a[j][i]) > EPS) break;

if (j == n) --rank;

else {

line_used[j] = true;

for (int p=i+1; p<m; ++p)

a[j][p] /= a[j][i]; for (int k=0; k<n; ++k)

if (k != j && abs (a[k][i]) > EPS) for (int p=i+1; p<m; ++p)

a[k][p] -= a[j][p] * a[k][i];

}

}

Вычисление определителя матрицы методом Гаусса

Пусть дана квадратная матрица A размером NxN. Требуется вычислить её определитель.

Алгоритм

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

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

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

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

Таким образом, мы можем с помощью алгоритма Гаусса вычислять определитель матрицы за O (N3).

Осталось только заметить, что если в какой-то момент мы не найдём в текущем столбце ненулевого элемента, то алгоритм следует остановить и вернуть 0.

Реализация

const double EPS = 1E-9;
int n;	> a (n, vector<double> (n));
vector < vector<double>
... чтение n и a ...
double det = 1;	{
for (int i=0; i<n; ++i)
int k = i;	j<n; ++j)
for (int j=i+1;
if (abs	(a[j][i]) > abs (a[k][i]))
	k = j;

if (abs (a[k][i]) < EPS) { det = 0;

break;

}

swap (a[i], a[k]); if (i != k)

det = -det; det *= a[i][i];

for (int j=i+1; j<=n; ++j) a[i][j] /= a[i][i];

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

if (j != i && abs (a[j][i]) > EPS) for (int k=i+1; k<n; ++k)

a[j][k] -= a[i][k] * a[j][i];

}

cout << det;

Интегрирование по формуле Симпсона

Требуется посчитать значение определённого интеграла:

Решение, описываемое здесь, было опубликовано в одной из диссертаций Томаса Симпсона (Thomas Simpson)

в 1743 г.

Формула Симпсона

Пусть — некоторое натуральное число. Разобьём отрезок интегрирования на равных частей:

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

Её можно найти аналитически, затем останется только проинтегрировать выражение для неё, и окончательно получаем:

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

Погрешность

Погрешность, даваемая формулой Симпсона, не превосходит по модулю величины:

Таким образом, погрешность имеет порядок уменьшения как .

Реализация

Здесь — некоторая пользовательская функция.

double a, b; // входные данные

const int N = 1000*1000; // количество шагов (уже умноженное на 2) double s = 0;

double h = (b - a) / N; for (int i=0; i<=N; ++i) {

double x = a + h * i;

s += f(x) * ((i==0 || i==N) ? 1 : ((i&1)==0) ? 2 : 4);

}

s *= h / 3;

Метод Ньютона (касательных) для поиска корней

Это итерационный метод, изобретённый Исааком Ньютоном (Isaak Newton) около 1664 г. Впрочем, иногда этот метод называют методом Ньютона-Рафсона (Raphson), поскольку Рафсон изобрёл тот же самый метод на несколько лет позже Ньютона, однако его статья была опубликована намного раньше.

Задача заключается в следующем. Дано уравнение:

F(X) = 0

Требуется решить это уравнение, точнее, найти один из его корней (предполагается, что корень существует). Предполагается, что F(X) непрерывна и дифференцируема на отрезке [A;B].

Алгоритм

Входным параметром алгоритма, кроме функции F(X), является также начальное приближение - некоторое X0, от которого алгоритм начинает идти.

Пусть уже вычислено Xi, вычислим Xi+1 следующим образом. Проведём касательную к графику функции F(X) в точке X = Xi, и найдём точку пересечения этой касательной с осью абсцисс. Xi+1 положим равным найденной точке, и повторим весь процесс с начала.

Нетрудно получить следующее выражение:

Xi+1 = Xi - F(Xi) / F'(Xi)

Интуитивно ясно, что если функция F(X) достаточно "хорошая", а Xi находится достаточно близко от корня, то Xi+1 будет находиться ещё ближе к искомому корню.

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

Применение для вычисления квадратного корня

Рассмотрим метод Ньютона на примере вычисления квадратного корня. Если подставить F(X) = sqrt(X), то после упрощения выражения получаем:

Xi+1 = (Xi + N / Xi) / 2

Первый типичный вариант задачи - когда дано дробное число N, и нужно подсчитать его корень с некоторой точностью EPS:

	double	n;
	cin >>	n;
	const double EPS = 1E-15;
	double	x = 1;
	for (;;) {
		double nx = (x + n / x) / 2;
		if (abs (x - nx) < EPS) break;
	}	x = nx;
		("%.15lf", x);
	printf

Другой распространённый вариант задачи - когда требуется посчитать целочисленный корень (для данного N

найти наибольшее X такое, что X2<=N). Здесь приходится немного изменять условие останова алгоритма, поскольку может случиться, что X начнёт "скакать" возле ответа. Поэтому мы добавляем условие, что если значение X

на предыдущем шаге уменьшилось, а на текущем шаге пытается увеличиться, то алгоритм надо остановить.

int n; cin >> n; int x = 1;

bool decreased = false; for (;;) {

int nx = (x + n / x) >> 1;

if (x == nx || nx > x && decreased) break; decreased = nx < x;

x = nx;

}

cout << x;

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

начального приближения число 2bits/2, где bits - количество битов в числе N. Вот код на языке Java, демонстрирующий этот вариант:

BigInteger n; // входные данные

BigInteger a = BigInteger.ONE.shiftLeft (n.bitLength() / 2); boolean p_dec = false;

for (;;) {

BigInteger b = n.divide(a).add(a).shiftRight(1);

if (a.compareTo(b) == 0 || a.compareTo(b) < 0 && p_dec) break; p_dec = a.compareTo(b) > 0;

a = b;

}

Например, этот вариант кода выполняется для числа 101000 за 60 миллисекунд, а если убрать улучшенный выбор начального приближения (просто начинать с 1), то будет выполняться примерно 120 миллисекунд.