Solução estacionária da EDP de Defeitos de Radiação por métodos numéricos
Forfatter
Davi
Last Updated
8 years ago
License
Creative Commons CC BY 4.0
Resumé
It's a work that discuss solutions of an specific equation using different numerical methods.
It's a work that discuss solutions of an specific equation using different numerical methods.
\documentclass[DIV=calc, paper=a4, fontsize=12pt, onecolumn]{scrartcl} % A4 paper and 12pt font size
\usepackage[brazil]{babel} % English language/hyphenation
\usepackage[protrusion=true,expansion=true]{microtype} % Better typography
\usepackage{amsmath,amsfonts,amsthm} % Math packages
\usepackage[svgnames]{xcolor} % Enabling colors by their 'svgnames'
\usepackage[hang, small,labelfont=bf,up,textfont=it,up]{caption} % Custom captions under/above floats in tables or figures
\usepackage{booktabs} % Horizontal rules in tables
\usepackage{fix-cm} % Custom font sizes - used for the initial letter in the document
\usepackage[utf8]{inputenc}
\usepackage{sectsty} % Enables custom section titles
\allsectionsfont{\usefont{OT1}{phv}{b}{n}} % Change the font of all section commands
\usepackage{graphicx}
\usepackage{fancyhdr} % Needed to define custom headers/footers
\pagestyle{fancy} % Enables the custom headers/footers
\usepackage{lastpage} % Used to determine the number of pages in the document (for "Page X of Total")
% Headers - all currently empty
\lhead{}
\chead{}
\rhead{}
% Footers
\lfoot{}
\cfoot{}
\rfoot{\footnotesize Página \thepage\ de \pageref{LastPage}} % "Page 1 of 2"
\renewcommand{\headrulewidth}{0.0pt} % No header rule
\renewcommand{\footrulewidth}{0.4pt} % Thin footer rule
\usepackage{lettrine} % Package to accentuate the first letter of the text
\newcommand{\initial}[1]{ % Defines the command and style for the first letter
\lettrine[lines=3,lhang=0.3,nindent=0em]{
\color{DarkGoldenrod}
{\textsf{#1}}}{}}
%----------------------------------------------------------------------------------------
% TITLE SECTION
%----------------------------------------------------------------------------------------
\usepackage{titling} % Allows custom title configuration
\newcommand{\HorRule}{\color{DarkGoldenrod} \rule{\linewidth}{1pt}} % Defines the gold horizontal rule around the title
\pretitle{\vspace{-90pt} \begin{flushleft} \HorRule \fontsize{25}{25} \usefont{OT1}{phv}{b}{n} \color{DarkRed} \selectfont} % Horizontal rule before the title
\title{Solução estacionária da EDP de Defeitos de Radiação por métodos numéricos} % Your article title
\posttitle{\par\end{flushleft}\vskip 0.5em} % Whitespace under the title
\preauthor{\begin{flushleft}\large \lineskip 0.5em \usefont{OT1}{phv}{b}{sl} \color{DarkRed}} % Author font configuration
\author{Davi Lazzari, } % Your name
\postauthor{\footnotesize \usefont{OT1}{phv}{m}{sl} \color{Black} % Configuration for the institution name
(232847) - Métodos Computacionais da Física C - UFRGS % Your institution
\par\end{flushleft}\HorRule} % Horizontal rule after the title
\date{Porto Alegre, 27 de novembro de 2015} % Add a date here if you would like one to appear underneath the title block
%----------------------------------------------------------------------------------------
\begin{document}
\maketitle % Print the title
\thispagestyle{fancy} % Enabling the custom headers/footers for the first page
%----------------------------------------------------------------------------------------
% ARTICLE CONTENTS
%----------------------------------------------------------------------------------------
\section*{Trabalho}
O trabalho consiste em comparar a solução da equação de difusão de Defeitos de Radiação para diferentes métodos numéricos e diferentes condições de contorno. O que se busca, mais especificamente, é comparar a solução estacionária da equação.
Equação da difusão de defeitos por radiação:
\begin{equation}
{{\partial \pmb{f} \over \partial t} = D{\partial^2 \pmb{f} \over \partial x^2} - S\pmb{f} + \pmb{s}(x)}
\end{equation}
\section*{Condições de Contorno: $f_0=0$ e $f_L=0$}
\subsection*{Método de Thomas}
O método de Thomas é um método numérico para encontrar-se uma solução na qual a transformação aplicada resulta em algo conhecido. Basicamente o que se faz é inverter uma matriz sem realmente precisar invertê-la.
Como busca-se a solução no ponto de equilíbrio, tem-se que o sistema encontra-se em um estado estacionário, logo, $f(x+\Delta X) \approx f(x) \to \partial f / \partial x = 0$. Assim:
\begin{equation}
{D{\partial^2 \pmb{f} \over \partial x^2} - S\pmb{f} + \pmb{s}(x) = 0 \to \left(D{\partial^2 \over \partial x^2} - S\right)\pmb{f} = -\pmb{s}(x)}
\end{equation}
o primeiro termo da equação é um operador de derivação subtraido de um sumidouro (\textit{sink}) aplicado em um vetor. O que se tem na verdade é o termo de fonte $\vec{s}(\vec{x})$ que é dado por uma gaussiana. Quando abre-se o termo para o espaço discreto, tem-se $\partial x \to \Delta x$ e a respectiva equação (omitindo-se o termo $-S$ e $D$):
\begin{equation}
{{\partial^2 f \over \partial x^2} \approx {f(x-\Delta x, t) -2f(x,t) + f(x+\Delta x,t) \over (\Delta x)^2}}
\end{equation}
que em forma matricial, apresenta elementos apenas na diagonal principal, na diagonal abaixo (elementos de $x-\Delta x$) e na diagonal acima da principal (elementos de $x+\Delta x$). Por escalonamento, resolvem-se o sistema, mantendo $\vec{f}(\vec{x})$ como incógnita em função de $\vec{s}(\vec{x})$.
E a solução da EDP é apresentada na \textit{Figura 1}.
\begin{figure}[!h]
\includegraphics[width=9cm,angle=-90,trim=0cm 0 0 0cm,clip]{thomas.eps}
\caption{Solução estacionária para a equação de Difusão de Defeitos por Radiação pelo Algoritmo de Thomas}
\label{fig1}
\end{figure}
O programa utilizado para gerar-se esta solução tem como condições para a fonte: $s_0=10$, $\sigma = 10$ e centrada em $x_0 = 100$; para o sumidouro utilizou-se $S=0.0002$ em uma rede de 400 pontos.
\subsection*{Foward Time Central Space (FTCS)}
Este método baseia-se na discretização do tempo e do espaço, assim o que ele faz não é uma solução única para o caso estacionário, mas sim a evolução temporal da solução. Discretiza-se o tempo então através do método de Euler, onde avança-se $\Delta t$ no tempo:
\begin{equation}
{{\partial f \over \partial t} \approx {f(x,t+\Delta t) -f(x,t) \over \Delta t}}
\end{equation}
Substitui-se então as \textit{Equações 3} e \textit{4} na \textit{Equação 1} e assumindo uma notação mais compacta, obtém-se:
\begin{equation}
{{f_{j}^{n+1} -f_{j}^{n} \over \Delta t} = D{f_{j-1}^{n} -2f_{j}^{n} + f_{j+1}^{n} \over (\Delta x)^2} -S\pmb{f} +\pmb{s}(x)}
\end{equation}
isolando-se o $f_{j}^{n+1}$
\begin{equation}
{f_{j}^{n+1} = f_{j}^{n} + {D\Delta t \over (\Delta x)^2} (f_{j-1}^{n} -2f_{j}^{n} + f_{j+1}^{n}) -S\pmb{f} +\pmb{s}(x)}
\end{equation}
Assim, para encontrar a solução estacionária, itera-se a equação para um tempo suficientemente grande. Iterando-se para tempo igual 30000, com as condições iniciais já citadas, obtém-se então a solução da \textit{Figura 2} que se assemelha muito com a obtida pelo método de Thomas.
\begin{figure}[!hp]
\includegraphics[width=9cm,angle=-90,trim=0cm 0 0 0cm,clip]{ftcs.eps}
\caption{Solução para a equação de Difusão de Defeitos por Radiação pelo Método FTCS para tempo = 30000 (estacionária)}
\label{fig2}
\end{figure}
\section*{Método de Crank-Nicolson}
O método de Crank-Nicolson baseia-se em uma organização algébrica das equações, de modo a se acompanhar a evolução temporal das soluções. Obtém-se a seguinte equação:
\begin{equation}
{\pmb{f}^{n+1} = \pmb{E}^{-1} ({\pmb{M}. \pmb{f}^{n} +dt.\pmb{s}})}
\end{equation}
Onde $\pmb{E}$ e $\pmb{M}$ são matrizes tridiagonais bem determinadas, com seus termos constantes dependentes de $\Delta t$, $\Delta x$, $D$ e $S$.
Nota-se então que o termo entre parenteses no lado esquerdo é um vetor, facilmente determinado pela multiplicação da matriz $\pmb{M}$ pelo vetor do tempo atual da iteração $\pmb{f}^{n}$ somado ao vetor de fonte. Definindo este vetor como um $\pmb{\psi}^{n}$, tem-se que:
\begin{equation}
{\pmb{f}^{n+1} = \pmb{E}^{-1}\pmb{\psi}^{n} \to \pmb{E}.\pmb{f}^{n+1} = \pmb{\psi}^{n}}
\end{equation}
como a matriz tridiagonal $\pmb{E}$ e o vetor $\pmb{\psi}^{n}$ são conhecidos, busca-se o vetor $\pmb{f}^{n+1}$; mas, observa-se ainda que a forma da equação é a mesma que fora resolvida na \textit{Equação 2}, ou seja, pode-se resolver a partir do método de Thomas.
Infelizmente, o programa para o cálculo da equação de Difusão de Defeitos por Radiação pelo método de Crank-Nicolson não funcionou como deveria. O que se teve foi um comportamento para temos iniciais próximo ao do método FTCS, porém, não se chegou em uma solução estacionária, o que ocorreu foi uma divergência por parte da equação.
\section*{Análise de Estabilidade}
Deve-se falar em estabilidade apenas para os dois métodos que demandavam uma evolução temporal, visto que para a solução estacionária apenas uma iteração resolvia a equação.
O método FTCS é um método muito bom pela facilidade de sua implementação e pela não obrigatoriedade no uso de matrizes e principalmente na inversão delas. Entretanto, o método demanda que o termo multiplicativo $k \equiv D\Delta t / (\Delta x)^2$ seja menor que $1/2$ e assim, a amplitude do módo de fourier associado não divirja, o que determina que $dt$ seja muito pequeno, e consequentemente, a convergência se dê de forma muito lenta.
Já para o método de Crank-Nicolson este problema não ocorre, ele é organizado de tal maneira em que o $dt$ seja basicamente um parâmetro livre. Assim, por mais que se demande muito mais tempo computacional para resolver as equações matriciais, os passos de tempo podem ser muito maiores, auxiliando numa rápida convergência do método.
Assim, pode-se comparar na \textit{Figura 3} as duas soluçõe obtidas pelos dois primeiros métodos. Percebe-se que o pico da solução para o método FTCS é menor que para a obtida pelo algoritmo de Thomas, supõe-se assim que essa diferênça seja devida ao tempo utilizado, não suficiente para se chegar no estado estacionário.
\begin{figure}[!h]
\begin{center}
\includegraphics[width=9cm,angle=-90,trim=0cm 0 0 0cm,clip]{juntas.eps}
\caption{Tipos de organização de discos para ocupação do espaço.}
\label{fig3}
\end{center}
\end{figure}
\section*{Condições de Contorno: $f_0=f_L$}
\subsection*{Método de Sherman-Morrison}
O método de Sherman-Morrison é um método que se aproxima do algoritmo de Thomas no sentido de funcionar para inverter uma matriz, com a grande diferença de que funciona para uma matriz qualqer, não apenas para matrizes tridiagonais.
A questão é que ao se usar condições de contorno periódicas, surgem dois termos nas diagonais opostas da matriz, e então há de se resolver o sistema para esta nova matriz.
O algoritmo de Sherman-Morrison baseia-se na consideração de uma matriz como uma combinação de matrizes com vetores, o que seria equivalente a um método de perturbação de uma matriz. Fala-se assim até em \textit{ranks} (ordens) de perturbação para a matriz.
Basicamente, a forma proposta pelo algoritmo é a seguinte:
\begin{equation}
{(\pmb{A}+\pmb{u}\pmb{v}^{T})^{-1} = \pmb{A}^{-1} - {\pmb{A}^{-1}\pmb{u}\pmb{v}^{T}\pmb{A}^{-1} \over 1+\pmb{v}^{T}\pmb{A}^{-1}\pmb{u}}}
\end{equation}
onde $\pmb{u}$ e $\pmb{v}$ são vetores quaisquer, $\pmb{A}$ é uma matriz inversível e $1+\pmb{v}^{T}\pmb{A}^{-1}\pmb{u} \neq 0$.
\end{document}