3 Глава Параллельные вычисления при численном решении уравнений математической физики

LU разложение матриц

Рассмотрим систему Ax = f . В результате применения метода Гаус- са мы получили систему Ux = y, где U есть верхняя треугольная матрица с единицами на главной диагонали. Вектор y зависит толь- ко от вектора f и матрицы A. Рассмотрим соотношение между y и f . В случае n = 3 мы имеем

y = (f − ⁽¹⁾ya )/a ,_{2 2 1 21}

^y₃ = (f₃ − y₁a₃₁ − y₂a⁽¹⁾)/a⁽²⁾.
^_y₁ = f₁/a₁₁,


22

Отсюда следует

32 33



f₁ = a₁₁y₁,

22
f₂ = y₁a₂₁ + y₂a⁽¹⁾,

32

33
^f₃ = y₁a₃₁ + y₂a⁽¹⁾ + y₃a⁽²⁾.

Понятно, что для произвольных n мы будем иметь

^  

a

( ⁻a , i ≥ j,

y₃ . (16)

y_n
Таким образом, f = Ly, где L – нижняя треугольная матрица с элементами

l_ij =
(j 1)
ij
0, i < j,

которая может быть рассчитана с использованием формул базового алгоритма (10)-(12). Так как

то получаем

Ax = f, Ux = y, Ly = f,

Ax = Ly = LUx.

Следовательно, мы имеем A = LU , где L есть нижняя треугольная матрица и U - верхняя треугольная матрица. Это представление LU разложения для матрицы A часто называется LU факторизацией матрицы A. Приведем пример для 3 × 3 матрицы

20 15 12

20 5 1/3

0 0 1

30 20 15^ = ^30 5 0

^{ }0 1 1

 _.
60 30 20 60 0 0  1 1/2 1/3

Предположим, что для системы Ax = f мы знаем матрицы L и U . Тогда мы сможем легко решить систему в два простых шага: первый - Ly = f дает решение для y, второй - Ux = y обеспечи- вает решение для x. Проблема состоит в получении матриц L и U . Матрица A имеет n² коэффициентов, тогда число коэффициентов в двух матрицах L и U определяется как n(n + 1), поэтому есть неко- торая недоопределенность. Это дает нам свободу выбора некоторых коэффициентов в L или U .
Укажем некоторые конкретные варианты:

•
Если метод исключения Гаусса косвенно выбирает l_ii = 1, то этот вариант называется Метод Дулитла.

•
Представление с условием u_ii = 1 известно как Метод Кроута.

•
Если A = A^T (т. е. A является симметричной матрицей) и A является также положительно-определенной матрицей (т. е. имеет положительные собственные значения), то мы можем выбрать u_ij = l_ji. Этот вариант называется Метод Холецкого. В этом случае будем иметь

A = LU, A^T = U^T L^T = LU, → L = U^T .
Заметим, что, строго говоря, LU разложение матрицы не все- гда возможно. Теорема утверждает, что если все главные миноры матрицы A не равны нулю, то A представляется как A = LU един- ственным образом. При этом диагональные элементы L отличны от нуля.
С другой стороны, все гауссовские преобразования могут быть записаны в матричной форме. Например, если n = 3, имеем
U = L₃L₂L₁A,

6. Метод гауссовых исключений с перестанов- кой строк

k+1,k+1
Что делать, если на одном шаге ведущий элемент a^(k) = 0?
Теорема утверждает, что если det A 0, существует матрица пере-
становок P , такая что матрица PA имеет все основные миноры, не равные нулю, следовательно PA = LU .

(k)

k+1,k+1
Основная идея исключения по Гауссу с перестановкой стро- кой заключается в том, что на каждом шаге процедуры исклю- чения неизвестных мы осуществляем перестановку строк с i, i = k + 1, ..., n, таким образом, что мы имеем старший элемент a^(k) максимальный по модулю. Мы ищем строку i с максимальным

| |
a_i,k+1 , i = k + 1, ..., n (см. (9)), а затем осуществляем перестановку между строками k + 1 и i. Наконец, получим
^{L = L}n^Ln−1^Pn−1,j_n−1^...L1^P1,j₁^A,

−
где P_k,jk с k = 1, ..., n 1, является матрицей перестановки строк. Например, если n = 3,
0 1 0 1 0 0 0 0 1




0 0 1

0 1 0

1 0 0
P₁₂ = ^1 0 0^ , P₂₃ = ^0 0 1^ , P₁₃ = ^0 1 0^ ,

обеспечивают перестановку между строками 1 и 2, 2 и 3, 1 и 3.

7. Численный расчет матрицы определителя и обратной матрицы

Используя LU разложение, легко вычислить определитель матри- цы
det A = ± det(PA) = ± det(LU ) = ± det L det U
= ± det L = ±l₁₁l₂₂...l_nn.
Если A⁻¹ является матрицей, обратной по отношению к A, то- гда AA⁻¹ = I, где I - единичная матрица. Обозначим X = A⁻¹. На самом деле мы имеем n² систем линейных алгебраических уравне-
ний

Σ
n
^aik^xkj ^{= δ}ij^,
k=1
или Ax^(j) = δ^(j), где x^(j) = (x_1j, x_2j, ..., , x_nj)^T и δ^(j) = (δ_1j, ..., δ_nj)^T ,
i, j = 1, ..., n. С помощью разложения A = LU мы должны сначала

решить систему Ly^(j) = δ^(j) для y^(j) = (y_1j, y_2j, ..., , y_nj)^T , и затем
_Ux(j) _{= y}(j) _{для x}(j)_.

(17)

l_ji =

_li−_i 1

a_ij −
i−1 k=1
^lik^ljk!



Σ
j = i + 1, i + 2, . . . n. (18)

Мы отмечаем, что (17) предполагает использование квадратно- го корня. Следовательно, факторизация будет существовать только для ограниченного класса матриц: это симметричные положитель- но определенные матрицы. Симметричная матрица A положитель- но определена, если выполняется условие (x, Ax) > 0 для всех нену- левых векторов x. Необходимое и достаточное условие для положи- тельной определенности состоит в том, что все собственные значе- ния положительны. Можно проверить матрицу на положительную определенность, пытаясь сделать вычисление факторизацией Хо- лецкого. На самом деле трудным шагом здесь является знание того, что матрица положительно определенная! Это предполагает поиск собственных значений и может занять больше времени, чем вре- мя на решение уравнения!

Глава 2. Методы Рунге-Кутты и другие методы

1.1 Методы Рунге-Кутты

Широкая категория методов, наиболее часто применяемых на практике для решения дифференциальных уравнений, известна под общим названием методы Рунге-Кутты. Различные методы этой категории требуют большего или меньшего объема вычислений и соответственно обеспечивают большую или меньшую точность.
Методы Рунге-Кутты обладают следующими отличительными свойствами:
• эти методы являются одноступенчатыми: чтобы найти значение
функции в точке у,+| нужна информация только о предыдущей точке (у„ х,);

• 
  • они согласуются с рядом Тейлора вплоть до членов порядка h^k, где степень к определяет порядок метода;
  
• 
  • эти методы не требуют вычисления производных отДх,у), а требуют вычисления самой функции.
  

Именно благодаря последнему свойству методы Рунге-Кутты более удобны для практических вычислений.

Метод Эйлера (метод Рунге-Кутты первого порядка)

Простейшим из численных методов решения дифференциальных уравнений является метод Этери. Это один из самых старых и широко известных методов. Метод Эйлера является сравнительно грубым методом решения дифференциальных уравнений, однако идеи, положенные в его основу, являются, по существу, исходными для очень широкого класса численных методов.
Пусть требуется найти приближенное решение дифференциального уравнения первого порядка 
с начальным условием 
т. е. необходимо решить задачу Коши.
В окрестности точки х₍₎ функциюу(х) разложим в ряд Тейлора:

который можно применить для приближенного определения искомой функции у(х). В точке х₍₎ + h при малых значениях h можно ограничиться двумя членами ряда (8.7), тогда

где 0(1г²) - бесконечно малая величина порядка /г². Заменим производную y'(x_Q), входящую в формулу (8.7), на правую часть уравнения (8.5):

Теперь приближенное решение в точке х_{ =x_Q + h можно вновь рассматривать как начальное условие и по формуле (8.9) найти значение искомой функции в следующей точке х₀=х_{ + h. В результате получен простейший алгоритм решения задачи Коши, который называется методом Эйлера, или методом ломаных.
Метод Эйлера можно представить в виде последовательного применения формул:
для точки 
Таким образом, формула Эйлера в общем случае имеет вид:

Название «метод ломаных» связано с его геометрической интерпретацией. Искомая функция >{л') заменяется ломаной линией, представляющей собой отрезки касательных к этой функции в узлах х₍₎, х_п. Выведем формулы на основе геометрических аналогий. Предположим, что нам известна точка (x₀, у₀) на искомой интегральной кривой (см. рис. 8.1).

Рис. 8.1
Через точку (х«, у₀) проведем касательную с тангенсом угла наклона:   Уравнение касательной имеет вид: 
Тогда в точке X] =л'₍₎ + h, с учетом (8.13), получим решение:

Ошибка решения в точке x=x_t показана в виде отрезка Д.
Формула (8.12) является методом Рунге-Кутты первого порядка, т. к. она согласуется с разложением в ряд Тейлора вплоть до членов порядка h^l.
Метод Эйлера имеет довольно большую погрешность вычисления: Д«0(И). Кроме того, он очень часто оказывается неустойчивым - малая ошибка (например, заложенная в исходных данных) увеличивается с ростом х. На рис. 8.2 приведена блок-схема метода Эйлера, в приложении - программа расчета.

Рис. 8.2. Блок-схема метода Эшера
Пример 8.1. Для химической реакции

Изменение концентраций веществ А и В можно описать следующими кинетическими уравнениями:

с начальными условиями:

Требуется получить зависимость изменения концентрации вещества А от времени, т. е. необходимо решить дифференциальное уравнение (решить задачу Коши).
Исходные данные:
G, о=1 моль/л;
Св. о⁼ 0;
к=0,2 с интервал интегрирования /=[0, 5].
Обозначим С_А =у, тогда

Примем величину шага h = 0,1. Решим данное уравнение методом Эйлера (8.12).