\documentclass[12pt]{article}
%\include{amsfonts}
\usepackage{amssymb,amsmath,latexsym,epsfig,euscript,multicol}
\usepackage{enumitem}
\usepackage[utf8x]{inputenc}
\usepackage[spanish]{babel} 

% Caracteres especiales
\def\A{\mathbb{A}}
\def\C{\mathbb{C}}
\def \N{\mathbb{N}}
\def \P{\mathbb{P}}
\def \Q{\mathbb{Q}}
\def \R{\mathbb{R}}
\def \Z{\mathbb{Z}}
\def \sen{\textrm{sen}}

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

\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}

\usepackage{chngcntr}

\counterwithin*{equation}{section}
\counterwithin*{ejer}{section}

\begin{document}


\begin{center}
 \textbf{Pautas para el trabajo práctico}
\end{center}

\begin{enumerate}
 \item Formar grupos de dos o tres personas.
 \item Enviar un e-mail a la dirección \texttt{ecn.octave@gmail.com} antes del lunes 7 de noviembre a las 13 horas con los datos de los integrantes del grupo (nombre y número de libreta) y con el número de enunciado elegido (que no podrá ser modificado).
 \item La defensa del trabajo es presencial, en un laboratorio u otro lugar a determinar, y será entre los días 15 y 22 de noviembre, inclusive, preferentemente en horario de la práctica de la materia. No hay que entregar código impreso. Se debe enviar por e-mail a la dirección \texttt{ecn.octave@gmail.com} el conjunto de archivos del programa que se pide en el enunciado, con al menos 24 horas de anticipación a la defensa. De ser necesario, se puede mostrar en papel todo aquello que complemente o explique la resolución. Haremos preguntas a cada integrante del grupo. 
\end{enumerate}


\newpage

\clearpage
\setcounter{page}{1}
\setcounter{ejer}{0}
\renewcommand{\theenumi}{\textbf{\arabic{enumi}}}





%----------------------------------------------------------
%Encabezado:

\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\Large{\sc Elementos de Cálculo Numérico / Cálculo Numérico}}}
 \vskip 0.2cm
 \centerline{\ttfamily Segundo Cuatrimestre 2016}
 \hrulefill

 \bigskip
 \centerline{\bf Trabajo Práctico N$^\circ$ 1: Curva Dosis-Respuesta.}
 \bigskip

 
 
Se tiene el siguiente sistema de reacciones químicas, que representa el sistema bacterial EnvZ-OmpR:
 
\[ \text{X}\overset{k_1}{\underset{k_2}{\rightleftarrows}} \text{XT}  \overset{k_3}{\longrightarrow}  \text{X}_p\]
\[ \text{X}_p+\text{Y}\overset{k_4}{\underset{k_5}{\rightleftarrows}} \text{X}_p\text{Y}  \overset{k_6}{\longrightarrow}  \text{X}+\text{Y}_p\]
\[ \text{XT}+\text{Y}_p\overset{k_7}{\underset{k_8}{\rightleftarrows}} \text{XT}\text{Y}_p  \overset{k_9}{\longrightarrow}  \text{XT}+\text{Y}_p\]
 
Se nota entre corchetes a la concentración de cada especie química, y se las ordena de la siguiente manera: $x_1=[\text{X}]$, $x_2=[\text{XT}]$, $x_3=[\text{X}_p]$, $x_4=[\text{Y}]$, $x_5=[\text{Y}_p]$, $x_6=[\text{X}_p\text{Y}]$, $x_7=[\text{XTY}_p]$. Las ecuaciones del sistema bajo cinética de acción de masas son:

\begin{align*}
 f_1=\frac{dx_1}{dt} &= -k_1x_1+k_2x_2+k_6x_6, \\
 f_2=\frac{dx_2}{dt} &= k_1x_1-(k_2+k_3)x_2-k_7x_2x_5+(k_8+k_9)x_7, \\
 f_3=\frac{dx_3}{dt} &= k_3x_2-k_4x_3x_4+k_5x_6,\\
 f_4=\frac{dx_4}{dt} &= -k_4x_3x_4+k_5x_6+k_9x_7, \\
 f_5=\frac{dx_5}{dt} &= k_6x_6-k_7x_2x_5+k_8x_7, \\
 f_6=\frac{dx_6}{dt} &= k_4x_3x_4-(k_5+k_6)x_6,\\
 f_7=\frac{dx_7}{dt} &= k_7x_2x_5-(k_8+k_9)x_7.\\
\end{align*}

Notemos que $f_1+f_2+f_3+f_6+f_7=0$ y $f_4+f_5+f_6+f_7=0$, por lo que existen constantes $C_1$ y $C_2$ tales que $x_1(t)+x_2(t)+x_3(t)+x_6(t)+x_7(t)=C_1$ y $x_4(t)+x_5(t)+x_6(t)+x_7(t)=C_2$ para todo $t\ge 0$. En particular, $x_1(0)+x_2(0)+x_3(0)+x_6(0)+x_7(0)=C_1$ y $x_4(0)+x_5(0)+x_6(0)+x_7(0)=C_2$.

Se desea estudiar la concentración en equilibrio de ciertas especies del sistema a medida que las cantidades (totales) $C_1$ y $C_2$ aumentan. Se proponen las siguientes constantes de reacción:
$k_1=0.1; k_2=0.2; k_3=0.01; k_4=0.01; k_5=0.01; k_6=0.1; k_7=0.02; k_8=0.2; k_9=0.05$.

\bej
Resolver el sistema numéricamente para distintas cantidades totales $C_1$, tomando como dato inicial al vector $x_1=C_1$, $x_5=C_2=100$, $x_i=0$ para todo $i\neq 1, \; 5$. 

Utilizar los siguientes métodos numéricos:
\begin{enumerate}
 \item El método multipaso Adams-Bashforth que se presenta a continuación. Los primeros términos se calculan utilizando el método de Euler visto en la Práctica 2.
 $$y_{n+3}-y_{n+2}=\frac{h}{12}(23f(y_{n+2})-16f(y_{n+1})+5f(y_{n})).$$
 \item El comando \verb+ode45+, o \verb+lsode+.
\end{enumerate}
En cada caso, avanzar hasta un tiempo $T$ lo suficientemente grande como para considerar el sistema en equilibrio. Graficar superpuestos los valores en equilibrio de $x_1$ y $x_5$ contra las distintas cantidades totales de $C_1$, para ambos métodos.  
\fej

\bej
Repetir el ejercicio anterior para distintas cantidades $C_2$, tomando como dato inicial al vector $x_1=C_1=100$, $x_5=C_2$, $x_i=0$ para todo $i\neq 1, \; 5$. Graficar superpuestos los valores en equilibrio de $x_1$ y $x_5$ contra las distintas cantidades totales de $C_2$, para ambos métodos.  
\fej

Notar que los valores en equilibrio de Y$_p$ no varían, sin importar las cantidades iniciales (totales). Se puede ver que $(1+\frac{k_8}{k_9})(f_1+f_5+f_7)+ \frac{k_8}{k_9}f_2=0$ por lo que, en equilibrio, se tiene  que $k_3(1+\frac{k_8}{k_9})x_2-k_7x_2x_5=0$, y si $x_2$ en equilibrio es positivo, entonces $x_5=\frac{k_3(1+\frac{k_8}{k_9})}{k_7}$. Esta constante es independiente de los valores de $C_1$ y $C_2$. Este fenómeno se lo conoce como ``robustez de concentración absoluta'' en la especie Y$_p$ y fue descripto en  \cite{sf10}.



%%=========================================
\begin{thebibliography}{1}
%%=========================================
\bibitem{sf10}
{\sc{G. Shinar y M. Feinberg}},
\emph{Structural sources of robustness in biochemical reaction networks},
{Science} 327(5971), (2010), pp. 1389--1391.
\end{thebibliography}

\newpage

\clearpage
\setcounter{page}{1}
\setcounter{ejer}{0}
\renewcommand{\theenumi}{\textbf{\arabic{enumi}}}
%----------------------------------------------------------
%Encabezado:

\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\Large{\sc Elementos de Cálculo Numérico / Cálculo Numérico}}}
 \vskip 0.2cm
 \centerline{\ttfamily Segundo Cuatrimestre 2016}
 \hrulefill

 \bigskip
 \centerline{\bf Trabajo Práctico N$^\circ$ 2: Oscilaciones de una membrana cuadrada.}
 \bigskip

El problema de la oscilación de una membrana cuadrada fija por sus bordes puede ser descripto, al menos para pequeñas elongaciones, por la ecuación de ondas con condición de contorno:
\begin{equation*}
 \left\lbrace\begin{array}{ll}
              u_{tt}=c^2\Delta u & (x,y) \in \Omega=[-L,L]\times[-L,L]\\
              & u|_{\partial \Omega}=0
             \end{array}\right.
\end{equation*}
Donde $\Delta u = \frac{\partial^2u}{\partial x^2}+\frac{\partial^2u}{\partial y^2}$. Los datos iniciales adecuados son, por ejemplo, la posici\'on inicial de la membrana y su velocidad inicial.\\[10pt]
Para resolver este problema numéricamente, se discretiza el dominio con una malla $(x_i , y_j) = (i/N, j/N), 0 \leq i, j \leq N$ , y el tiempo $t_n = nk$ hasta un $T$ m\'aximo.

Considerando $u^n_{i,j}$, la aproximación de $u(x_i, y_j , t_n)$ (el valor de la función en el nodo $(x_i, y_j , t_n)$ de la malla), se propone el siguiente esquema, que proviene de la discretización usual de la derivada segunda.
\begin{equation}\label{eq:discret}
 u^{n+1}_{i,j}-2u^{n}_{i,j}+u^{n-1}_{i,j}=\dfrac{c^2k^2}{h^2}[u^{n+1}_{i-1,j}+u^{n+1}_{i+1,j}+u^{n+1}_{i,j-1}-4u^{n+1}_{i,j}+u^{n+1}_{i,j+1}]
\end{equation}
donde $h = 2L/N$ . Por lo tanto, suponiendo conocidos la posición de la membrana a tiempo $0$ y a tiempo $−k$, se puede calcular la posición para cualquier $t_n$ posterior. La solución a tiempo cero se conoce por el dato inicial. La condición inicial en la derivada respecto del tiempo permite plantear una ecuación adicional
\[u^{0}_{i,j}-u^{-1}_{i,j}=k\, u_t(x_i , y_j , t_0 )\]
de la cual despejar el valor de $u^{-1}_{i,j}$. Reemplazando este valor en la ecuación \eqref{eq:discret}, para $n = 0$, se obtiene el valor de $u^1$.

\bej Escribir un programa que resuelva el problema de la membrana, con dato de contorno $0$, posición inicial
\[u_0(x, y) = \sin\left(\tfrac{\pi}{2L} (x+L)\right) \sin\left(\tfrac{\pi}{2L} (y+L)\right) \qquad (x, y) \in [-L, L] \times [-L, L],\]
velocidad inicial cero, $L=1/2$ y $c=1$.  Comparar con la solución exacta
\[u(x, y, t) = \cos\left(\tfrac{\pi}{\sqrt{2}L} t\right) \sin\left(\tfrac{\pi}{2L} (x+L)\right) \sin\left(\tfrac{\pi}{2L} (y+L)\right).\]
\fej

\bej Se considera ahora una nueva condición inicial. Escribir un programa que resuelva el problema de la membrana, con dato de contorno $0$, posición inicial
\[u_0(x, y) = -e^{-\alpha(x^2+y^2)} \qquad (x, y) \in [-L, L] \times [-L, L],\]
velocidad inicial cero, $\alpha=12$, $L=5/2$ y $c=10$. ¿Cómo varía la dinámica si $\alpha \in [1,20]$? 
\fej

\bej  Usando los comandos \verb+drawnow+ y \verb+pause(+$\cdot$\verb+)+ obtener una \emph{película} mostrando la evolución de las soluciones numéricas obtenidas en los dos ítems anteriores.
\fej

\newpage

\clearpage
\setcounter{page}{1}
\setcounter{equation}{0}
\setcounter{ejer}{0}
\renewcommand{\theenumi}{\textbf{\arabic{enumi}}}
%----------------------------------------------------------
%Encabezado:

\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\Large{\sc Elementos de Cálculo Numérico / Cálculo Numérico}}}
 \vskip 0.2cm
 \centerline{\ttfamily Segundo Cuatrimestre 2016}
 \hrulefill

 \bigskip
 \centerline{\bf Trabajo Práctico N$^\circ$ 3: Resonancia de la cuerda vibrante.}
 \bigskip

 
El problema de la oscilación de una cuerda fija en sus extremos, con un forzante $f$, puede ser descripto, al menos para pequeñas elongaciones, por
la ecuación de ondas con condición de contorno:
\begin{equation}\label{eq:osc}
 \left\lbrace\begin{array}{rl}
              u_{tt}(x,t)= \alpha u_{xx}(x,t) +f(x,t) & x \in \Omega=[0,1], t>0\\
               u(0)=u(1) & =0
             \end{array}\right.
\end{equation}
Los datos iniciales adecuados son, por ejemplo, la posición inicial de la cuerda y su velocidad inicial. Sin p\'erdida de generalidad supondremos que $\alpha = 1$. 

Resolviendo el problema con $f=0$, por el método de separación de variables, se buscan soluciones de la forma
$U (t, x) = T (t)Z(x)$; de modo que
\[\dfrac{T''(t)}{T(t)}=\dfrac{Z''(x)}{Z(x)}=\lambda,\]
para alguna constante $\lambda$. Usando las condiciones de contorno, tenemos que
\[Z(0)=Z(1)= 0\]
es decir que $\lambda$ es un autovalor de la ecuación de Laplace con datos de
contorno Dirichlet homogéneo
\begin{equation}\label{eq:Z}
 \left\lbrace\begin{array}{ll}
              Z'' -\lambda Z =0 & x \in [0,1]\\
               Z(0)=Z(1) &=0
             \end{array}\right.
\end{equation}
Se sabe que la ecuación \eqref{eq:Z} tiene una sucesión de autovalores negativos $\{\lambda_n\}_n$ con correspondientes autofunciones $Z_n(x)$ y que $\lambda_n \to -\infty$ cuando $n \to +\infty$. La correspondiente función $T_n$ es una solución de la ecuación
\[T'' − \lambda_n T = 0\]
Si $\lambda_n = -\omega_n^2$ resulta que $T_n(t) = a_n \cos(\omega_nt) + b_n \sin(\omega_nt)$, y de este modo, las soluciones de \eqref{eq:osc} son de la forma
\begin{equation}\label{eq:u}
 u(t,x)=\underset{n=1}{\overset{\infty}\sum}(a_n\cos(\omega_nt)+b_n\sin(\omega_nt))Z_n(x)
\end{equation}
donde los valores $a_n$ y $b_n$ se obtienen a partir de las condiciones iniciales. A los autovalores $\lambda_n$ se los llama los modos \textit{normales de oscilación}.


\bej Calcular los autovalores y las autofunciones del laplaciano. Para ello, considerar una discretización del dominio $\Omega = [0, 1]$ con una malla $x_i = i/N$, $0 \leq i\leq N$ . 

Llamando $h= 1/N$ se obtiene el siguiente esquema, sobre la malla introducida, para aproximar la ecuación \eqref{eq:Z}. 
\[\dfrac{z_{i+1}-2z_{i}+z_{i-1}}{h^2}=\lambda z_i \quad i=1,\dots,N-1.\]
Reescribir esto como un problema lineal $Az = \lambda z$. Y resolver numéricamente. Graficar los 9 primeros autovectores (aproximaciones de las autofunciones). (Consultar la documentación de los comandos \verb+eig+, \verb+eigs+ y \verb+subplot+).
\fej


Finalmente, se desea resolver el problema \eqref{eq:osc}, considerando datos iniciales: 
$$\begin{array}{cc} u(x,0) = g(x),\\
  u_t(x,0) = 0.\end{array},$$
y forzante dado por una función periódica de la forma: $\cos(\omega t)v(x)$. Cuando la frecuencia de oscilación $\omega$ se corresponde con un autovalor del laplaciano se produce el fenómeno de resonancia: como la frecuencia del forzante coincide con una frecuencia propia de la cuerda, la oscilación se hace cada vez más grande.

Discretizando la ecuación en un esquema implícito se obtiene:
$$\frac{u^{n+1}_{j}-2u^{n}_j+u^{n-1}{j}}{dt^2} = \frac{u_{j+1}^{n+1}-2u^{n+1}_j+u^{n+1}_{j-1}}{h^2} + f_{j}^n,$$
donde $u_j^n$ es la aproximación de $u(x_j,t_n)$ y $f^n_j =f(x_j,t_n)$. 
La solución a tiempo cero se conoce por el dato inicial. La condición inicial en la derivada respecto del tiempo permite plantear una ecuación adicional
\[u^{0}_j-u^{-1}_j=dt \,u_t(x_j , t_0 )\]
de la cual despejar el valor de $u^{-1}_j$.

\bej Teniendo en cuenta las condiciones de contorno, escribir el sistema lineal que debe resolverse en cada iteración para calcular $u^{n+1}$ en función de $u^n$, $u^{n-1}$ y $f^{n}$. Implementar un programa que resuelva el problema y grafique la evolución de la solución en función del tiempo, aprovechando los comandos \verb+pause+ y \verb+drawnow+.

Resolver considerando $f(x,t) = \cos(\sqrt{|\lambda|}t)v(x)$, donde $\lambda$  y $v$ son un autovalor del laplaciano y su correspondiente autovector, calculados con el programa del ejercicio anterior. En principio considerar $g=0$. Observar qué sucede, por ejemplo, cuando se toma $g$ igual al primer autovector del laplaciano (multiplicado por $0.1$) y $v$ el tercer autovector. Resolver también con dato inicial $g(x) = 0.1x(1-x)$.
\fej


\newpage

\clearpage
\setcounter{page}{1}
\setcounter{equation}{0}
\setcounter{ejer}{0}
\renewcommand{\theenumi}{\textbf{\arabic{enumi}}}

%----------------------------------------------------------
%Encabezado:

\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\Large{\sc Elementos de Cálculo Numérico / Cálculo Numérico}}}
 \vskip 0.2cm
 \centerline{\ttfamily Segundo Cuatrimestre 2016}
 \hrulefill

 \bigskip
 \centerline{\bf Trabajo Práctico N$^\circ$ 4: Obtención de la curva braquistocrona.}
 \bigskip

 
El problema de la braquistocrona fue propuesto en 1696 por Johann Bernoulli. Consiste en encontrar la curva, digamos la forma del tobogán, que une dos puntos de manera que los cuerpos que caen por ella lo hagan en el menor tiempo posible. Una primera conjetura que se podría hacer es que la braquistocrona es una recta, que es la curva más corta que une dos puntos. La solución correcta, sin embargo, es un arco de cicloide, la curva descrita por un punto en el borde de una moneda que rueda a lo largo de una regla. Concretamente es un arco que parte de uno de sus puntos picudos. 

El objetivo de este trabajo práctico es realizar una aproximación numérica de la curva. Para esto, fijemos primero unas ideas del problema. Digamos que inicialmente nuestra partícula se halla en reposo en el punto $(0,1)$ donde $1$ es la altura en metros y debe llegar al punto $(L,0)$, moviéndose sobre el vínculo $y = y(x)$ sometida únicamente a la fuerza de gravedad. Consideramos que todo rozamiento resulta despreciable.  

Se sabe que unas de las características de este movimiento es la conservación de la energía, por lo que se puede afirmar que $$ E = \frac{m}{2} (\dot{x}^2 + \dot{y}^2 ) + mgy$$

Donde $m$ es la masa de la partícula, $g = 9.8 \; mts/seg^2$ es la aceleración de la gravedad y $\dot{x} = \frac{dx}{dt} \; , \; \dot{y} = \frac{dy}{dt} $.

Luego como $E$ es constante, en particular para $t=0$ se tiene que $E = mg$. Observemos que además vale que
$$ \dot{y} =  \frac{dy}{dt} = \frac{dy}{dx} \frac{dx}{dt} = y' \; \dot{x} $$

Concluimos que 

$$ mg = \frac{m}{2} \dot{x}^2 (1 + {y'}^2) + mgy $$

De la ecuación de arriba podemos despejar $\dot{x}$ llegando a que

$$ \dot{x} = \frac{dx}{dt} = \sqrt{\frac{2g(1-y(x))}{1+y'(x)^2}} $$

Finalmente 

$$ dt = \sqrt{\frac{1+{y'(x)}^2}{2g(1-y(x))}} \; dx $$

El tiempo que tardará nuestra partícula en recorrer toda la curva $y = y(x)$ para llegar de un extremo a otro será entonces igual a 
$$ \emph{F}(y) = \int_{0}^{L}  \sqrt{\frac{1+{y'(x)}^2}{2g(1-y(x))}} \; dx $$ 

Donde se asume que $\emph{F}: C^2(0,L) \rightarrow \R $. 

Para resolver este problema consideramos una discretización del intervalo $[0,L]$ dado por $x_0=0,x_1=h,\dots,x_n=L-h,x_{n+1}=L$, para algún parámetro $h$. Notamos $y_i$ a la aproximación de la curva $y(x)$ en el punto $x_i$. Sabemos que $y_0=1$ e $y_{n+1}=0$. Las incógnitas del problema son $y_1,\dots,y_n$. Dada esta discretización, si asumimos que $y$ es lineal en cada intervalo $[x_i,x_{i+1}]$, se tiene que:
\begin{align*}
F(y) &= \sum_{i=0}^{n} \int_{x_i}^{x_{i+1}} \sqrt{\frac{1+\Big(\frac{y_{i+1}-y_i}{h}\Big)^2}{2g(1-y(x))}}dx = \sum_{i=0}^{n} \int_{x_i}^{x_{i+1}} \sqrt{1+\big(\frac{y_{i+1}-y_i}{h}\big)^2}(2g)^{-\frac{1}{2}} (1-y(x))^{-\frac{1}{2}}dx \\
&= \sum_{i=0}^{n}\sqrt{1+\Big(\frac{y_{i+1}-y_i}{h}\Big)^2}(2g)^{-\frac{1}{2}}\frac{2h}{y_{i+1}-y_i}\big(\sqrt{1-y_i}-\sqrt{1-y_{i+1}}\big) \\
&= -(2/g)^{\frac{1}{2}}\sum_{i=0}^{n}\sqrt{\Big(\frac{h}{y_{i+1}-y_i}\Big)^2+1}\big(\sqrt{1-y_i}-\sqrt{1-y_{i+1}}\big)
\end{align*}

\bej 
Implementar una función \verb+f_braq+ que reciba un vector $y$ de cierta longitud $n$ y calcule el valor del funcional evaluado en $y$: $F(y)$.
\fej


Para resolver el problema, y estimar la braquistocrona, utilizaremos el método de Newton para minimizar $F$. El método de Newton está dado por:
$$y^{k+1} = y^k - H(y^k)^{-1}\nabla F(y^k),$$
donde $H(y)$ es la matriz hessiana de $F$ evaluada en $y$.

\bej Implementar un programa que reciba como input una función $f$ y un punto $z$ y calcule el vector gradiente de $f$ evaluado en $z$, y un programa que calcule el hessiano de $f$ evaluado en $z$. Ambos utilizando diferencias forward.
\fej

\bej
Finalmente, implementar un programa que aplique el método de Newton al funcional $F$ con dato inicial dado por la recta que une los puntos $(0,1)$ y $(L,0)$. Resolver para $L=\pi/2$. Utilizando los comandos \verb+pause+ y \verb+drawnow+, graficar la solución para cada iteración de Newton. Comparar con la solución exacta, que viene dada por: $x = (\theta-\sin(\theta))/2,\quad
y = 1-(1-\cos(\theta))/2;$ para $\theta\in[0,\pi]$.
\fej

\newpage

\clearpage
\setcounter{page}{1}
\setcounter{equation}{0}
\setcounter{ejer}{0}
\renewcommand{\theenumi}{\textbf{\arabic{enumi}}}

%----------------------------------------------------------
%Encabezado:

\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\Large{\sc Elementos de Cálculo Numérico / Cálculo Numérico}}}
 \vskip 0.2cm
 \centerline{\ttfamily Segundo Cuatrimestre 2016}
 \hrulefill

 \bigskip
 \centerline{\bf Trabajo Práctico N$^\circ$ 5: Ajuste de formas geométricas.}
 \bigskip
 
 En muchas aplicaciones es necesario resumir un conjunto de datos con una primitiva geométrica. Consideremos, por ejemplo, el caso de un círculo. Imaginemos la siguiente situación. Un robot provisto de un palpador puede obtener una serie de puntos sobre una circunferencia de eje cilíndrico, que puede estar afectado por desgaste o mala alineación. Se desea entonces calcular el radio y la posición del centro de su base a los datos medidos.
 
 Una idea simple es la siguiente. Un círculo de centro $(x_0,y_0)$ y radio $r$ está descripto por la ecuación $(x-x_0)^2+(y-y_0)^2=r^2$. Desarrollando: $x^2-2xx_0+x_0^2+y^2-2yy_0+y_0^2=r^2 \quad \Leftrightarrow \quad 2xx_0+2yy_0+(r^2-x_0^2-y_0^2)=x^2+y^2$.
 
 Si disponemos de $N$ datos $(x_i, y_i)$, $i=1,\dots, N$, esperaríamos encontrar valores $x_0, y_0$ y $r$ tales que se cumplen las $N$ ecuaciones 
 \[2x_ix_0+2y_iy_0+(r^2-x_0^2-y_0^2)=x_i^2+y_i^2.\]
 Este problema no tiene solución, pero podemos pensarlo en el sentido de cuadrados mínimos: Hallar $A, B, C$ tales que $x_iA+y_iB+C=x_i^2+y_i^2$, $i=1,\dots, N$ (en el sentido de que $\sum_{i=1}^N (x_iA+y_iB+C-(x_i^2+y_i^2))^2$ sea mínima). Este problema siempre tiene solución, y luego podemos interpretar: $A=2x_0$ ($\Rightarrow x_0=A/2$), $B=2y_0$ ($\Rightarrow y_0=B/2$), $C=r^2+x_0^2+y_0^2$ ($\Rightarrow r=\sqrt{C-\frac{A^2}{4}-\frac{B^2}{4}}$).
 
 \bej 
  Implementar esta idea en Octave, de modo que acepte una matriz de $N$ pares $(x_i,y_i)$ y devuelva los parámetros de la circunferencia.
 \fej
 
 \bej
  Chequear el programa anterior con datos generados según:
  \begin{enumerate}[label=\alph*)]
   \item círculos completos, sin ruido
   \item círculos completos, con ruido aleatorio
   \item un arco de circunferencia de amplitud $0<\alpha\leq 2\pi$. ¿Qué pasa para $\alpha$ chico (del orden de $\pi/4$, por ejemplo)? 
  \end{enumerate}

 \fej
 
 El abordaje anterior es lo que se llama ``ajuste algebraico''. Son los cuadrados de los residuos de la representación algebraica de la figura geométrica los que se minimizan, y no la distancia cuadrática a la curva. 
 
 \bej
   Para implementar un ``ajuste geométrico'', se desea minimizar 
   \[\varepsilon_g(x_0,y_0,r)=\sum_{i=1}^Nd_i^2(x_0,y_0,r),\]
   donde $d_i(x_0,y_0,r)$ es la distancia del punto $(x_i,y_i)$ al círculo de parámetros $x_0, y_0, r$.
   \begin{enumerate}[label=\alph*)]
    \item Para resolver el problema, y estimar los valores $x_0, y_0, r$, utilizaremos el método de Newton para minimizar $\varepsilon_g$. El método de Newton está dado por:
    $$v^{k+1} = v^k - H(v^k)^{-1}\nabla \varepsilon_g(v^k),$$ 
    donde $H(v)$ es la matriz hessiana de $\varepsilon_g$ evaluada en $v$. Se puede tomar como valor inicial, por ejemplo $v^0=[v_1^0, v_2^0,v_ 3^0]'$, donde $(v_1^0, v_2^0)$ es el punto promedio de los datos, y $v_3^0$ es el promedio de las distancias a ese centro. \\    
    Implementar un programa que reciba como input una función $f$ y un punto $z$ y calcule el vector gradiente de $f$ evaluado en $z$, y un programa que calcule el hessiano de $f$ evaluado en $z$. Ambos utilizando diferencias forward.
    \item Chequear el comportamiento del ajuste geométrico, con los valores de $x_0$, $y_0$ y $r$ por el método de Newton, sobre datos similares a los anteriores.
   \end{enumerate}
 \fej




 
 \newpage

\clearpage
\setcounter{page}{1}
\setcounter{equation}{0}
\setcounter{ejer}{0}
\renewcommand{\theenumi}{\textbf{\arabic{enumi}}}

%----------------------------------------------------------
%Encabezado:

\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\Large{\sc Elementos de Cálculo Numérico / Cálculo Numérico}}}
 \vskip 0.2cm
 \centerline{\ttfamily Segundo Cuatrimestre 2016}
 \hrulefill

 \bigskip
 \centerline{\bf Trabajo Práctico N$^\circ$ 6: Reducción del ancho de banda de matrices ralas.}
 \bigskip
 
 Una matriz rala es una donde hay muchos ceros. Ejemplos típicos provienen de la resolución de ecuaciones diferenciales en derivadas parciales en varias dimensiones, de la representación de grafos, etc.
 
 Cuando se debe resolver un sistema $Ax=b$, con $A$ una matriz rala, las técnicas de eliminación directa (eliminación gaussiana, descomposición LU, etc.) tienen el problema de que la matriz se llena de elementos no nulos. Y esto puede tornar el problema intratable.
 
 Una posible solución consiste en reordenar filas (lo que no cambia la solución) y/o reordenar las columnas (lo que sólo introduce permutaciones en las componentes de la solución). Este procedimiento mantiene la cantidad de ceros.
 
 Si se lograra, por ejemplo, que todos los elementos por debajo de una cierta ``envolvente'' $i_m(j)$, donde $j$ es el índice de la columna e $i_m(j)$ es el mayor índice de fila para el cual hay un elemento no nulo en la columna $j$, claramente no es necesario seguir triangulando hacia abajo dicha columna, y eso acota las posibilidades de llenado.
 
 \bej
  Fabricar matrices ralas con el comando \verb#sprand#. Calcular las descomposiciones LU de las matrices. Utilizar el comando \verb#spy# para ver la estructura de las matrices y de sus descomposiciones. Estudiar los comandos \verb#sparse# y \verb#full#.
 \fej
 
 \bej
  Pensar reordenamientos de filas que achiquen la envolvente, y verificar cómo se ve la matriz reordenada. Verificar también cómo quedan sus descomposiciones LU.
 \fej
 
 \bej
  Para Matlab: investigar y utilizar los comandos \verb#colamd#, \verb#amd#, \verb#symamd#, \verb#colperm#, \verb#symrcm#. Ver, utilizando el comando \verb#spy#, cómo quedan los ejemplos y calcular las descomposiciones LU de las matrices reordenadas. Comparar.
 \fej
 
 
 
 
 
\end{document}




