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

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


\begin{document}

%----------------------------------------------------------
%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 Práctica N$^\circ$ 2: Ecuaciones Diferenciales: Problemas de valores iniciales.}
 \bigskip


%----------------------------------------------------------
% 1
\bej
 Se considera la siguiente ecuación diferencial:
$$
\left\{
\begin{array}{rl}
y'(t) = & 2y(t) - 5\sen(t)\\
y(0) =  & 1
\end{array}
\right.
$$
cuya solución exacta es la función $y(t) = 2\sen(t) +
\cos(t)$. 
\begin{enumerate}[label=\alph*)]
\item Escribir la iteración del método de Euler para esta ecuación. 
\item Calcular el error de truncado local. 
\item ¿Qué paso $h$ debe elegirse para que el error al estimar $y\big(\frac{\pi}{2}\big)$ sea menor que $10^{-2}$?
\end{enumerate}
\fej


%----------------------------------------------------------
% 2
\bej Hacer el mapa de curvas integrales en la región $[0, 10]
\times [0, 10]$ de la ecuación diferencial
$$y'(t) = (y(t) - 5).(\cos^2(t) - 0.5),$$
graficando simultáneamente, para $k = 0, 1, \dots, 10$, la
solución que se obtiene utilizando el método de Euler con paso $h = 0.01$ y con condición inicial $$y(0) = k.$$
\fej

%----------------------------------------------------------
% 3
\bej Considerar el problema $y'=-2ty$, $y(0)=1$, con $t\ge 0$.
\begin{enumerate}[label=\alph*)]
\item Determinar una cota, en términos de $h$, para el error
cometido si se usa el método de Euler para calcular $y(1)$.

\item ¿Cómo deberí­a tomar $h$ si se desea que el error
cometido sea menor que $10^{-2}$?

\item Calcular la solución en $t=1$ usando el valor de $h$
obtenido en el item previo, y verificar las estimaciones previstas comparando con la solución exacta.
\end{enumerate}
\fej

%----------------------------------------------------------
% 4
\bej \label{ej:eulersimple}  Escribir un programa que implemente el método de Euler explí­cito para resolver ecuaciones de la forma 
$$\begin{array}{c} y' = f(t,y)\\ y(t_0) = y_0,\end{array}$$ tomando como parámetros la función $f$, los tiempos inicial y final $t_0$ y $t_f$, el paso $h$ y el dato inicial $y_0$; y arrojando como resultados el vector $t=(t_0,t_0+h,\dots,t_f)$ y la solución $y$. 
\fej

%----------------------------------------------------------
% 5
\bej Se quiere estimar, aplicando el método de Euler, el valor de $e$ como $y(1)$ donde $y(t)$ es solución de $y'=y,\ y(0)=1$. Hallar un paso $h$ de modo que el error cometido resulte menor que $10^{-3}$. Realizar el mismo trabajo para el método de Taylor de orden 2.
\fej


%----------------------------------------------------------
% 6
\bej Se quiere verificar numéricamente el orden de convergencia de los métodos de Euler y Taylor de orden 2. Para ello: resolver numéricamente la ecuación del ejercicio anterior, en el intervalo $[0,1]$ con ambos métodos, tomando $h=0.1$, $h=0.0625$, $h=0.05$, $h=0.025$ y $h=0.01$. Para cada $h$ calcular el error que se comete al aproximar $y(1)$: $e_h = |y(1)-y_N|$. Graficar $\log(e_h)$ en función de $\log(h)$. ¿Qué se espera ver? ¿El resultado es consistente con el esperado?
\fej


%----------------------------------------------------------
% 7
\bej Considerar el problema $\left\{\begin{array}{ll}
y'=\lambda y \\
y(0)=y_0 \end{array}\right.$.
 
\begin{enumerate}[label=\alph*)]
\item Verificar que el método de Euler con paso $h$ genera la
sucesión:
\[ y_i=(1+\lambda h)^iy_0 \ \ i=0,1,\ldots.  \]
\item  Para $\lambda <0$, determinar para qué valores de $h$
ocurre que $y_i \rightarrow 0$ cuando $i \to \infty$. Comparar la con la solución exacta. 
\item Resolver usando el programa del Ejercicio \ref{ej:eulersimple} para distintos valores de $\lambda$ ($\lambda=1,10,50,100$) y comparar con la solución exacta. ¿Qué sucede?
\item Implementar un programa que resuelva el problema usando el método de Euler implí­cito y repetir los experimentos del item anterior.
\end{enumerate}
\fej

%----------------------------------------------------------
% 8
\bej Considerar la ecuación: 
$$\left\{\begin{array}{rcl} y'(t) & = &\dfrac{e^{-y}}{t}\\ y(1) &= &0.\end{array}\right.$$
\begin{enumerate}[label=\alph*)]
\item Probar que $0\le y(t)\le t$ para $t\ge 1$. 
\item Escribir la iteración dada para esta ecuación por el método de Euler. Probar que la solución numérica resultará creciente. 
\item Calcular el error de truncado del método de Euler aplicado a la ecuación.
\item Dar un valor de paso $h$ que garantice que el error de la estimación numérica de $y(2)$ sea menor que $10^{-3}$.
\end{enumerate}
\fej


%----------------------------------------------------------
% 9
\bej Probar que una ecuación de orden $n$:
$$y^{(n)} = f\big(t,y,y',\dots,y^{(n-1)}\big),$$ 
se puede escribir como un sistema de $n$ ecuaciones de primer orden. Mostrar que un problema de valores iniciales para la primera se transforma en un problema de valores iniciales para el sistema.
\fej

%----------------------------------------------------------
% 10
\bej
Modificar el programa del Ejercicio \ref{ej:eulersimple} para que acepte ecuaciones vectoriales: la solución $y$ deberá ser una matriz de $m\times n$, donde $m$ es el número de pasos temporales y $n$ la cantidad de variables del problema. De este modo, la fila $i$ de $y$ corresponderá al valor de la solución en todas sus variables a tiempo $t_i$ 
\fej

%----------------------------------------------------------
% 11
\bej\label{ej:particula}
La trayectoria de una partí­cula que se mueve en el plano está dada por las curva $(y_1(t), y_2(t))$, donde las funciones $y_1, y_2$ son la solución del siguiente sistema de ecuaciones diferenciales:
\[
\begin{array}{lcl}
y_1'(t)&=&-y_2(t)\\
y_2'(t)&=&y_1(t)-y_2(t)
\end{array}.
\]
Resolver este sistema en el intervalo $[0,20]$ con el método de Euler utilizando paso $h = 0.05$ y graficar la trayectoria de la partí­cula, sabiendo que en tiempo $t= 0$ se encontraba en el punto $(1,-1)$. 
\fej

%----------------------------------------------------------
% 12
\bej
Considerar el método de Euler modificado:
$$y_{i+1} = y_i + h~f\left(t_i+\frac{h}{2},y_i+\frac{h}{2}f(t_i,y_i)\right).$$
Probar que el error de truncado es $O(h^2)$. 
\fej


%----------------------------------------------------------
% 13
\bej\label{ej:RK4}
Implementar un programa que resuelva sistemas de la forma:
$$y' = f(t,y),$$ utilizando el método de Runge Kutta de orden 4 dado por:
$$y_{i+1} = y_i + \frac{h}{6}\Big[K_1 +2K_2 +2K_3 +K_4\Big],$$
donde:
$$
\begin{array}{c}
K_1 = f(t_i,y_i),\\
K_2 = f\big(t_i+\frac{h}{2},y_i+\frac{h}{2}K_1\big),\\
K_3 = f\big(t_i+\frac{h}{2},y_i+\frac{h}{2}K_2\big),\\
K_4 = f\big(t_i+h,y_i+hK_3\big).
\end{array}
$$
Utilizar este método para resolver nuevamente el Ejercicio \ref{ej:particula}. Comparar la solución con la obtenida con el método de Euler.
\fej



%----------------------------------------------------------
% 14
\bej\label{ej:tirooblicuo}{\bf Tiro oblicuo:}
Un proyectil de masa $m$ se arroja desde un punto del plano $(x_0,y_0)$, con una velocidad inicial dada por el vector $(v^x_0,v^y_0)$. La trayectoria del proyectil se rige por las ecuaciones dadas por la segunda ley de Newton:
$$\begin{array}{rcl}
m\ddot{x} &=& -\gamma v^x\\ 
m\ddot{y} &=& -mg -\gamma v^y,
\end{array}$$
donde $g$ es la aceleración gravitatoria $g=9.81\frac{\textrm{m}}{\textrm{s}^2}$, y $\gamma$ es una constante de rozamiento con el medio en que se realiza el lanzamiento. 
Formular el problema en forma de sistema de orden uno. 

Tomando $m=10 \textrm{Kg}$ y $\gamma=0.2 \frac{\textrm{Kg}}{\textrm{s}}$, y suponiendo que el proyectil se lanza desde $30$ metros de altura con una velocidad inicial horizontal de $40\frac{\textrm{m}}{\textrm{s}}$, ¿qué distancia recorre antes de tocar el piso?

Hacer un programa que permita responder esta pregunta, utilizando el método de Euler modificado para resolver el sistema. 
\fej

%----------------------------------------------------------
% 15
\bej\label{pendulo}{\bf Péndulo:} Se considera el problema del péndulo:
$$\left\{\begin{array}{ll}
\ddot{\theta}(t)=-A\sen( \theta(t)) \mbox{ en } [0,T]\\
\dot{\theta}(0) = v_0 \\
\theta(0)=\theta_0 
\end{array}\right.$$
donde $\theta$ representa el ángulo que forma la vara del péndulo con la vertical.

\begin{enumerate}[label=\alph*)]
\item Formular el problema como un sistema de ecuaciones de orden uno.
\item Utilizar el método de Euler modificado, con paso $h=0.05$ para obtener una aproximación de la soluci\'on y graficarla.  
\item Graficar la solución que se obtiene al utilizar método de Runge Kutta del Ejercicio \ref{ej:RK4}.
\end{enumerate}

Pueden considerarse, a modo de ejemplo, los valores $A=7$, $T=10$, $\theta_0=\frac{\pi}{4}$, $v_o=0$.
\fej


%----------------------------------------------------------
% 16
\bej
Se considera el oscilador de Van der Pol:
\[\dfrac {d^{2}x}{dt^{2}}-\mu (1-x^{2})\dfrac{dx}{dt}+x=0,\]
donde $x$ es la posición, en función del tiempo $t$, y $\mu$ es un parámetro que indica la no linealidad y el amortiguamiento.
\begin{enumerate}[label=\alph*)]
\item Resolver numéricamente el sistema para $\mu=0.1$ y valor inicial $x=1$, $\dfrac{dx}{dt}=1$  utilizando \verb+ode45+. ¿Qué se observa?
\item Resolver numéricamente el sistema para $\mu=0.1$ y $\mu=1.4$ y valores iniciales $x=0$, $\dfrac{dx}{dt}=0.1$ y $x=3$, $\dfrac{dx}{dt}=3$  utilizando \verb+ode23s+. Si $y=\frac{dx}{dt}$, graficar en todos los casos $x$ e $y$ en función del tiempo $t$, e $y$ en función de $x$ (el diagrama de fases).
\end{enumerate}
(Observación: En \verb+Octave+ los comandos \verb+ode45s+ y \verb+ode23s+ se encuentran en el paquete \verb+odepkg+.)
\fej

%----------------------------------------------------------
% 17
\bej
Considerar el problema de la propagación de una llama:
$$\left\{\begin{array}{ll}
\dot{y}=y^2-y^3 \mbox{ para } t\in[0,2/\delta]\\
y(0) = \delta. \\
\end{array}\right.$$
Al prender un fósforo, la bola de fuego crece rápidamente hasta que alcanza un tama\~no crítico. Luego, permanece en ese tama\~no porque la cantidad de oxígeno que se consume por la combustión en el interior de la bola se compensa con la cantidad disponible en la superficie. El radio de la bola de fuego se representa con $y(t)$ y los términos $y^2$ e $y^3$ provienen de la superficie y el volumen. El parámetro crítico es el radio inicial $\delta$.
\begin{enumerate}[label=\alph*)]
\item Resolver numéricamente el sistema para $\delta=0.0001$ y un error relativo de $10^{-4}$,  utilizando \verb+ode45+. ¿Qué se observa? (Usar, por ejemplo, los comandos \verb+tic+ y \verb+toc+.)
\item Repetir el ítem anterior con el comando \verb+ode23s+.
\end{enumerate}

\fej
\end{document}