Прямые методы решения СЛАУ.#
Помимо известного вам метода Гаусса и метода Жордана-Гаусса (с выбором главного элемента) рассмотрим ещё пару методов.
Метод прогонки используется для матриц с небольшим числом диагоналей, а метод вращений обладает большей устойчивостью, чем метод Гаусса.
Метод прогонки#
Метод прогонки (TDMA)
Рассмотрим трехдиагональную матрицу \(A\) и систему \(Ax = f\). Тогда, определяя \(\alpha_k, \beta_k\) следующим реккурентным образом:
Находим \(x_i\), используя найденные \(\alpha_k, \beta_k\):
Доказательство
\(\Delta\) Общая форма:
Представим решение в следующем реккурентном виде:
Подставляя в общую форму, имеем:
Учитывая независимость \(\alpha_i, \beta_i\) и коэффициентов матрицы от \(x_k\), можем продиффиринцировать обе части уравнения по \(x_{i+1}\):
Получаем рекуррентные формулы, по которым находим \(\alpha_k, \, \beta_k\):
Идя по убыванию индексов, находим \(x_i\):
Show code cell source
import numpy as np
import time
import matplotlib.pyplot as plt
from tqdm import tqdm
def TDMA(A, B, C, F): # Прогонка
# A, B, C, F - просто числовые последовательности, НЕ матрицы
A = np.insert(A, 0, 0)
n = len(F)
x = np.zeros(n)
alpha = np.zeros(n)
beta = np.zeros(n)
alpha[1] = -C[0]/B[0]
beta[1] = F[0]/B[0]
for i in range(2, n):
alpha[i] = -C[i-1]/(A[i-1]*alpha[i-1] + B[i-1])
beta[i] = (F[i-1] - A[i-1]*beta[i-1])/(A[i-1]*alpha[i-1] + B[i-1])
x[n-1] = (F[n-1] - A[n-1]*beta[n-1])/(B[n-1] + A[n-1]*alpha[n-1])
for i in range(n-1, 0, -1): # помним, что последнее значение не кушает
x[i - 1] = alpha[i]*x[i] + beta[i]
return x
def gaussian_elimination(matrix, rhs): # Gaussian Elimination
n = len(rhs)
for i in range(n):
max_row = i + np.argmax(np.abs(matrix[i:, i]))
matrix[[i, max_row]] = matrix[[max_row, i]]
rhs[[i, max_row]] = rhs[[max_row, i]]
pivot = matrix[i, i]
matrix[i] = matrix[i] / pivot
rhs[i] = rhs[i] / pivot
for j in range(i+1, n):
factor = matrix[j, i]
matrix[j] -= factor * matrix[i]
rhs[j] -= factor * rhs[i]
x = np.zeros(n)
for i in range(n-1, -1, -1):
x[i] = rhs[i] - np.sum(matrix[i, i+1:] * x[i+1:])
return x
# Sizes and sample size
sizes = range(10, 100, 10)
sample_size = 1000
time_results_tdma = []
accuracy_results_tdma = []
time_results_gaussian = []
accuracy_results_gaussian = []
for size in tqdm(sizes, desc="Processing sizes"):
total_time_tdma = 0
total_error_tdma = 0
total_time_gaussian = 0
total_error_gaussian = 0
for _ in range(sample_size):
A = np.random.rand(size-1)
B = np.random.rand(size)
C = np.random.rand(size-1)
F = np.random.rand(size)
exact_x = np.random.rand(size)
F = np.zeros(size)
F[0] = B[0] * exact_x[0] + C[0] * exact_x[1]
for i in range(1, size-1):
F[i] = A[i-1] * exact_x[i-1] + B[i] * exact_x[i] + C[i] * exact_x[i+1]
F[-1] = A[-1] * exact_x[-2] + B[-1] * exact_x[-1]
start_time = time.time()
computed_x_tdma = TDMA(A, B, C, F)
end_time = time.time()
total_time_tdma += (end_time - start_time)
total_error_tdma += np.linalg.norm(computed_x_tdma - exact_x)
matrix = np.zeros((size, size))
np.fill_diagonal(matrix, B)
np.fill_diagonal(matrix[1:], A)
np.fill_diagonal(matrix[:, 1:], C)
start_time = time.time()
computed_x_gaussian = gaussian_elimination(matrix.copy(), F.copy())
end_time = time.time()
total_time_gaussian += (end_time - start_time)
total_error_gaussian += np.linalg.norm(computed_x_gaussian - exact_x)
time_results_tdma.append(total_time_tdma / sample_size)
accuracy_results_tdma.append(total_error_tdma / sample_size)
time_results_gaussian.append(total_time_gaussian / sample_size)
accuracy_results_gaussian.append(total_error_gaussian / sample_size)
# Plotting results
fig, ax1 = plt.subplots()
ax1.set_xlabel('Matrix Size (n)')
ax1.set_ylabel('Time (s)', color='tab:blue')
ax1.plot(sizes, time_results_tdma, label="TDMA Time", color='tab:red')
ax1.plot(sizes, time_results_gaussian, label="Gaussian Time", color='tab:cyan')
ax1.set_yscale("log")
ax1.tick_params(axis='y', labelcolor='tab:blue')
ax2 = ax1.twinx()
ax2.set_ylabel('Accuracy (Mean Error)', color='tab:red')
ax2.plot(sizes, accuracy_results_tdma, label="TDMA Accuracy", color='tab:red', linestyle='--')
ax2.plot(sizes, accuracy_results_gaussian, label="Gaussian Accuracy", color='tab:cyan', linestyle='--')
ax2.tick_params(axis='y', labelcolor='tab:red')
ax1.legend(loc='upper left')
ax2.legend(loc='upper right')
fig.tight_layout()
plt.title('Comparison of TDMA and Gaussian Elimination')
plt.show()
Processing sizes: 100%|███████████████████████████████████████████████████████████████████| 9/9 [00:32<00:00, 3.66s/it]
Метод вращений#
Коэффициенты матрицы при выполнении метода Гаусса могут очень быстро расти. Этого недостатка лишён метод вращений.
Как и метод Гаусса, метод вращений нужен для приведения системы к треугольному виду с последующим решением.
Выбираются 2 неотрицательных числа \(s_1, c_1\). Обозначим первую строку \(a_1\), а вторую \(a_2\). Заменим первую строку на строку \(c_1a_1 + s_1a_2\), а вторую на строку \(-s_1a_1 + c_1a_2\).
На числа \(c_1\) и \(s_1\) накладываются условия:
\(-s_1a_{11} + c_1a_{21} = 0\)
\(s_1^2 + c_2^2 = 1\)
Откуда \(c_1 = \frac{a_{11}}{\sqrt{a_{11}^2 + a_{21}^2}}\) и \(s_1 = \frac{a_{21}}{\sqrt{a_{11}^2 + a_{21}^2}}\)
Данные числа являются синусом и косинусом некоторого угла.
Сначала зануляем первый столбец, потом переходим ко второму и так далее.
Обратный ход абсолютно аналогичен обратному ходу в методе Гаусса.
Уменьшение невязки повторным решением.#
При решении прямыми методами мы получали некоторое решение \(x\) с погрешностью. Введем понятие невязки:
Рассмотрим точное решение \(x^*\). Тогда погрешность решения \(e = x^* - x\).
Запишем получившуюся систему:
, тогда \(A(x - x^*) = r\).
Перепишем систему :
Решая эту систему относительно \(e\), получим новое решение \(\hat x\), которое будет ближе к точному. Таким образом можно достигать заданного уровня точности.