File contents

\documentclass[12pt]{article}
\usepackage{graphicx,amsmath,amsfonts,amssymb,epsfig,euscript,enumerate}
\usepackage[T1]{fontenc}
\usepackage[utf8x]{inputenc}

\newtheorem{ejer}{Ejercicio}
\newcommand{\bej}{\begin{ejer}\rm}
\newcommand{\fej}{\end{ejer}}

\newcommand{\gO}{\Omega}
\newcommand{\R}{\mathbb{R}}
\newcommand{\C}{\mathbb{C}}
\def\dt{\Delta t}
\def\dx{\Delta x}

\topmargin-2cm \vsize 29.5cm \hsize 21cm
\setlength{\textwidth}{16.75cm}\setlength{\textheight}{23.5cm}
\setlength{\oddsidemargin}{0.0cm}
\setlength{\evensidemargin}{0.0cm}

\begin{document}
\centerline{{\small Universidad de Buenos Aires - Facultad de Ciencias Exactas y Naturales - Depto. de Matemática}}

 \vskip 0.2cm
 \hrulefill
 \vskip 0.2cm

 \centerline{{\bf\Huge {\sc Análisis Numérico}}}
 \vskip 0.2cm
 \centerline{\ttfamily Segundo Cuatrimestre 2019}
 \hrulefill

 \bigskip
 \centerline{\bf TP  N$^\circ$ 2 - Tema 1: Campo eléctrico}
 \bigskip


\section{El problema:}

Considerar el  problema de Poisson:

Dada $f:\R^2\longrightarrow\R$, se quiere hallar $u:\R^2\longrightarrow\R$ tal que:
\begin{equation}\label{eq}
-\Delta u = f(x,y) \quad\quad (x,y)\in \Omega = B(0,1)
\end{equation}
con condiciones de Dirichlet homogéneas. 

Físicamente $f$ representa la densidad de carga eléctrica, mientras que la incógnita $u$ es el potencial eléctrico. 
Nos interesa calcular tanto $u$ como el campo eléctrico generado por la distribución de cargas $f$, que viene dado por $E:\R^2\longrightarrow\R^2$,
$$E = -\nabla u$$

El objetivo de este trabajo es implementar un algoritmo de elementos finitos para resolver este problema y estimar el orden de convergencia. 

\section{Formulación Discreta}
Dada una triangulacion $\mathcal{T}_h$ de $\Omega$, considerar el espacio $V_h$ de funciones continuas y lineales en cada triángulo. Llamemos $\{\phi_i\}_i$ a la base nodal de $V_h$. La discretización de la formulación débil del problema conduce a la necesidad de calcular la matriz  de rigidez $R$ y el vector $b$, dados por:
$$R_{i,j} =\int_\Omega\nabla\phi_i\nabla\phi_j,\quad b_i = \int_{\Omega}f\phi_i$$

El problema discreto se reduce a resolver $R\cdot U = b,$ donde $U$ son los coeficientes de $u_h$.
En nuestro caso nos interesa también calcular $\nabla u$, para lo cual necesitaremos una función que tome como datos $U$ y la triangulación y calcule, en cada nodo, el valor de $\nabla u_h$. Para ello, se sugiere seguir una lógica similar a la de la construcción de la matriz de rigidez. Es decir: hacer un ciclo sobre triángulos y calcular el gradiente de $u_h$ en cada nodo localmente. Estos valores se suman y al finalizar el valor acumulado sobre cada nodo se divide por el número de triángulos incidentes sobre el nodo. Así, se calcula, en cada nodo el promedio de los gradientes de $u_h$ sobre los triángulos incidentes. 

\section{Error}
Se quiere también estimar el error de aproximación en $L^2$. Para ello es necesario calcular:
$$\left(\int_{\Omega_h}(u-u_h)^2\right)^\frac{1}{2} = \sum_{T} \left(\int_T (u-u_h)^2\right)^2.$$
Donde $\Omega_h$ es el polígono mallado, que aproxima a $\Omega$. Es importante observar que hay otro error involucrado en la resolución, que es el error de aproximación entre $\Omega_h$ y $\Omega$. Utilizando elementos lineales y siendo $\partial\Omega$ suave, esto no afectará el orden de convergencia, por lo que no lo tomaremos en cuenta. Sin embargo, cuando se trabaja con elementos de mayor orden es necesario incluir este error en el análisis. 

Para calcular el error local, utilizaremos la cuadratura:
$$\int_T f(x) dx \sim |T|(-27f(q_0)+25f(q_1)+25f(q_2)+25f(q_3))/48,$$
donde los nodos pueden calcularse como: 
$$\renewcommand\arraystretch{1.3}\begin{pmatrix} q_0 \\ q_1 \\ q_2 \\ q_3\end{pmatrix} = 
\begin{pmatrix}\frac13 & \frac13 \\ \frac15 &\frac35 \\ \frac15 &\frac15 \\ \frac35 &\frac15\end{pmatrix}\cdot \begin{pmatrix} p_1-p_0\vspace{1pt}&\vspace{2pt}\vdots&p_2-p_0\vspace{1pt} \end{pmatrix} + \begin{pmatrix}p_0\\ p_0\\ p_0 \\ p_0\end{pmatrix}.$$
siendo $p_0,p_1,p_2$ los vértices del triángulo. 

Es importante tener en cuenta que como $u_h$ es lineal en $T$, los valores $u_h(q_i)$ pueden calcularse sencillamente:
$$\renewcommand\arraystretch{1.3}\begin{pmatrix} u_h(q_0) \\ u_h(q_1) \\ u_h(q_2) \\ u_h(q_3)\end{pmatrix} = 
\begin{pmatrix}\frac13 & \frac13 \\ \frac15 &\frac35 \\ \frac15 &\frac15 \\ \frac35 &\frac15\end{pmatrix}\cdot \begin{pmatrix} \vspace{2pt}u_h(p_1)-u_h(p_0)&\vdots& \vspace{2pt}u_h(p_2)-u_h(p_0) \end{pmatrix} + \begin{pmatrix}u_h(p_0)\\ u_h(p_0)\\ u_h(p_0)\\ u_h(p_0)\end{pmatrix}.$$

\section{Datos}
Resolver \eqref{eq} para: 
\begin{itemize}
\item $f(x,y) = 4$, que da por solución exacta: $u(x,y) = 1-x^2-y^2$
\item para una fuente puntual $f$ tal que $\int f v = v(0,0)$, que da por solución exacta: $u(x,y) = \frac{1}{2\pi}\log(\sqrt{x^2+y^2}).$ 
\end{itemize}

Resolver para distintos valores de $h$ (por ejemplo: $h=2^{-2},2^{-3},2^{-4}$) y graficar, en escala logarítmica, el error en función de $\frac{1}{N^2}$, donde $N$ el número de vértices de la malla. 

\subsection{Comandos útiles:}
Pueden resultar útiles las siguientes funciones: 
\begin{itemize}
\item \verb+Delaunay+, de la biblioteca \verb+scipy.spatial+ para la construcción de la triangulación. Alternativamente, puede usarse \verb+Triangulation+ de \verb+matplotlib.tri+. 
\item \verb+csr_matrix+ de la biblioteca \verb+scipy.sparse+ para el ensamblado de la matriz de rigidez. 
\item \verb+spsolve+ o \verb+cg+ de \verb+scipy.sparse.linalg+ para la resolución de sistemas.
\item \verb+plot_trisurf+, de \verb+matplotlib.pyplot+ para el gráfico del potencial $u$ y \verb+quiver+ para el gráfico de $-\nabla u$. Preferentemente, graficar ambos en una misma figura, con dos subplots.  
\item \verb+loglog+ de \verb+matplotlib.pyplot+ para gráficos en escala logarítmica. 
\end{itemize}
\end{document}