RGR.tex

\documentclass[14pt]{article}

\usepackage[table]{xcolor}
\usepackage[a4paper, left=2cm, right=2cm, bottom=3cm]{geometry}
\usepackage[T1, T2A]{fontenc}
\usepackage[english, main=ukrainian]{babel}
\usepackage{tempora}
\usepackage{amssymb,amsfonts,amsmath,amsthm}   
\usepackage{sectsty}
\usepackage{graphicx}
\usepackage{caption}
\usepackage{tikz}
\usepackage{wrapfig}
\usepackage{mathtools}
\usepackage{algpseudocode}
\usepackage{algorithm}
\usepackage{cmap}                                                        
\usepackage[utf8]{inputenc}                                    
\usepackage{mathtext}
\usepackage{hyperref}
\usepackage{url}
\usepackage{listings}
\usepackage{subcaption}

\lstset{
  basicstyle=\small\ttfamily,
  columns=fullflexible,
  frame=single, % Add a frame around the code listing
}

\hypersetup{
    colorlinks,
    citecolor=black,
    filecolor=black,
    linkcolor=black,
    urlcolor=black
} 

\graphicspath{ {./Images/} }

\setcounter{MaxMatrixCols}{30}
\begin{document}

\begin{titlepage}
	\centering 
	\fontsize{14pt}{16pt}\selectfont
	НАЦІОНАЛЬНИЙ ТЕХНІЧНИЙ УНІВЕРСИТЕТ УКРАЇНИ\\
	\text{«КИЇВСЬКИЙ ПОЛІТЕХНІЧНИЙ ІНСТИТУТ імені Ігоря СІКОРСЬКОГО»}\\
	\vspace{0.1cm}
	Фізико-технічний інститут\\
	\vspace{5cm}
	
	КУРСОВА РОБОТА\\
	з кредитного модуля «Методи обчислень»\\
	на тему:\\
	«ОБЧИСЛЮВАЛЬНЕ РОЗВ’ЯЗАННЯ ДИФЕРЕНЦІАЛЬНИХ \\
	РІВНЯНЬ У ЧАСТИННИХ ПОХІДНИХ»\\
	Варіант №\underline{   9   } 
	\vfill

	\raggedleft 
	\fontsize{12pt}{14pt}\selectfont
	Виконав:\\
	студент 3 курсу ФТІ\\
	Хоменко О.В.\\
	Перевірила:\\
	Стьопочкіна І. В.\\
	Оцінка:\\
	
	\vfill
	\centering
	Київ-2023
	
\end{titlepage}

\fontsize{12pt}{16pt}\selectfont
\tableofcontents


\pagebreak
\section{Постановка задачі.}

\fontsize{12pt}{16pt}\selectfont
Процес дифузії речовини, яка розповсюджується із точкового джерела інтенсивності ${Q=10~\text{г/с}}$, розташованого в центральній точці області. Розміри області $h \times h$ м.\\
Рівняння еволюції речовини:\\


\begin{equation}\label{eq:1}
\frac{\partial u}{\partial t} = 
C \sum_{i=1}^{2}\frac{\partial^2 u}{\partial z_{i}^{2}} + 
c_{1} \frac{\partial u}{\partial z_{1}} + Q \cdot \delta \left(z_1 - \frac{h}{2}\right) \cdot \delta \left(z_2 - \frac{h}{2} \right)
\end{equation}

\begin{equation}\label{eq:2}
\left.\frac{\partial u}{\partial z_{i}}\right |_{\text{г}} = 0;
\end{equation}

\begin{equation}\label{eq:3}
\left. u \right|_{t=0} = 0;
\end{equation}

$$\text{де } C=50; ~ c_{1}=1~\text{м/с}; ~ h=100~\text{м}; ~Q = 10~\text{г/c};$$

\centering
Зробити короткостроковий прогноз (5 хв.) поведінки речовини.

\vspace{1 cm}
\raggedright
Рівняння такого типу застосовуються при описі расповсюдження деякої речовини внаслідок дифузії. В даному випадку воно еволюційне, параболічне та лінійне(оскільки коефіцент $C$ - постійна величина). Варто зауважити що джерело в нас точкове, тому в рівнянні \eqref{eq:1} фігурує дельта функція.

\pagebreak

\section{Огляд та аналіз існуючих методів чисельного розв'язання ДРПЧ.}

Існує багато методів ров'язання ДРПЧ параболічного типу. Розглянемо деякі з них: \\

\begin{enumerate}
	\item Метод скінченних різниць
	\item Метод скінченних елементів
	\item Метод ліній
	\item Метод функцій Гріна	
\end{enumerate}

\fontsize{14pt}{16pt}\selectfont
\textbf{Метод скінченних різниць} \\

\fontsize{12pt}{16pt}\selectfont
Цей метод використовує регулярні сітки. Область протікання процесу апроксимується прямокутниками, оскільки метод застосвується переважно в областях із прямолінійними границями. В зв’язку з цим МСР частіше використовується для аналізу задач з прямолінійними границями областей визначення функцій. Традиційно розв‘язувані методом МСР практичні задачі – це задачі дифузії та тепломасообміну рідин та газів.

\vspace{1 cm}
\fontsize{14pt}{16pt}\selectfont
\textbf{Метод скінченних елементів} \\

\fontsize{12pt}{16pt}\selectfont
Цей метод використовує розбиття з урахуванням геометричних особливостей області. На елементи розбивають внутрішні області, при цьому алгоритм розбиття будується так, щоб елементи задовольняли деяким обмеженням,
наприклад сторони трикутників не дуже відрізнялись по довжині і т.д. Тому областю використання МСЕ є коло задач з областю визначення з довільною геометрією, наприклад - розрахунок на міцність деталей та вузлів будівельних конструкцій, авіаційних і космічних апаратів, тепловий розрахунок двигунів і т.д. Але МСЕ набуває досить складної реалізації у випадку еволюційної задачі. Цей метод краще пристосований до застосування у випадку стаціонарних задач.

\vspace{1 cm}
\fontsize{14pt}{16pt}\selectfont
\textbf{Метод ліній} \\

\fontsize{12pt}{16pt}\selectfont
Основна ідея цього методу полягає в тому, щоб замінити просторові(граничні значення) похідних в ДРЧП алгебраїчними наближеннями. Після того, як це буде зроблено, просторові похідні більше не вказуються явно в термінах просторових незалежних змінних. Таким чином, фактично залишається лише змінна початкового значення, як правило, час у фізичній задачі. Іншими словами, лише з однією незалежною змінною, що залишилася, ми маємо систему звичайних диференціальних рінянь, яка наближає вихідне ДРЧП. Таким чином, завдання полягає в тому, щоб сформулювати наближену систему ЗДР. Як тільки це буде зроблено, ми можемо застосувати будь-який алгоритм розв'язання системи ЗДР для початкового значення, щоб обчислити наближене чисельне розв’язання ДРЧП. Таким чином, однією з головних особливостей МЛ є використання існуючих і загалом добре відомих чисельних методів для розв'язання систем ЗДР.

\vspace{1 cm}
\fontsize{14pt}{16pt}\selectfont
\textbf{Метод функцій Гріна} \\

\fontsize{12pt}{16pt}\selectfont
Початкові та граничні умови замінюються системою точкових джерел, і задача розв’язується для кожного точкового джерела. Повний
розв’язок вихідної задачі отримується в результаті сумування розв’язків для елементарних джерел.

\vspace{1cm}

Оскільки дана задача еволюційна та обмеження задані у вигляді прямокутної форми, то одразу можна сказати що МСЕ не підходить. Метод ліній найчастіше застосовується для непрямокутних областей. Метод функцій Гріна потребує знання функції Гріна, яку знайти в деяких випадках доволі таки складно. Тому був обраний метод скінченних різниць. МСР, як вже було зазначено вище, дуже часто застосовується при розв‘язувані задач дифузії та тепломасообміну рідин та газів.

\pagebreak
\section{Дослідження умов застосування обраного методу.}

\subsection{Дискретизація рівняння}

Будемо використовувати наступні позначення.\\
$0 \leq x \leq h$, $0 \leq y \leq h$ - діапазони зміни просторових координати \\ 
$0 \leq t \leq \infty$ - діапазон зміни часу \\

\vspace{0.5 cm}

Розіб’ємо відрізок для x та y - $[0,h]$ на $n$ вузлів, перший вузол співпадає з нулем, $n\text{-й}$ співпадає із $h$. Позначимо літерою $q$ індекс, який показує номер вузла на відрізку $x$, $r$ індекс, який показує номер вузла на відрізку $y$, \: $q,r=0,\dots n-1$. Позначимо відстань між сусідніми просторовими вузлами $\Delta x, \Delta y$ відповідно. Введемо індекс дискретного часу $p$, $p=0, \dots m$.
Позначимо відстань між сусідніми моментами часу $\Delta p$.



\begin{figure}[!h]
	\centering
	\includegraphics[scale=0.5]{Figure.png}
	\caption{Сітка дискретизації}
	\label{fig:method}
\end{figure}

\pagebreak

\subsection{Огляд різницевих схем}


\textbf{Можливі різницеві схеми} \\

\begin{enumerate}
	\item Пряма схема Ейлера
	\item Зворотня схема Ейлера
	\item Схема Кранка-Ніколсона
\end{enumerate}

Усі ці схеми можна отримати за допомогою так званого $\theta$-правила. Для рівняння вигляду: \\
$$\frac{\partial u}{\partial t} = G(u),$$ де $G(u)$ - деякий просторовий диференціальний оператор, $\theta$-правило виглядає так:
$$\frac{u_{i}^{n+1} - u_i^{n}}{\Delta t} = \theta G(u_i^{n+1}) + (1- \theta)G(u_i^n)$$

Тоді при різних значеннях $\theta$ будуть різні схеми: \\

\begin{itemize}
	\item $\theta = 0$ дає пряму схему Ейлера в часі
	\item $\theta = 1$ дає зворотню схему Ейлера в часі
	\item $\theta = \frac{1}{2}$ дає схему Кранка-Ніколсона в часі
\end{itemize}

\vspace{1cm}
В нашому випадку це правило буде виглядати наступним чином. \\
Запишемо скінченно-різницеві апроксимації для частинних похідних $u_{t},\: u_{x},\: u_{xx},\: u_{yy}$:

\begin{align*}
\frac{\partial u}{\partial t} &= \frac{u_{q,r}^{p+1} - u_{q,r}^{p}}{\Delta t} \\[5pt]
\frac{\partial u}{\partial x} &= \theta \cdot \frac{u_{q+1,r}^{p+1} - u_{q,r}^{p+1}}{\Delta x} + \left(1 - \theta \right) \cdot \frac{u_{q+1,r}^{p} - u_{q,r}^{p}}{\Delta x} \\[5pt]
\frac{\partial^2 u}{\partial x^2} &= \theta \cdot \frac{u_{q-1, r}^{p+1}-2u_{q, r}^{p+1} + u_{q+1,r}^{p+1}}{\Delta x^2} + \left(1 - \theta \right) \cdot \frac{u_{q-1, r}^{p}-2u_{q, r}^{p} + u_{q+1,r}^{p}}{\Delta x^2} \\[5pt]
\frac{\partial^2 u}{\partial y^2} &= \theta \cdot \frac{u_{q, r-1}^{p+1}-2u_{q, r}^{p+1} + u_{q,r+1}^{p+1}}{\Delta y^2} + \left( 1 - \theta \right) \cdot \frac{u_{q, r-1}^{p}-2u_{q, r}^{p} + u_{q,r+1}^{p}}{\Delta y^2}
\end{align*}

\pagebreak
\subsection{Оцінка схем}

Рівняння (\ref{eq:1}) в усіх точках, при умові що нема джерела та речовина ще не досягла стінок, буде мати вигляд: \\

\begin{equation}\label{eq:4}
\frac{\partial u}{\partial t} = C \sum_{i=1}^{2}\frac{\partial^2 u}{\partial z_{i}^{2}} + c_{1}
\frac{\partial u}{\partial z_{1}};
\end{equation}

Рівняння (\ref{eq:4}) має розв'язок у вигляді: \\
$$u(x,y,t) = Ae^{-C k^2 t}e^{i\left(k_x \cdot (x+c_1 t) + k_y y\right)}$$
та відповіний дискретний розв'язок у вигляді: \\
$$u_{q,r}^{p} = A \xi^{p} e^{i\left(k_x (q \Delta x + c_1 p \Delta t) + k_y r \Delta y \right)}$$


\textbf{Пряма схема Ейлера}: \\

\begin{align*}
\frac{u_{q,r}^{p+1} - u_{q,r}^{p}}{\Delta t} &= 
C \cdot \left(\frac{u_{q-1, r}^{p}-2u_{q, r}^{p} + u_{q+1,r}^{p}}{\Delta x^2} + \frac{u_{q, r-1}^{p}-2u_{q, r}^{p} + u_{q,r+1}^{p}}{\Delta y^2} \right) \\
&\qquad + c_1 \cdot \left( \frac{u_{q+1,r}^{p} - u_{q,r}^{p}}{\Delta x} \right)
\end{align*}

\vspace{1cm}
\begin{align*}
& \frac{A \xi^p e^{i\left(k_x (q \Delta x + c_1 p \Delta t) + k_y r \Delta y   \right)} \cdot \left(\xi e^{i k_x c_1 \Delta t} - 1\right)}{\Delta t} = \\[1pt] 
& C \cdot A \xi^p e^{i\left(k_x (q \Delta x + c_1 p \Delta t)+ k_y r \Delta y   \right)} \cdot \left(\frac{e^{-i k_x \Delta x} - 2 + e^{i k_x  \Delta x}}{\Delta x^2} + \frac{e^{-i k_y  \Delta y} - 2 + e^{i k_y  \Delta y}}{\Delta y^2} \right) + \\[1pt]
& \frac{c_1 \cdot A \xi^p e^{i\left(k_x (q \Delta x + c_1 p \Delta t)+ k_y r \Delta y \right)} \left( e^{i k_x \Delta x} - 1 \right)}{\Delta x}
\end{align*}

Після спрощень отримуємо: \\
$$
\frac{\xi e^{i k_x c_1 \Delta t} - 1}{\Delta t} = 
C \cdot \left(\frac{e^{-i k_x \Delta x} - 2 + e^{i k_x \Delta x}}{\Delta x^2} + \frac{e^{-i k_y \Delta y} - 2 + e^{i k_y \Delta y}}{\Delta y^2} \right) +
\frac{c_1 \cdot \left( e^{i k_x \Delta x} - 1 \right)}{\Delta x}
$$


Або ж: \\
\begin{align*}
&\frac{\xi e^{i k_x c_1 \Delta t} - 1}{\Delta t} = 
-C \frac{4}{\Delta x^2} \sin^{2}\left(\frac{k_x \Delta x}{2}\right)
-C \frac{4}{\Delta y^2} \sin^{2}\left(\frac{k_y \Delta y}{2}\right)+ \\[2pt]
&+ \frac{c_1 \cdot \left(\cos \left(k_x \Delta x \right) + i \sin \left(k_x \Delta x \right)- 1 \right)}
{\Delta x}
\end{align*}

Зробимо заміну: \\

$$p_x = \frac{k_x \Delta x}{2}, \: p_y = \frac{k_y \Delta y}{2}$$

Тоді:\\
$$
\frac{\xi e^{i k_x c_1 \Delta t} - 1}{\Delta t} = 
-C \frac{4}{\Delta x^2} \sin^{2}\left(p_x \right)
-C \frac{4}{\Delta y^2} \sin^{2}\left(p_y \right)+
\frac{c_1 \cdot \left(\cos \left(2 p_x \right) + i \sin \left(2 p_x \right)- 1 \right)}
{\Delta x}
$$

$$\xi e^{i k_x c_1 \Delta t} = 1 -C \frac{4 \Delta t}{\Delta x^2} \sin^{2}\left(p_x \right)
-C \frac{4 \Delta t}{\Delta y^2} \sin^{2}\left(p_y \right)+
\frac{c_1 \Delta t \cdot \left(\cos \left(2 p_x \right) + i \sin \left(2 p_x \right)- 1 \right)}
{\Delta x}$$

Модуль $\xi$ буде досягати максимума в точках $p_x = \frac{\pi}{2}, p_y = \frac{\pi}{2}\: $ \\
Тобто: \\
$$\xi_{max} = 1 - C \frac{4 \Delta t}{\Delta x^2} 
-C \frac{4 \Delta t}{\Delta y^2}-2
\frac{c_1 \Delta t}
{\Delta x}$$

Для стійкості розв'язку потрібно щоб $\left| \xi_{max} \right| \le 1$, тобто: \\


$$-1 \le -C \frac{2 \Delta t}{{\Delta x}^2} 
-C \frac{2 \Delta t}{{\Delta y}^2}-
\frac{c_1 \Delta t}{\Delta x} \le 0$$

Для спеціального, але расповсюдженого випадку $\Delta x = \Delta y = h$ \\

$$-1 \le -C \frac{4 \Delta t}{h^2} - \frac{c_1 \Delta t}{h} \le 0$$

Звідси: \\

$$\Delta t \le \frac{h^2}{c_1 h + 4C}$$

\vspace{2cm}
\textbf{Зворотня схема Ейлера:} \\

\begin{align*}
\frac{u_{q,r}^{p+1} - u_{q,r}^{p}}{\Delta t} &= 
C \cdot \left(\frac{u_{q-1, r}^{p+1}-2u_{q, r}^{p+1} + u_{q+1,r}^{p+1}}{\Delta x^2} + \frac{u_{q, r-1}^{p+1}-2u_{q, r}^{p+1} + u_{q,r+1}^{p+1}}{\Delta y^2} \right) \\
&\qquad + c_1 \cdot \left( \frac{u_{q+1,r}^{p+1} - u_{q,r}^{p+1}}{\Delta x} \right)
\end{align*}



\vspace{1cm}
\begin{align*}
& \frac{A \xi^p e^{i\left(k_x (q \Delta x + c_1 p \Delta t) + k_y r \Delta y \right)} \cdot \left(\xi e^{i k_x c_1 \Delta t} - 1\right)}{\Delta t} = \\[1pt] 
& C \cdot A \xi^{p+1} e^{i\left(k_x (q \Delta x + c_1 (p+1) \Delta t) + k_y r \Delta y \right)} \cdot \left(\frac{e^{-i k_x \Delta x} - 2 + e^{i k_x  \Delta x}}{\Delta x^2} + \frac{e^{-i k_y  \Delta y} - 2 + e^{i k_y  \Delta y}}{\Delta y^2} \right) + \\[1pt]
& \frac{c_1 \cdot A \xi^{p+1} e^{i\left(k_x (q \Delta x + c_1 (p + 1) \Delta t) + k_y r \Delta y \right)} \left( e^{i k_x \Delta x} - 1 \right)}{\Delta x}
\end{align*}



Після спрощень отримуємо: \\
\begin{align*}
& \frac{\xi e^{i k_x c_1 \Delta t} - 1}{\Delta t} = \\[1pt]
& C \cdot \xi e^{i k_x c_1 \Delta t} \left(\frac{e^{-i k_x \Delta x} - 2 + e^{i k_x \Delta x}}{\Delta x^2} + \frac{e^{-i k_y \Delta y} - 2 + e^{i k_y \Delta y}}{\Delta y^2} \right) +
\frac{c_1 \xi e^{i k_x c_1 \Delta t} \cdot \left( e^{i k_x \Delta x} - 1 \right)}{\Delta x}
\end{align*}

Зробимо заміну: $\hat{\xi} = \xi e^{i k_x c_1 \Delta t}$\\[3pt]
Тоді: \\
\begin{align*}
&\frac{1 - {\hat{\xi}}^{-1}}{\Delta t} = 
-C \frac{4}{\Delta x^2} \sin^{2}\left(\frac{k_x \Delta x}{2}\right)
-C \frac{4}{\Delta y^2} \sin^{2}\left(\frac{k_y \Delta y}{2}\right)+ \\[2pt]
&+ \frac{c_1 \cdot \left(\cos \left(k_x \Delta x \right) + i \sin \left(k_x \Delta x \right)- 1 \right)}
{\Delta x}
\end{align*}

Зробимо заміну:
$$p_x = \frac{k_x \Delta x}{2}, \: p_y = \frac{k_y \Delta y}{2}$$


Тоді:\\
$$
\frac{1 - {\hat{\xi}}^{-1}}{\Delta t} = 
-C \frac{4}{\Delta x^2} \sin^{2}\left(p_x \right)
-C \frac{4}{\Delta y^2} \sin^{2}\left(p_y \right)+
\frac{c_1 \cdot \left(\cos \left(2 p_x \right) + i \sin \left(2 p_x \right)- 1 \right)}
{\Delta x}
$$

$$1 - {\hat{\xi}}^{-1} = -C \frac{4 \Delta t}{\Delta x^2} \sin^{2}\left(p_x \right)
-C \frac{4 \Delta t}{\Delta y^2} \sin^{2}\left(p_y \right)+
\frac{c_1 \Delta t \cdot \left(\cos \left(2 p_x \right) + i \sin \left(2 p_x \right)- 1 \right)}
{\Delta x}$$

$$\hat{\xi} = \left(1 + C \frac{4 \Delta t}{\Delta x^2} \sin^{2}\left(p_x \right)
+C \frac{4 \Delta t}{\Delta y^2} \sin^{2}\left(p_y \right)-
\frac{c_1 \Delta t \cdot \left(\cos \left(2 p_x \right) + i \sin \left(2 p_x \right)- 1 \right)}
{\Delta x} \right)^{-1}$$


$$\left|\xi \right| = \frac{1}{\sqrt{\left(1 + C \frac{4 \Delta t}{\Delta x^2} \sin^{2}\left(p_x \right)
+C \frac{4 \Delta t}{\Delta y^2} \sin^{2}\left(p_y \right) 
- c_1 \frac{\Delta t}{\Delta x} \cdot \left(\cos\left(2p_x\right)-1 \right) \right)^2 + \left( \frac{c_1 \Delta t \sin \left(2 p_x \right)}{\Delta x}\right)^2}}$$


Максимум модуля $\xi$ буде досягатись коли вираз в знаменнику буде мінімальним. А це буде відбуватись, коли $p_x = 0, \: p_y = 0$ та значення буде дорівнювати 1. Тому $\left| \xi \right| \le 1$ завжди.

\pagebreak
\textbf{Схема Кранка-Ніколсона}: \\

\begin{align*}
\frac{u_{q,r}^{p+1} - u_{q,r}^{p}}{\Delta t} &= 
\frac{C}{2} \cdot \left(\frac{u_{q-1, r}^{p+1}-2u_{q, r}^{p+1} + u_{q+1,r}^{p+1}}{\Delta x^2} + \frac{u_{q-1, r}^{p}-2u_{q, r}^{p} + u_{q+1,r}^{p}}{\Delta x^2} \right. \\
&\qquad + \frac{u_{q, r-1}^{p+1}-2u_{q, r}^{p+1} + u_{q,r+1}^{p+1}}{\Delta y^2} + \frac{u_{q, r-1}^{p}-2u_{q, r}^{p} + u_{q,r+1}^{p}}{\Delta y^2}\bigg) \\
&\qquad + \frac{c_1}{2} \cdot \left(\frac{u_{q+1,r}^{p+1} - u_{q,r}^{p+1}}{\Delta x} + \frac{u_{q+1,r}^{p} - u_{q,r}^{p}}{\Delta x} \right)
\end{align*}

\vspace{1cm}
\begin{align*}
& \frac{A \xi^p e^{i\left(k_x (q \Delta x + c_1 p \Delta t) + k_y r \Delta y \right)} \cdot \left(\xi e^{i k_x c_1 \Delta t} - 1\right)}{\Delta t} = \\[1.5pt] 
& \frac{C}{2} \cdot A \xi^{p} e^{i\left(k_x (q \Delta x + c_1 p \Delta t) + k_y r \Delta y \right)} \left( \xi e^{i k_x c_1 \Delta t} + 1\right) \left( \frac{e^{-i k_x \Delta x} - 2 + e^{i k_x  \Delta x}}{\Delta x^2} + \frac{e^{-i k_y \Delta y} - 2 + e^{i k_y  \Delta y}}{\Delta y^2} \right) + \\[1.5pt]
& + \frac{c_1 \cdot A \xi^{p} e^{i\left(k_x (q \Delta x + c_1 p \Delta t) + k_y r \Delta y \right)} \cdot \left( \xi e^{i k_x c_1 \Delta t} + 1\right) \left( e^{i k_x \Delta x} - 1 \right)}{2 \Delta x}
\end{align*}

Зробимо заміну: $\hat{\xi} = \xi e^{i k_x c_1 \Delta t}$\\[3pt]

Після спрощень отримуємо: \\

\begin{align*}
\frac{\left(\hat{\xi} - 1\right)}{\Delta t} = 
\frac{C}{2} \cdot \left( \hat{\xi} + 1\right)\cdot \left( \frac{e^{-i k_x \Delta x} - 2 + e^{i k_x  \Delta x}}{\Delta x^2} + \frac{e^{-i k_y \Delta y} - 2 + e^{i k_y  \Delta y}}{\Delta y^2} \right) + \\[1pt]
\frac{ c_1 \left( \hat{\xi} + 1 \right) \cdot \left( e^{i k_x \Delta x} - 1 \right)  }{2 \Delta x}
\end{align*}

Або ж: \\

\begin{align*}
& \frac{\left(\hat{\xi} - 1\right)}{\Delta t} = \\[1.5pt] 
& \left( \hat{\xi} + 1\right)\cdot \left(-2C \left( \frac{ \sin^{2}\left(\frac{k_x \Delta x}{2}\right)}{\Delta x^2} + \frac{\sin^{2}\left(\frac{k_y \Delta y}{2}\right)}{\Delta y^2} \right) + \frac{c_1 \cdot \left(\cos \left(k_x \Delta x \right) + i \sin \left(k_x \Delta x \right)- 1 \right)}{2 \Delta x} \right)
\end{align*}


Зробимо заміну: \\
$$p_x = \frac{k_x \Delta x}{2}, \: p_y = \frac{k_y \Delta y}{2}$$

Тоді: \\
\begin{align*}
\frac{\left(\hat{\xi} - 1\right)}{\Delta t} = 
\left( \hat{\xi} + 1\right)\cdot \left(-2C \left( \frac{ \sin^{2}\left(p_x \right)}{\Delta x^2} + \frac{\sin^{2}\left( p_y \right)}{\Delta y^2} \right) + \frac{c_1 \cdot \left(\cos \left(2 p_x \right) + i \sin \left(2 p_x \right)- 1 \right)}{2 \Delta x} \right)
\end{align*}

\begin{align*}
\frac{\hat{\xi} - 1}{\hat{\xi} + 1} = 
\Delta t \cdot \left(-2C \left( \frac{ \sin^{2}\left(p_x \right)}{\Delta x^2} + \frac{\sin^{2}\left( p_y \right)}{\Delta y^2} \right) + \frac{c_1 \cdot \left(\cos \left(2 p_x \right) + i \sin \left(2 p_x \right)- 1 \right)}{2 \Delta x} \right)
\end{align*}

Тоді: \\
\begin{align*}
\hat{\xi} = \frac{1 -2C \Delta t \cdot \left( \frac{ \sin^{2}\left(p_x \right)}{\Delta x^2} + \frac{\sin^{2}\left( p_y \right)}{\Delta y^2} \right) + \frac{\Delta t \cdot c_1 \cdot \left(\cos \left(2 p_x \right) + i \sin \left(2 p_x \right)- 1 \right)}{2 \Delta x}}{1 +2C \Delta t \cdot \left( \frac{ \sin^{2}\left(p_x \right)}{\Delta x^2} + \frac{\sin^{2}\left( p_y \right)}{\Delta y^2} \right) - \frac{\Delta t \cdot c_1 \cdot \left(\cos \left(2 p_x \right) + i \sin \left(2 p_x \right)- 1 \right)}{2 \Delta x}}
\end{align*}


\begin{align*}
|\xi| = \sqrt{\frac{\left(1 - 2C \Delta t \cdot \left( \frac{ \sin^{2}\left(p_x \right)}{\Delta x^2} + \frac{\sin^{2}\left( p_y \right)}{\Delta y^2} \right) + \Delta t c_1 \cdot \frac{ \cos \left(2 p_x \right) - 1}{2 \Delta x} \right)^2 + \left(\frac{\Delta t \cdot c_1 \cdot \sin \left(2 p_x \right)}{2 \Delta x} \right)^2}{\left(1 + 2C \Delta t \cdot \left( \frac{ \sin^{2}\left(p_x \right)}{\Delta x^2} + \frac{\sin^{2}\left( p_y \right)}{\Delta y^2} \right) - \Delta t c_1 \cdot \frac{ \cos \left(2 p_x \right) - 1}{2 \Delta x} \right)^2 + \left(\frac{\Delta t \cdot c_1 \cdot \sin \left(2 p_x \right)}{2 \Delta x} \right)^2}}
\end{align*}

Максимум досягається у точках $p_x=0, p_y=0$ та максимальне значення модуля дорівнює 1. Отже, ця схема є стійкою завжди.

\pagebreak
\subsection{Вибір різницевої схеми}
Надалі будемо використовувати схему Кранка-Ніколсона, оскільки вона стійка та її точність більша у порівнянні з зворотною схемою Ейлера. Ця схема також є неявною.

Отже, в нас буде таке рівняння : \\

\begin{equation}
\begin{aligned}\label{eq:5}
\frac{u_{q,r}^{p+1} - u_{q,r}^{p}}{\Delta t} &= 
\frac{C}{2} \cdot \left(\frac{u_{q-1, r}^{p+1}-2u_{q, r}^{p+1} + u_{q+1,r}^{p+1}}{\Delta x^2} + \frac{u_{q-1, r}^{p}-2u_{q, r}^{p} + u_{q+1,r}^{p}}{\Delta x^2} \right. \\
&\qquad + \frac{u_{q, r-1}^{p+1}-2u_{q, r}^{p+1} + u_{q,r+1}^{p+1}}{\Delta y^2} + \frac{u_{q, r-1}^{p}-2u_{q, r}^{p} + u_{q,r+1}^{p}}{\Delta y^2}\bigg) \\
&\qquad + \frac{c_1}{2 \Delta x} \cdot \left(u_{q+1,r}^{p+1} - u_{q,r}^{p+1} + u_{q+1,r}^{p} - u_{q,r}^{p} \right)
\end{aligned}
\end{equation}


із граничними \\


\begin{align*}
& u_{q,r+1}^{p+1} - u_{q,r}^{p+1} = u_{q,r+1}^{p} - u_{q,r}^{p} && \text{при } r=0 ; ~ q=1, \ldots, n-2 \\
& u_{q,r}^{p+1} - u_{q,r-1}^{p+1} = u_{q,r-1}^{p} - u_{q,r}^{p} && \text{при } r=n-1 ; ~ q=1, \ldots, n-2 \\
& u_{q+1,r}^{p+1} - u_{q,r}^{p+1} = u_{q,r}^{p} - u_{q+1,r}^{p} && \text{при } q=0 ; ~ r=1, \ldots, n-2 \\
& u_{q,r}^{p+1} - u_{q-1,r}^{p+1} = u_{q-1,r}^{p} - u_{q,r}^{p} && \text{при } q=n-1 ; ~ r=1, \ldots, n-2
\end{align*}

\vspace{1cm}

та початковими умовами \\
$$u_{q,r}^{1} = 0, \:\: q,r=0,\dots n-1 $$

\vspace{1cm}

Перенесемо всі невідомі значення u з верхнього часового шару (з індексом t+1) в ліву частину рівняння (\ref{eq:5}) і одержемо: \\


\begin{equation*}
\begin{aligned}
&u_{q,r}^{p+1} - \frac{C \Delta t}{2} \cdot \left( \frac{u_{q-1, r}^{p+1}-2u_{q, r}^{p+1} + u_{q+1,r}^{p+1}}{\Delta x^2} + \frac{u_{q, r-1}^{p+1}-2u_{q, r}^{p+1} + u_{q,r+1}^{p+1}}{\Delta y^2} \right) - \\
&- \frac{c_1 \Delta t}{2 \Delta x} \cdot \left( u_{q+1,r}^{p+1} - u_{q,r}^{p+1} \right) = \frac{C \Delta t}{2} \cdot \left(\frac{u_{q-1, r}^{p}-2u_{q, r}^{p} + u_{q+1,r}^{p}}{\Delta x^2}
+ \frac{u_{q, r-1}^{p}-2u_{q, r}^{p} + u_{q,r+1}^{p}}{\Delta y^2} \right)+\\
&+ \frac{c_1 \Delta t}{2 \Delta x} \cdot \left(u_{q+1,r}^{p} - u_{q,r}^{p} \right) + u_{q,r}^{p}
\end{aligned}
\end{equation*}

\pagebreak
Згрупуємо невідомі значення та отримуємо кінцеве рівняння: \\ 

\begin{equation}
\begin{aligned}\label{eq:6}
&u_{q,r}^{p+1} \left(1 + F_x + F_y + \frac{c_1 \Delta t}{2 \Delta x} \right) - u_{q+1, r}^{p+1} \cdot \frac{\Delta t}{2 \Delta x} \left( \frac{C}{\Delta x} + c_1  \right) - \\[3pt] 
& - u_{q-1,r}^{p+1} \cdot \frac{F_x}{2} - u_{q,r-1}^{p+1} \cdot \frac{F_y}{2} - u_{q,r+1}^{p+1} \cdot \frac{F_y}{2} = \\[3pt]
& \frac{F_x}{2} \cdot \left(u_{q-1, r}^{p}-2u_{q, r}^{p} + u_{q+1,r}^{p}\right)
+ \frac{F_y}{2} \cdot \left(u_{q, r-1}^{p}-2u_{q, r}^{p} + u_{q,r+1}^{p}\right)+ \\[3pt]
&+ \frac{c_1 \Delta t}{2 \Delta x} \cdot \left(u_{q+1,r}^{p} - u_{q,r}^{p} \right) + u_{q,r}^{p}
\end{aligned}
\end{equation}


\vspace{0.5cm}
$F_x = \frac{C \Delta t}{{\Delta x}^2},
F_y = \frac{C \Delta t}{{\Delta y}^2} \text{ - числа Фур'є сітки}$
\vspace{1cm}

Якщо $p$ - фіксоване, а $q$ та $r$ змінюється від 1 до $n-2$, співвідношення (\ref{eq:6}) разом з граничними умовами повністю визначає усі $n^2 - 4$ невідомі. Точки $u_{0,0}, u_{0,n-1}, u_{n-1,0}, u_{n-1,n-1} $ неможливо визначити, оскільки граничні умови у вигляді Фон-Неймана не дозволяють це зробити (немає похідних у кутах). Тому для цього виродженого випадку покладемо що вони будуть рівні півсумі сусідніх вузлів.

\vspace{1cm}
В нашому випадку, загальна кількість невідомих - $n^2$. Щоб їх знайти, потрібно сформувати рівняння у матричному вигляді($Ax=b$). Вектор розв'язку $x$ може мати тільки один індекс, тому потрібно зробити деяке відображення двовимірних індексів в одновимірний індекс. Задамо наступне відображення $m$:

$$m\left(q, r\right) \mapsto r \cdot n + q$$
\vspace{1 cm}


\pagebreak


Розглянемо наприклад n = 5. Тоді сітка буде розміру $5 \times 5$
\vspace{0.5cm}

\begin{tikzpicture}[scale=2]
    % Draw horizontal lines
    \foreach \x in {0,1,...,4} {
        \draw (\x cm,0cm) -- (\x cm,4cm);
    }
    % Draw vertical lines
    \foreach \y in {0,1,...,4} {
        \draw (0cm,\y cm) -- (4cm,\y cm);
    }
    % Draw coordinates
    \foreach \r in {0,1,...,4} {
        \foreach \q in {0,1,...,4} {
        		\pgfmathsetmacro{\coord}{int(\q+\r*5)}
            \filldraw (\q cm, \r cm) circle (1pt) node[anchor=south west, xshift=-3pt, yshift=-3pt, font=\small] {(\q,\r):{\coord}};
        }
    }
\end{tikzpicture}

\vspace{0.5cm}
Запишемо структуру матриці коефіцентів для системи рівнянь.\\
Позначимо $\cdot$ - нульовий коефіцент, $\bullet$ - ненульовий. Варто зауважити, що вона отримується як матриця коефіцентів рівнянь для кожного послідовного вузла): \\

\[
\begin{psmallmatrix}
\bullet & \bullet & \cdot & \cdot & \cdot & \bullet & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\
\cdot & \bullet & \cdot & \cdot & \cdot & \cdot & \bullet & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\
\cdot & \cdot & \bullet & \cdot & \cdot & \cdot & \cdot & \bullet & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\
\cdot & \cdot & \cdot & \bullet & \cdot & \cdot & \cdot & \cdot & \bullet & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\
\cdot & \cdot & \cdot & \bullet & \bullet & \cdot & \cdot & \cdot & \cdot & \bullet & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\
\cdot & \cdot & \cdot & \cdot & \cdot & \bullet & \bullet & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\
\cdot & \bullet & \cdot & \cdot & \cdot & \bullet & \bullet & \bullet & \cdot & \cdot & \cdot & \bullet & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\
\cdot & \cdot & \bullet & \cdot & \cdot & \cdot & \bullet & \bullet & \bullet & \cdot & \cdot & \cdot & \bullet & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\
\cdot & \cdot & \cdot & \bullet & \cdot & \cdot & \cdot & \bullet & \bullet & \bullet & \cdot & \cdot & \cdot & \bullet & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\
\cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \bullet & \bullet & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\
\cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \bullet & \bullet & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\
\cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \bullet & \cdot & \cdot & \cdot & \bullet & \bullet & \bullet & \cdot & \cdot & \cdot & \bullet & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\
\cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \bullet & \cdot & \cdot & \cdot & \bullet & \bullet & \bullet & \cdot & \cdot & \cdot & \bullet & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\
\cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \bullet & \cdot & \cdot & \cdot & \bullet & \bullet & \bullet & \cdot & \cdot & \cdot & \bullet & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\
\cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \bullet & \bullet & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\
\cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \bullet & \bullet & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\
\cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \bullet & \cdot & \cdot & \cdot & \bullet & \bullet & \bullet & \cdot & \cdot & \cdot & \bullet & \cdot & \cdot & \cdot \\
\cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \bullet & \cdot & \cdot & \cdot & \bullet & \bullet & \bullet & \cdot & \cdot & \cdot & \bullet & \cdot & \cdot \\
\cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \bullet & \cdot & \cdot & \cdot & \bullet & \bullet & \bullet & \cdot & \cdot & \cdot & \bullet & \cdot \\
\cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \bullet & \bullet & \cdot & \cdot & \cdot & \cdot & \cdot \\
\cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \bullet & \cdot & \cdot & \cdot & \cdot & \bullet & \bullet & \cdot & \cdot & \cdot \\
\cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \bullet & \cdot & \cdot & \cdot & \cdot & \bullet & \cdot & \cdot & \cdot \\
\cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \bullet & \cdot & \cdot & \cdot & \cdot & \bullet & \cdot & \cdot \\
\cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \bullet & \cdot & \cdot & \cdot & \cdot & \bullet & \cdot \\
\cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \bullet & \cdot & \cdot & \cdot & \bullet & \bullet \\
\end{psmallmatrix}
\]

\vspace{1cm}
Одержана матриця є сильно розрідженою, і має стрічкову структуру з шириною 5. Зберігати її в пам'яті повністю не має сенсу, оскільки це буде займати $O(n^2)$ пам'яті. З того, що для вдалої візуалізації та моделювання потрібно брати $n > 70$, робимо висновок що потрібен інший підхід для зберігання таких матриць.

\pagebreak

\subsection{Ефективне зберігання розріджених матриць стрічкового вигляду}

Розглянемо два методи зберігання таких матриць: CSR та DIA.\\[1cm]

\textbf{DIA}\\
Цей метод базується на тому, що зберігаються лише певні діагоналі елементи розрідженної матриці у вигляді двовимірного масиву та масив зсувів.\\[3pt]

Приклад: \\

\[
A = 
\begin{pmatrix}
\cellcolor{blue!25} 1.0 & 0 & 0 & \cellcolor{purple!25} 1.6 & 0 \\
0 & \cellcolor{blue!25} 2.0 & 0 & 0 & \cellcolor{purple!25} 1.7 \\
0 & 0 & \cellcolor{blue!25} 3.0 & 0 & 0 \\
0 & 0 & 0 & \cellcolor{blue!25} 4.0 & 0 \\
\cellcolor{orange!25} 1.5 & 0 & 0 & 0 & \cellcolor{blue!25} 2.0 \\
\end{pmatrix}
\]

Дані: \[
\begin{array}{|c|c|c|c|c|}
\hline
1.5 & 0 & 0 & 0 & 0 \\
\hline
1.0 & 2.0 & 3.0 & 4.0 & 5.0 \\
\hline
0 & 0 & 0 & 1.6 & 1.7 \\
\hline
\end{array}
\]

Зсуви: \[
\begin{array}{|c|c|c|}
\hline
-4 & 0 & 3\\
\hline
\end{array}
\]

Таким чином вже потрібно зберігати лише ненульові діагоналі. Тому використовується пам'ять як $O(n \cdot m)$, де $m$ - кількість ненульових діагоналей.

\vspace{1cm}

\textbf{CSR}\\
Цей метод можна зустріти в декількох варіаціях. Ми будемо розглядати такий, який реалізовано в пакеті scipy.\\
Compressed Sparse Row - формат який представляє розріджену матрицю використовуючи три масиви: data, indices, indptr. Масив data зберігає усі ненульові елементи матриці в послідовному порядку(зліва направо, зверху вниз). Масив indices зберігає індекси стовпців. Масив indptr(index pointer) вказує на індекс в масиві індексів на перший ненульовий елемент у рядку. В кінці indptr вказується кількість ненульових елементів матриці.

\vspace{1cm}

Приклад: \\
\[
A = 
\begin{pmatrix}
1.0 & 0 & 0 & 1.6 & 0 \\
0 &  2.0 & 0 & 0 &  1.7 \\
0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 4.0 & 0 \\
1.5 & 0 & 0 & 0 & 2.0 \\
\end{pmatrix}
\]

\[
\begin{array}{|c|c|c|c|c|c|c|c|}
\hline
\text{data} & 1.0 & 1.6 & 2.0 & 1.7 & 4.0 & 1.5 & 2.0 \\
\hline
\text{indices} & 0 & 3 & 1 & 4 & 3 & 0 & 4 \\
\hline
\text{indptr} & 0 & 2 & 4 & 4 & 5 & 7 &  \\
\hline
\end{array}
\]

За допомогою такого представлення можна досить ефективно шукати ненульові елементи матриці. Пам'ять використовується як $O(m)$, де $m$ - кількість ненульових елементів матриці.

\vspace{1cm}

\textbf{CSC}\\
Цей метод зберігання розрідженної матриці є повним аналогом методу CSR, тільки замість рядків ми використовуємо стовпчики.

\pagebreak
\subsection{Огляд методів розв'язання СЛАР}

В даній роботі будуть розглянуті декілька як прямих, так й ітераційних методів розв'язання СЛАР. Будуть розглянути наступні прямі методи: LU - декомпозиції,  та такі ітеративні: метод Якобі, метод Гаусса-Зейделя, метод релаксації, метод спряжених градієнтів.


\vspace{1cm}
\textbf{Метод LU декомпозиції:}

В нашому випадку можна помітити, що матриця $A$ фактично ніяк не змінюється з часом, в той час як ветор $b$ змінюється. Головна ідея в використанні LU-декомпозиції це відокремити фазу факторизації від фактичної фази розв'язання. На етапі факторизації потрібна лише матриця $A$, в той час як для фактичного розв'язання використовується розкладена матриця $A$ та права частина для вирішення лінійної системи. Складність LU-факторизації займає $O(n^3)$. Як тільки ми знайдемо цю факторизацію(матриці L, U), складність розв'язання становить $O(n^2)$, оскільки маємо два рівняння з трикутними матрицями $(Ly=b, Ux=y)$. Варто зауважити, що не завжди LU декомпозиція існує. Тому застосовують часткове впорядкування до матриці $A$ за допомогою матриці впорядкування $P$. Такий метод називають вже PLU декомпозиція, але в усіх бібліотеках це також називають LU декомпозиція. Надалі будемо застосовувати метод розв'язання розріджених та щільних систем лінійних рівнянь який реалізовано в пакеті scipy.
На даному етапі вже можемо порівняти два підходи. А саме: розв'язання системи лінійних рівнянь з використанням LU декомпозиції для щільної матриці, та для матриці що зберігається у форматі CSС. 



\begin{figure}[!h]
  \centering
  \includegraphics[scale=0.7] {LU_sparse_dense_comparison.png}
  \caption{Порівняння розріджених та щільних матриць}
  \label{fig:lu_comparison}
\end{figure}

\pagebreak

\textbf{Ітеративні методи:}\\
Поки що ми показали як розв'язувати СЛАР за допомогою методу LU декомпозиції. Інший підхід полягає в тому, що велику матрицю коефіцентів можна й не зберігати, а застосовувати ітеративні методи для розв'язання СЛАР.

\vspace{1cm}

\textbf{Метод Якобі:\\}
Щоб пояснити цей метод, нагадаємо кінцеве рівняння яке ми отримали:

\begin{equation*}
\begin{aligned}\
&u_{q,r}^{p+1} \left(1 + F_x + F_y + \frac{c_1 \Delta t}{2 \Delta x} \right) - u_{q+1, r}^{p+1} \cdot \frac{\Delta t}{2 \Delta x} \left( \frac{C}{\Delta x} + c_1  \right) - \\[3pt] 
& - u_{q-1,r}^{p+1} \cdot \frac{F_x}{2} - u_{q,r-1}^{p+1} \cdot \frac{F_y}{2} - u_{q,r+1}^{p+1} \cdot \frac{F_y}{2} = \\[3pt]
& \frac{F_x}{2} \cdot \left(u_{q-1, r}^{p}-2u_{q, r}^{p} + u_{q+1,r}^{p}\right)
+ \frac{F_y}{2} \cdot \left(u_{q, r-1}^{p}-2u_{q, r}^{p} + u_{q,r+1}^{p}\right)+ \\[3pt]
&+ \frac{c_1 \Delta t}{2 \Delta x} \cdot \left(u_{q+1,r}^{p} - u_{q,r}^{p} \right) + u_{q,r}^{p}
\end{aligned}
\end{equation*}

Введемо номер ітерації $h$. Тоді $u_{q,r}^{p+1, h}$ буде апроксимацією $u_{q,r}^{p+1}$ на $h$ ітерації.
Суть метода Якобі полягає в тому, щоб знайти $u_{q,r}^{p+1, h}$ через вже відомі значення $u_{q \pm 1,r \pm 1}^{p+1, h-1}$ з минулої ітерації.\\
Спочатку виразимо $u_{q,r}^{p+1}$ як невідоме з нашого рівняння:\\

\begin{equation*}
\begin{aligned}\
&u_{q,r}^{p+1} = \left(1 + F_x + F_y + \frac{c_1 \Delta t}{2 \Delta x} \right)^{-1} \left(u_{q+1, r}^{p+1} \cdot \frac{\Delta t}{2 \Delta x} \left( \frac{C}{\Delta x} + c_1  \right) + \right. \\[3pt] 
& + u_{q-1,r}^{p+1} \cdot \frac{F_x}{2} + u_{q,r-1}^{p+1} \cdot \frac{F_y}{2} + u_{q,r+1}^{p+1} \cdot \frac{F_y}{2} + \\[3pt]
& + \frac{F_x}{2} \cdot \left(u_{q-1, r}^{p}-2u_{q, r}^{p} + u_{q+1,r}^{p}\right)
+ \frac{F_y}{2} \cdot \left(u_{q, r-1}^{p}-2u_{q, r}^{p} + u_{q,r+1}^{p}\right)+ \\[3pt]
&+ \left. \frac{c_1 \Delta t}{2 \Delta x} \cdot \left(u_{q+1,r}^{p} - u_{q,r}^{p} \right) + u_{q,r}^{p} \right)
\end{aligned}
\end{equation*}

В ітеративній формі це буде виглядати наступним чином:\\

\begin{equation}
\begin{aligned}\label{eq:7}
&u_{q,r}^{p+1, h+1} = \left(1 + F_x + F_y + \frac{c_1 \Delta t}{2 \Delta x} \right)^{-1} \left(u_{q+1, r}^{p+1, h} \cdot \frac{\Delta t}{2 \Delta x} \left( \frac{C}{\Delta x} + c_1  \right) + \right. \\[3pt] 
& + u_{q-1,r}^{p+1, h} \cdot \frac{F_x}{2} + u_{q,r-1}^{p+1, h} \cdot \frac{F_y}{2} + u_{q,r+1}^{p+1, h} \cdot \frac{F_y}{2} + \\[3pt]
& + \frac{F_x}{2} \cdot \left(u_{q-1, r}^{p}-2u_{q, r}^{p} + u_{q+1,r}^{p}\right)
+ \frac{F_y}{2} \cdot \left(u_{q, r-1}^{p}-2u_{q, r}^{p} + u_{q,r+1}^{p}\right)+ \\[3pt]
&+ \left. \frac{c_1 \Delta t}{2 \Delta x} \cdot \left(u_{q+1,r}^{p} - u_{q,r}^{p} \right) + u_{q,r}^{p} \right)
\end{aligned}
\end{equation}

З граничними умовами:
\begin{align*}
& u_{q,r}^{p+1, h+1} = u_{q,r+1}^{p+1, h} - u_{q,r+1}^{p} + u_{q,r}^{p} && \text{при } r=0 ; ~ q=1, \ldots, n-2 \\
& u_{q,r}^{p+1, h} = u_{q,r-1}^{p+1, h} + u_{q,r-1}^{p} - u_{q,r}^{p} && \text{при } r=n-1 ; ~ q=1, \ldots, n-2 \\
& u_{q,r}^{p+1, h+1} = u_{q+1,r}^{p+1, h} - u_{q,r}^{p} + u_{q+1,r}^{p}  && \text{при } q=0 ; ~ r=1, \ldots, n-2 \\
& u_{q,r}^{p+1, h+1} = u_{q-1,r}^{p+1, h} + u_{q-1,r}^{p} - u_{q,r}^{p} && \text{при } q=n-1 ; ~ r=1, \ldots, n-2
\end{align*}



Варто також зауважити, що $u_{q,r}^{p+1, 0} = u_{q,r}^{p} ~~ \forall q, r \in \overline{0, n-1}$.

Зазвичай в ітеративних методах використовують так звану релаксацію. Це означає що нова апроксимація це зважене середнє розв'язку запропонованого алгоритмом та минулої апроксимації. Тобто:
$$u_{q, r}^{p+1,h+1} = \omega u_{q, r}^{p+1,*} + (1 - \omega) u_{q, r}^{p + 1, h}$$

Критерій завершення ітерацій можна задавати декількома способами:
$$\max_{q, r} |u_{q,r}^{p+1,h+1} - u_{q,r}^{p+1,h}| \le \mathcal{E}$$

Або за допомогою норми $l_2$:
$$\left(\Delta x \Delta y \sum_{q, r} (u_{q, r}^{p+1,h+1} - u_{q, r}^{p+1,h})^2 \right)^\frac{1}{2} \le \mathcal{E}$$

Також користуються методом залишків для того, щоб перевіряти збіжність. Якщо знайшли розв'язок на $h+1$ ітерації, то залишок 
буде мати вигляд:
\begin{equation*}
\begin{aligned}\
&R_{q, r} = u_{q,r}^{p+1, h+1} \left(1 + F_x + F_y + \frac{c_1 \Delta t}{2 \Delta x} \right) - u_{q+1, r}^{p+1, h+1} \cdot \frac{\Delta t}{2 \Delta x} \left( \frac{C}{\Delta x} + c_1  \right) - \\[3pt] 
& - u_{q-1,r}^{p+1, h+1} \cdot \frac{F_x}{2} - u_{q,r-1}^{p+1, h+1} \cdot \frac{F_y}{2} - u_{q,r+1}^{p+1, h+1} \cdot \frac{F_y}{2} - \\[3pt]
& - \frac{F_x}{2} \cdot \left(u_{q-1, r}^{p}-2u_{q, r}^{p} + u_{q+1,r}^{p}\right)
- \frac{F_y}{2} \cdot \left(u_{q, r-1}^{p}-2u_{q, r}^{p} + u_{q,r+1}^{p}\right)- \\[3pt]
&- \frac{c_1 \Delta t}{2 \Delta x} \cdot \left(u_{q+1,r}^{p} - u_{q,r}^{p} \right) - u_{q,r}^{p}
\end{aligned}
\end{equation*}

Потім можна робити ітерації доти, поки:
$$\left(\Delta x \Delta y \sum_{q,r} R_{q, r}^2 \right)^{\frac{1}{2}} \le \mathcal{E}$$

На практиці будемо використовувати 1 критерій.\\[3pt]

Метод Якобі для нормальної збіжності потребує діагонально - панівної матриці. В нашому випадку це виконується, а тому можна реалізувувати цей метод.\\
Однак, можна показати що при $\omega > 1$ цей метод є нестійким. До того ж, вн досить повільно збігається до точного розв'язку. Тому використовують інші методи, як наприклад метод Гаусса-Зейделя, який є модифікацією метода Якобі.

\vspace{1cm}

\textbf{Метод Гаусса-Зейделя:\\}

З методу Якобі можна помітити, що при обчисленні значення $u_{q,r}^{p+1,h+1}$ ми вже обраховували значення $u_{q-1,r}^{p+1, h+1}$ та $u_{q,r-1}^{p+1, h+1}$. Тому будемо використовувати саме їх замість $u_{q-1,r}^{p+1, h}$ та $u_{q,r-1}^{p+1,h}$ відповідно. Тоді формула (\ref{eq:7}) набуває вигляду:

\begin{equation}
\begin{aligned}\label{eq:8}
&u_{q,r}^{p+1, h+1} = \left(1 + F_x + F_y + \frac{c_1 \Delta t}{2 \Delta x} \right)^{-1} \left(u_{q+1, r}^{p+1, h} \cdot \frac{\Delta t}{2 \Delta x} \left( \frac{C}{\Delta x} + c_1  \right) + \right. \\[3pt] 
& + u_{q-1,r}^{p+1, h+1} \cdot \frac{F_x}{2} + u_{q,r-1}^{p+1, h+1} \cdot \frac{F_y}{2} + u_{q,r+1}^{p+1, h} \cdot \frac{F_y}{2} + \\[3pt]
& + \frac{F_x}{2} \cdot \left(u_{q-1, r}^{p}-2u_{q, r}^{p} + u_{q+1,r}^{p}\right)
+ \frac{F_y}{2} \cdot \left(u_{q, r-1}^{p}-2u_{q, r}^{p} + u_{q,r+1}^{p}\right)+ \\[3pt]
&+ \left. \frac{c_1 \Delta t}{2 \Delta x} \cdot \left(u_{q+1,r}^{p} - u_{q,r}^{p} \right) + u_{q,r}^{p} \right)
\end{aligned}
\end{equation}

Варто зауважити, що цей метод збігається коли матриця $A$ є симетричною та додатньо визначеною, або є діагонально панівною. В даній задачі вона діагональна панівна, тому буде гарна збіжність.\\
Тут також можна використовувати релаксацію. Тобто:
$$u_{q, r}^{p+1,h+1} = \omega u_{q, r}^{p+1,*} + (1 - \omega) u_{q, r}^{p + 1, h}$$

При $0 < \omega \le 2$ метод Гаусса-Зейделя називають методом релаксації (Successive over-relaxation method).

Наразі вже можемо порівняти два методи з точки зору максимальної кількості ітерацій та часу виконання для різного розміру сітки n:


\begin{figure}[!h]
  \centering
  \includegraphics[scale=0.7] {Jacobi_Seidel_comparison.png}
  \caption{Порівняння методу Якобі та Гаусса-Зейделя}
  \label{fig:lu_comparison}
\end{figure}


\textbf{Інші методи розв'язання СЛАР:\\}
В даному пункті розглянемо два більш загальні методи, які насправді є методами підпростору Крилова.
\vspace{1cm}

\textbf{BiCGSTAB:\\}
Оскільки в даній задачі ми маємо несиметричну матрицю, то будемо розгладати модифікований алгоритм методу біспряжених градієнтів, а саме: BiCGSTAB. На відміну від методу спряжених градієнтів, cтабілізований метод біспряжених градієнтів збігається швидше та більш плавно. Наведемо його алгоритм для розв'язання системи $Ax=b$:\\

\vspace{2cm}
\begin{algorithm}[H]
	\fontsize{13pt}{18pt}\selectfont
	\caption*{\Large BiCGSTAB Algorithm}
	\begin{algorithmic}[1]
        \State $r_0 = Ax_0 - b$
        \State Обрати такий $\hat{r_0}$, що $(\hat{r_0}, r_0) \neq 0$
        \State $\rho = \alpha = \omega_0 = 1$
        \State $p_0 = r_0$
        	\State $j = 0$ \Comment{вказує на номер ітерації}
        
        \While{$||r|| > \varepsilon$}  % Виконувати поки виконується умова
            \State $\alpha_j = \frac{(r_j, \hat{r_0})}{(Ap_j, \hat{r_0})}$
            \State $s_j = r_j - \alpha_jAp_j$
            \State $\omega_j = \frac{(As_j, s_j)}{(As_j, As_j)}$
            \State $x_{j+1} = x_j + \alpha_j p_j + \omega_j s_j$
            \State $r_{j+1} = s_j - \omega_jAs_j$
            \State $\beta_j = \frac{(r_{j+1}, \hat{r_0})}{(r_j, \hat{r_0})} \frac{\alpha_j}{\omega_j}$
            \State $p_{j+1} = r_{j+1} + \beta_j(p_j - \omega_jAp_j)$
        \EndWhile
    
    \end{algorithmic}
\end{algorithm}


\pagebreak
\textbf{GMRes}

Цей метод є узагальненням методу мінімальної нев'язки, що апроксимує розв'язок вектором з мінімальної нев'язкою.
Наведемо його алгоритм для розв'язання системи $Ax=b$:\\

\vspace{2cm}
\begin{algorithm}[H]
		\fontsize{13pt}{18pt}\selectfont
        \caption*{\Large GMRes Algorithm}
        \begin{algorithmic}[1]     
            \State $r_0 = Ax_0 - b$
            \State $v_1 = \frac{r_0}{||r_0||}$
            
            \For{$j = 1$ to $m$}
                \State $h_{ij} = (Av_j, v_i) \quad i \in \{ 1 \cdots j \}$
                \State $\hat{v}_{j+1} = Av_j - \sum\limits_{i=1}^j h_{ij}v_i$
                \State $h_{j+1;j} = ||\hat{v}_{j+1}||$
                \State $v_{j+1} = \frac{\hat{v}_{j+1}}{h_{j+1;j}}$
            \EndFor
            
            \State $y_m = \arg \min ||||r_0|| e_1 - \hat{H}y||$
            \State $x_m = x_0 + V_m y_m$
            
            \If{$||b-Ax_m|| > \varepsilon$}
                \State $x_0 = x_m$
                \State $v_1 = \frac{r_m}{||r_m||}$
                \State Перейти до кроку 4
            \EndIf
            
        \end{algorithmic}
\end{algorithm}

\pagebreak
\section{Опис програмної реалізації.}

Оскільки метод Гаусса-Зейделя та Якобі досить повільно збігаються до точного розв'язку, то надалі для розв'язання нашої задачі будемо використовувати методи BiCGSTAB та GMRes. Ці алгоритми вже реалізовані в пакеті scipy, тому будемо використовувати їх звідти. Також будемо зберігати розріджену матрицю у форматі CSR для більш ефективної реалізації.\\

\textbf{Програма:}

\lstinputlisting[language=Python]{bicgstab_solver.py}

Були розглянуті декілька розмірів сіток, а саме:
$n = \{51, 71, 101, 201, 301\}$ та застосовані алгоритми BiCGSTAB та GMRes з точністю $\mathcal{E} = 10^{-5}$ на проміжках часу в 1, 10, 50, 100, 150, 200, 250, 300 секунд.

Були отримані наступні результати для сітки $101 \times 101$:\\


\begin{figure}[ht]
    \centering
        \begin{subfigure}[b]{0.5\textwidth}
            \includegraphics[scale=0.6]{1_bicgstab_101x101.png}
            \caption{1 секунда}
        \end{subfigure}%
        \begin{subfigure}[b]{0.5\textwidth}
            \includegraphics[scale=0.6]{10_bicgstab_101x101.png}
            \caption{10 секунда}
        \end{subfigure}
\end{figure}
\begin{figure}[ht]
    \ContinuedFloat
        \begin{subfigure}[b]{0.5\textwidth}
            \includegraphics[scale=0.6]{50_bicgstab_101x101}
            \caption{50 секунда}
        \end{subfigure}%
        \begin{subfigure}[b]{0.5\textwidth}
            \includegraphics[scale=0.6]{100_bicgstab_101x101}
            \caption{100 секунда}
        \end{subfigure}
        \begin{subfigure}[b]{0.5\textwidth}
            \includegraphics[scale=0.6]{150_bicgstab_101x101.png}
            \caption{150 секунда}
        \end{subfigure}%
        \begin{subfigure}[b]{0.5\textwidth}
            \includegraphics[scale=0.6]{200_bicgstab_101x101}
            \caption{200 секунда}
        \end{subfigure}
        \begin{subfigure}[b]{0.5\textwidth}
            \includegraphics[scale=0.6]{250_bicgstab_101x101}
            \caption{250 секунда}
        \end{subfigure}%
        \begin{subfigure}[b]{0.5\textwidth}
            \includegraphics[scale=0.6]{300_bicgstab_101x101}
            \caption{300 секунда}
        \end{subfigure}
        \caption{Результати розв'язку методом BiCGStab на сітці 101х101}
        \label{fig_bicgstab_101}
\end{figure}
\clearpage

\begin{figure}[ht]
    \centering
        \begin{subfigure}[b]{0.5\textwidth}
            \includegraphics[scale=0.6]{1_gmres_101x101.png}
            \caption{1 секунда}
        \end{subfigure}%
        \begin{subfigure}[b]{0.5\textwidth}
            \includegraphics[scale=0.6]{10_gmres_101x101.png}
            \caption{10 секунда}
        \end{subfigure}
\end{figure}
\begin{figure}[ht]
    \ContinuedFloat
        \begin{subfigure}[b]{0.5\textwidth}
            \includegraphics[scale=0.6]{50_gmres_101x101}
            \caption{50 секунда}
        \end{subfigure}%
        \begin{subfigure}[b]{0.5\textwidth}
            \includegraphics[scale=0.6]{100_gmres_101x101}
            \caption{100 секунда}
        \end{subfigure}
        \begin{subfigure}[b]{0.5\textwidth}
            \includegraphics[scale=0.6]{150_gmres_101x101.png}
            \caption{150 секунда}
        \end{subfigure}%
        \begin{subfigure}[b]{0.5\textwidth}
            \includegraphics[scale=0.6]{200_gmres_101x101}
            \caption{200 секунда}
        \end{subfigure}
        \begin{subfigure}[b]{0.5\textwidth}
            \includegraphics[scale=0.6]{250_gmres_101x101}
            \caption{250 секунда}
        \end{subfigure}%
        \begin{subfigure}[b]{0.5\textwidth}
            \includegraphics[scale=0.6]{300_gmres_101x101}
            \caption{300 секунда}
        \end{subfigure}
        \caption{Результати розв'язку методом GMRes на сітці 101х101}
        \label{fig_gmres_101}
\end{figure}
\clearpage

Точний розв'язок використовуючи LU декомпозицію:

\begin{figure}[ht]
    \centering
        \begin{subfigure}[b]{0.5\textwidth}
            \includegraphics[scale=0.6]{1_explicit_101x101.png}
            \caption{1 секунда}
        \end{subfigure}%
        \begin{subfigure}[b]{0.5\textwidth}
            \includegraphics[scale=0.6]{10_explicit_101x101.png}
            \caption{10 секунда}
        \end{subfigure}
\end{figure}
\begin{figure}[ht]
    \ContinuedFloat
        \begin{subfigure}[b]{0.5\textwidth}
            \includegraphics[scale=0.6]{50_explicit_101x101}
            \caption{50 секунда}
        \end{subfigure}%
        \begin{subfigure}[b]{0.5\textwidth}
            \includegraphics[scale=0.6]{100_explicit_101x101}
            \caption{100 секунда}
        \end{subfigure}
        \begin{subfigure}[b]{0.5\textwidth}
            \includegraphics[scale=0.6]{150_explicit_101x101.png}
            \caption{150 секунда}
        \end{subfigure}%
        \begin{subfigure}[b]{0.5\textwidth}
            \includegraphics[scale=0.6]{200_explicit_101x101}
            \caption{200 секунда}
        \end{subfigure}
        \begin{subfigure}[b]{0.5\textwidth}
            \includegraphics[scale=0.6]{250_explicit_101x101}
            \caption{250 секунда}
        \end{subfigure}%
        \begin{subfigure}[b]{0.5\textwidth}
            \includegraphics[scale=0.6]{300_explicit_101x101}
            \caption{300 секунда}
        \end{subfigure}
        \caption{Результати розв'язку методом LU декомпозиції на сітці 101х101}
        \label{fig_explicit_101}
\end{figure}
\clearpage

\pagebreak

\section{Висновки}
Зробивши реалізації декількох як прямих так й ітераційних методів можемо зробити висновки який метод є кращим. По перше, зрозуміло що матриці потрібно зберігати у розрідженій формі. Таким чином ми зекономимо як пам'ять так й пришвидшуємо розв'язок. Прямий метод, а саме метод LU декомпозиції із зберіганням розрідженної матриці працює швидше всіх методів. Він може здаватися найкращим за ітеративні, але це не завжди так. Іноді потрібно розв'язати задачу лише з певною точністю і коли ми будемо мати дуже великі сітки метод LU декомпозиції буде накопичувати великі похибки обчислень. Стосовно ітеративних методів можна зауважити, що метод Якобі працює повільніше за метод Гаусса-Зейделя. З іншого боку вони разом працюють в рази повільніше аніж методи BiCGSTAB та GMRes. Тому, коли мають дуже великі сітки використовують саме їх. Також при знаходженні розв'язку було помічено що GMRes працює повільніше за BiCGSTAB. А в BiCGSTAB може розходитись розв'язок при великому розмірі сітки. 

\pagebreak

\section{Список використаних джерел}

\begin{thebibliography}{}

\bibitem{fdm-book} Hans Petter Langtangen. "Numerical Methods for Partial Differential Equations: Finite Difference and Finite Volume Methods." Available at: \url{https://hplgit.github.io/fdm-book/doc/pub/diffu/pdf/diffu-4print.pdf}

\bibitem{crank1947} Crank, J.; Nicolson, P. (1947). "A practical method for numerical evaluation of solutions of partial differential equations of the heat conduction type."

\bibitem{thomas1995} Thomas, J. W. (1995). Numerical Partial Differential Equations: Finite Difference Methods. Texts in Applied Mathematics. Vol. 22. Berlin, New York: Springer-Verlag. ISBN 978-0-387-97999-1.


\bibitem{numerical-methods-book} Machas, L. "Numerical Methods for Engineers." Available at: \url{https://www.math.hkust.edu.hk/~machas/numerical-methods-for-engineers.pdf}

\bibitem{scipedia1} Tinoco-Guerrero, J. M., et al. (2022). "Numerical Methods for Solving Partial Differential Equations in Finance." Available at: \url{https://www.scipedia.com/public/Tinoco_Guerrero_et_al_2022b}

\bibitem{scholarpedia} Wikipedia. "Method of Lines." Available at: \url{http://www.scholarpedia.org/article/Method_of_lines\#Initial_and_Boundary_Conditions}

\bibitem{libretexts} LibreTexts. "Sparse Matrices." Available at: \url{https://phys.libretexts.org/Bookshelves/Mathematical_Physics_and_Pedagogy/Computational_Physics_(Chong)/08\%3A_Sparse_Matrices/8.02\%3A_Sparse_Matrix_Formats}

\bibitem{superlu} NERSC. "SuperLU." Available at: \url{https://portal.nersc.gov/project/sparse/superlu/}

\bibitem{scipy-sparse} Scipy Lectures. "Solving sparse linear systems with SciPy." Available at: \texttt{https://scipy-lectures.org/advanced/scipy\_sparse/solvers.html}
\bibitem{golubvanloan} Golub, G. H., \& Van Loan, C. F. (2013). Matrix Computations. Johns Hopkins University Press.

\bibitem{itermethbook} Saad, Y. (2003). Iterative Methods for Sparse Linear Systems (2nd ed.). Society for Industrial and Applied Mathematics.
\end{thebibliography}

\pagebreak
\section{Додатки}

\subsection{Функції які реалізують різні методи}
\lstinputlisting[language=Python]{lu_dense_solver.py}

\lstinputlisting[language=Python]{lu_sparse_solver.py}

\lstinputlisting[language=Python]{jacobi_solver.py}

\lstinputlisting[language=Python]{gauss_seidel_solver.py}

\lstinputlisting[language=Python]{gmres_solver.py}

\lstinputlisting[language=Python]{bicgstab_solver.py}

\end{document}