Метод на прогонката за решаване на системи линейни уравнения с тридиагонална матрица. Числено решаване на гранична задача за ОДУ

Дефиниция (матрица с доминиращ главен диагонал)

Казваме, че квадратната матрицата $A = (a_{ij})_{nxn}$ е с доминиращ главен диагонал по редове (по стълбове) ако е изпълнено, че:

(1)
\begin{eqnarray} & & | a_{i,i}| > \sum_{j = 1, j \neq i}^n |a_{i,j}|, \ i = 1, 2, \dots, n \\ & & \bigg(| a_{j,j}| > \sum_{i = 1, i \neq j}^n |a_{i,j}|, \ j = 1, 2, \dots, n \bigg) \end{eqnarray}

Теорема

Всяка матрица с доминиращ главен диагонал е неособена. (т.е. $\mathrm{det} A \neq 0$ )

Метод на прогонката

Всъщност този метод е най-логичното нещо което се прави когато трябва да се реши система от следния тип:

(5)
\begin{array} {|ccccccccccl} b_1x_1 &+& c_1x_2 &&&&&& &=& d_1 \\ a_2x_1 &+& b_2x_2 &+& c_2x_3 &&&& &=& d_2 \\ &\ddots & & \ddots & & \ddots & & & & & \vdots \\ & &a_kx_{k-1} &+&b_kx_k &+& c_kx_{k+1} && &=& d_k \\ & & &\ddots& & \ddots & & \ddots & & & \vdots \\ &&&& a_{n-1}x_{n-2} &+& b_{n-1}x_{n-1} &+& c_{n-1}x_n &=& d_{n-1} \\ &&&&&& a_nx_{n-1} &+& b_nx_n &=& d_n \\ \end{array}

Ако от първото уравнение, в което има две неизвестни, изразим едната и заместим във второто, тогава и второто ще стане уравнение с две неизвестни.. и пак изразяваме, заместваме.. и т.н. Тогава от предпоследното уравнение изразяваме предпоследната неизвестна, заместваме в последното уравнение и съответно то се получава линейно уравнение на една неизвестна (което е решимо). И после, в обратен ред намираме всички други неизвестни. Но стига думи, ето малко код:

Прав ход на прогонката

$x_k = \alpha_k x_{k+1} + \beta_k, \ k = 1, 2, \dots, n - 1$
$\displaystyle{\alpha_1 = - \frac{c_1}{b_1}, \ \beta_1 = \frac{d_1}{b_1}}$

Ако сме изразили $\alpha_{k-1}$ и $\beta_{k-1}$, тогава $x_{k-1} = \alpha_{k-1}x_k + \beta_{k-1}$
Заместваме $x_{k-1}$ в $k$-тото уравнение и получаваме
$a_k(\alpha_{k-1}x_k + \beta_{k-1}) + b_kx_k + c_kx_{k+1} = d_k$
$\Longrightarrow (a_k\alpha_{k-1} + b_k) x_k = -c_k x_{k+1} + d_k - a_k\beta_{k-1}$
от където $\displaystyle{\alpha_k = - \frac{c_k}{a_k\alpha_{k-1} + b_k}}$ и $\displaystyle{\beta_k = \frac{d_k - a_k\beta_{k-1}}{a_k\alpha_{k-1} + b_k}}$
Всичко това става за $k = 2, 3, \dots, n-1$

Тогава при $k = n-1$
$x_{n-1} = \alpha_{n-1}x_n + \beta_{n-1}$
И заместваме в последното уравнение
$a_n(\alpha_{n-1}x_n + \beta_{n-1}) + b_nx_n = d_n$
$\Longrightarrow (a_n\alpha_{n-1})x_n = d_n - a_n\beta_{n-1}$
$\displaystyle{\Longrightarrow x_n = \frac{d_n - a_n\beta_{n-1}}{a_n\alpha_{n-1} + b_n} = \beta_n}$ (т.е формулата за $x_n$ се оказа същата като за $\beta_n$).

Обратен ход на прогонката
С обратния ход просто намираме всички неизвестни в обратен ред (разбира се, след като сме намерили послената неизвестна променлива)
$x_k = \alpha_kx_{k+1} + \beta_k, \ k =n-1, n-2, \dots, 1$

Забележка: Ако на някое място се наложи да делим на нула, именно $a_k \alpha_{k-1} + b_k = 0$ то намираме $x_{k+1} = \frac{d_k - a_k \beta_{k-1}}{c_k}$ и като заместим в долното уравнение се получават две триагонални матрици (един вид - разцепили сме на $k+1$вия ред - отгоре и отдолу остават тридиагонални и пускаме рекурсивно за 2те части - по точно - за след $k + 1$ пускаме първи и втори ход на прогонката, а преди $k + 1$ пускаме обратен ход).

Използване на кубичен сплайн за числено решаване на гранична задача за ОДУ от 2-ри ред

Нека имаме следната гранична задача за ОДУ1 от 2-ри ред

(6)
\begin{array} {|cccc} y'' - q(x)y &=& f(x) \\ y(a) &=& \alpha \\ y(b) &=& \beta \\ \end{array}

Тук $x \in [a,b]$, а функциите $f(x)$ и $g(x)$ са дадени.

Ще намерим приближено решение $y(x)$ на системата с помощта на кубичен сплайн.
Първо - харесваме си $n+1$ възела от интервала $[a,b]$
$a = x_0 < x_1 < \dots < x_n = b$

От сега нататък ще записваме за по-кратко ето така:
$f(x_i) = f_i, \ g(x_i) = g_i, \ y(x_i) = y_i,\ \ i = 0, 1, \dots, n$

Освен това, правим и следното означение:
$\Delta_i = x_{i+1} - x_i, \ i = 0, 1, \dots, n-1$
Т.е $\Delta_i$ е дължината на $i$-тия интервал.

Търсим кубичен сплайн $s(x)$:
$s(x)|_{(x_i, x_{i+1})} = p_i \in \Pi_3, \ i=0,1, \dots, n-1$
като ще предполагаме, че $s(x)$ удовлетворява диференциалното уравнение във възлите $\{x_i\}_{i=0}^n$, т.е

(7)
\begin{array} {|ccccc} &p''_i(x_i) - q_ip_i(x_i) &=& f_i \\ &p''_i(x_{i+1}) - q_{i+1}p_i(x_{i+1}) &=& f_{i+1} \end{array}

Да означим $p''_i(x_i) = M_i$ и $p''_i(x_{i+1}) = M_{i+1}$.
Тъй като $p_i(x) \in \Pi_3$, то $p''_i(x) \in \Pi_1$

(8)
\begin{align} p''_i(x) = \frac{x_{i+1} - x}{\Delta_i} M_i + \frac{x - x_i}{\Delta_i} M_{i+1} \end{align}

Сега, с елементарно интегриране получаваме

(9)
\begin{align} p'_i(x) = - \frac{(x - x_{i+1})^2}{2 \Delta_i} M_i + \frac{(x - x_i)^2}{2 \Delta_i} M_{i+1} + A_i \end{align}

където $A_i$ е константа.

И съответно с още едно интегриране:

(10)
\begin{align} p_i(x) = - \frac{(x - x_{i+1})^3}{6 \Delta_i} M_i + \frac{(x - x_i)^3}{6 \Delta_i} M_{i+1} + A_i(x - x_i) + B_i \end{align}

където $B_i$ е друга константа.

Константите $A_i$ и $B_i$ намираме от условията: $p_i(x_i) = y_i$ и $p_i(x_{i+1}) = y_{i+1}$

(11)
\begin{align} p_i(x_i) = \frac{\Delta_i^2}{6} M_i + B_i = y_i \Longrightarrow B_i = y_i - \frac{\Delta_i^2}{6} M_i \end{align}
(12)
\begin{align} p_i(x_{i+1}) = \frac{\Delta_i^2}{6} M_{i+1} + A_i \Delta_i + y_i - \frac{\Delta_i^2}{6} M_i = y_{i+1} \Longrightarrow A_i = \frac{y_{i+1} - y_i}{\Delta_i} - \frac{\Delta_i}{6} (M_{i+1} - M_i) \end{align}

Така получаваме

(13)
\begin{align} p_i(x) = - \frac{(x - x_{i+1})^3}{6 \Delta_i} M_i + \frac{(x - x_i)^3}{6 \Delta_i} M_{i+1} + \Big(\frac{y_{i+1} - y_i}{\Delta_i} - \frac{\Delta_i}{6} (M_{i+1} - M_i)\Big)(x - x_i) + y_i - \frac{\Delta^2}{6} M_i \end{align}

От условието за непрекъснатост знаем, че $p'_{i-1}(x_i) = p'_i(x_i), \ i = 1, 2, \dots, n-1$

(14)
\begin{eqnarray} p'_i(x) &=& -\frac{(x - x_{i+1})^2}{2 \Delta_i} M_i + \frac{(x - x_i)^2}{2 \Delta_i} M_{i+1} + \frac{y_{i+1} - y_i}{\Delta_i} - \frac{\Delta_i}{6} (M_{i+1} - M_i) \\ p'_{i-1}(x) &=& \frac{(x - x_i)^2}{2 \Delta_{i-1}} M_{i-1} + \frac{(x - x_{i-1})^2}{2 \Delta_{i-1}} M_i + \frac{y_i - y_{i-1}}{\Delta_{i-1}} - \frac{\Delta_{i-1}}{6} (M_i - M_{i - 1}) \end{eqnarray}

Тогава

(15)
\begin{eqnarray} p'_i(x_i) &=& -\frac{\Delta_i}{2}M_i + \frac{y_{i+1} - y_i}{\Delta_i} - \frac{\Delta_i}{6}(M_{i+1} - M_i) \\ p'_{i-1}(x_i) &=& \frac{\Delta_{i-1}}{2} M_{i} + \frac{y_i - y_{i-1}}{\Delta_{i-1}} - \frac{\Delta_{i-1}}{6} (M_i - M_{i-1}) \end{eqnarray}

Приравняваме и получаваме:

(16)
\begin{eqnarray} & &-\frac{\Delta_i}{2} M_i + \frac{y_{i+1} - y_i}{\Delta_i} - \frac{\Delta_i}{6} (M_{i+1} - M_i) = \frac{\Delta_{i - 1}}{2} M_i + \frac{y_i - y_{i-1}}{\Delta_{i-1}} - \frac{\Delta_{i-1}}{6} (M_i - M_{i-1}) \\ & &-\Big( \frac{\Delta_i}{3} + \frac{\Delta_{i-1}}{3}\Big) M_i - \frac{\Delta_i}{6} M_{i+1} - \frac{\Delta_{i-1}}{6} M_{i-1} + \frac{y_{i+1} - y_i}{\Delta_i} - \frac{y_i - y_{i-1}}{\Delta_{i-1}} = 0 \end{eqnarray}

Сега се сещаме, че всъщност $M_i = s''(x_i) = q_iy_i + f_i$
Заместваме и се получава следната боза:

(17)
\begin{align} -\Big( \frac{\Delta_i}{3} + \frac{\Delta_{i-1}}{3}\Big) (q_iy_i + f_i) - \frac{\Delta_i}{6} (q_{i+1}y_{i+1} + f_{i+1}) - \frac{\Delta_{i-1}}{6} (q_{i - 1}y_{i-1} + f_{i-1}) + \frac{y_{i+1} - y_i}{\Delta_i} - \frac{y_i - y_{i-1}}{\Delta_{i-1} = 0} \end{align}

Малко да разместим и да вземем предвид началните условия на граничната задача и получаваме следната система:

(18)
\begin{array} {|cccccc} &y_{i-1}\Big( \frac{1}{\Delta_{i-1}} - \frac{\Delta_{i-1}}{6}q_{i-1}\Big) - y_i \Big( \frac{1}{\Delta_{i-1}} + \frac{1}{\Delta_i} + \frac{\Delta_{i-1} + \Delta_i}{3} q_i\Big) + y_{i+1}\Big( \frac{1}{\Delta_i} - \frac{\Delta_i}{6} q_{i+1}\Big) = \frac{\Delta_{i-1}}{6} f_{i-1} + \frac{\Delta_{i-1} + \Delta_i}{3} f_i + \frac{\Delta_i}{6} f_{i+1} \\ &i = 1, 2, \dots, n-1 \\ &y_0 = \alpha \\ &y_n = \beta \end{array}
2

Получената система е с $n-1$ уравнения и $n-1$ неизвестни.
Забелязваме, че матрицата на системата е с доминиращ главен диагонал, т.е има решение (и то единствено), което намираме по метода на прогонката.

Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-NonCommercial-ShareAlike 3.0 License