\documentclass[10pt]{article}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage{graphicx}
\usepackage{amsmath, amssymb}
\usepackage{verbatim}
\usepackage{hyperref}
\usepackage{showlabels}

\usepackage[spanish]{babel}


\begin{document}
\begin{titlepage}

%\raggedright
\centering{\includegraphics[width=0.2\textwidth]{URJC-Logo-PNG-1.eps}\par}
\vspace{1cm} %\raggedleft
\begin{center} {\bfseries\Large C\'ODIGOS EN OCTAVE/MATLAB
\par} \vspace{0.5cm}
{\bfseries\Large SEMINARIOS Y PR\'ACTICAS\par} \vspace{0.5cm}
{\bfseries\Large M\'ETODOS MATEM\'ATICOS APLICADOS\par}
\vspace{0.5cm} {\bfseries\Large A LA INGENIER\'IA\par}
%\vspace{0.5cm} {\bfseries\Large \textit{EJERCICIOS Y PROBLEMAS
%RESUELTOS}
%\par}
 \vspace{1cm}

{\scshape Asignaturas:\par}
 {\scshape \large  M\'etodos Matem\'aticos
aplicados a la Ing. de Materiales}\\ {\scshape \large  M\'etodos
Matem\'aticos aplicados a la Ing. de la Energ\'{\i}a}\\ {\scshape
\large M\'etodos Num\'ericos (m\'odulo I) en el M\'aster en Ing.
Industrial
\par} \vspace{0.5cm}
%{\scshape\Large PARA M\'ASTER Y GRADOS EN INGENIER\'IA \par}
\vspace{1cm}

%{\itshape\Large M\'etodos Matem\'aticos aplicados a la
%Ingenier\'{\i}a  \par}

%\vfill {\Large Autores: \par}
{\Large A. I. Mu\~noz Montalvo, A. Nolla y E. Schiavi\par}

{\Large Septiembre 2022
\par}\vspace{1cm}

{\large{\copyright 2022. Autores: A. I. Mu\~noz
Montalvo, A. Nolla, E. Schiavi.\\
Algunos derechos reservados.\\
Este documento se distribuye bajo la licencia internacional
Creative Commons Attribution-ShareAlike 4.0 International
License.\\ Disponible en:
http://creativecommons.org/licenses/by-sa/4.0/}}\vspace{0.5cm}

{\large{Publicado en: https://burjcdigital.urjc.es}}
\end{center}

\end{titlepage}

\tableofcontents

\newpage

\section{C\'ODIGOS EN OCTAVE/MATLAB }

 \vspace{1cm}

{\bf{Observaciones:}} Todo aquello escrito en cursiva y que no
vaya despu\'es del s\'{\i}mbolo prompt del sistema, \texttt{$>$},
debe ser interpretado como un comentario o explicaci\'on, y no
como una l\'{\i}nea de comandos.\vspace{0.5cm}

 La mayor parte de los c\'odigos son adaptaciones de las funciones del libro
{\em ``C\'alculo cient\'ifico con MATLAB y Octave''} de A.
Quarteroni, F. Saleri, que se pueden obtener
https://mox.polimi.it/qs/.

\vspace{1.5cm}



\subsection{ C\'odigos para resolver ecuaciones no
lineales }\vspace{1cm}

\subsubsection{M\'etodo de bisecci\'on}

\texttt{$>$
function[sol,itera]=metbiseccion(fecu,a,b,errorper,maxitera)}

  \textit{Este c\'odigo encuentra una aproximaci\'on de una ra\'{\i}z en el intervalo
  [a,b], de la ecuaci\'on fecu=0. Los argumentos son:
  \begin{itemize}
  \item errorper: el error permitido.
  \item maxitera: el n\'umero m\'aximo de iteraciones que se permiten
  realizar.
   \item sol: la soluci\'on num\'erica aproximada obtenida.
  \item itera: el n\'umero de iteraciones que ha realizado para obtener
  sol.
  \end{itemize}}

\texttt{$>$ x = [a, (a+b)*0.5, b]; fx = fecu(x);}

\textit{Comprobamos que el intervalo elegido satisface el teorema
de Bolzano, fecu(a)*fecu(b)$<$0, si no lo verifica, aparece un
mensaje de error. Tambi\'en puede ocurrir que a o b sean ya
ra\'{\i}ces.}

\texttt{$>$ if fx(1)*fx(3) > 0}

  \texttt{$>$ error([' El signo de la funci\'on en los
  extremos',...}

  \texttt{ 'del intervalo (a,b) tiene que ser distinto']);}

\texttt{$>$ elseif fx(1) == 0}

  \texttt{$>$  sol = a; itera = 0;}

  \texttt{$>$ return}

\texttt{$>$ elseif fx(3) == 0}

    \texttt{$>$ sol = b; itera = 0;}

    \texttt{$>$ return}


\texttt{$>$ end}

\texttt{$>$ itera = 0; I = (b - a)*0.5;}

\texttt{$>$ while I >= errorper \&\& itera < maxitera}

 \texttt{$>$  itera = itera + 1;}
   \texttt{$>$ if fx(1)*fx(2) <  0}

     \texttt{$>$ x(3) = x(2);}

      \texttt{$>$ x(2) = x(1)+(x(3)-x(1))*0.5;}

      \texttt{$>$ fx = fecu(x);}

      \texttt{$>$ I = (x(3)-x(1))*0.5;}

   \texttt{$>$ elseif fx(2)*fx(3) < 0}

      \texttt{$>$ x(1) = x(2);}

      \texttt{$>$ x(2) = x(1)+(x(3)-x(1))*0.5;}

     \texttt{$>$ fx = fecu(x);}

     \texttt{$>$ I = (x(3)-x(1))*0.5;}

   \texttt{$>$ else}

     \texttt{$>$  x(2) = x(find(fx==0)); I = 0;}

   \texttt{$>$ end}

\texttt{$>$ end}

  \texttt{$>$ if(itera==maxitera \&\& I > errorper)}

 \texttt{$>$ fprintf(['El m\'etodo no converge,'...}

  \texttt{$>$ 'se ha llegado al n\'umero m\'aximo de iteraciones
  ',...}

 \texttt{$>$  'sin estar por debajo del error m\'aximo
 permitido']);}

\texttt{$>$ end}

\texttt{$>$ sol = x(2);}

\vspace{1cm}

\subsection{M\'etodo del punto fijo}

\texttt{$>$
function[sol,itera]=metpuntofijo(g,x0,errorper,maxitera)}

\textit{Este m\'etodo iterativo encuentra una aproximaci\'on del
punto fijo de la funci\'on g. N\'otese que f(x)=0 $<=>$ g(x)=x.
Los argumentos son:
\begin{itemize}
\item g: funci\'on que define el esquema num\'erico.
\item x0: semilla para iniciar el esquema iterativo.
\item sol: soluci\'on num\'erica obtenida.
\item  errorper: el error m\'aximo permitido. El error se mide mediante el valor absoluto de la diferencia de dos valores
obtenidos en iteraciones consecutivas.
\item maxitera: el n\'umero m\'aximo de iteraciones permitido.
\item itera: el n\'umero de
iteraciones realizado para obtener sol, cumpliendo con el error
m\'aximo permitido.
\end{itemize}
}

\texttt{$>$ x = x0; diferencia= errorper +1; itera =0;}

\texttt{$>$ while itera <= maxitera \&\& diferencia >= errorper}

\texttt{$>$ gx = feval(g,x);} %x(i+1)=g(x(i)), xnueva=g(xanterior)

\texttt{$>$ diferenciait=x-gx;}%x(i+1)-x(i)

 \texttt{$>$ diferencia=abs(diferenciait);}

 \texttt{$>$ xnueva=gx;}

 \texttt{$>$ x=xnueva;}

\texttt{$>$ itera=itera+1;}

 \texttt{$>$ end}

\texttt{$>$ if itera >= maxitera}

\texttt{$>$ fprintf('Se ha llegado al n\'umero m\'aximo de
iteraciones',...}

  \texttt{$>$        'sin estar por debajo del m\'aximo error permitido');}

 \texttt{$>$ end}

\texttt{$>$ sol = x;}

\texttt{$>$ return}

\texttt{$>$ endfunction}

\vspace{1cm}


\subsubsection{M\'etodo Aitken}

\texttt{$>$ function
[sol,itera]=metodoaitken(g,x0,errorper,maxitera)}

\textit{Este c\'odigo utiliza el m\'etodo de Aitken para obtener
una aproximaci\'on num\'erica al punto fijo de una funci\'on g,
utilizando como semilla x0. Los argumentos son:
\begin{itemize}
\item g: funci\'on que define el esquema num\'erico.
\item x0: semilla para iniciar el esquema iterativo.
\item sol: soluci\'on num\'erica obtenida.
\item  errorper: el error m\'aximo permitido. El error se mide mediante el valor absoluto de la diferencia de dos valores
obtenidos en iteraciones consecutivas.
\item maxitera: el n\'umero m\'aximo de iteraciones permitido.
\item itera: el n\'umero deiteraciones realizado para obtener sol, cumpliendo con el error
m\'aximo permitido.
\end{itemize}}

\texttt{$>$ x = x0; difer = errorper + 1; itera = 0;}

\texttt{$>$ while itera < maxiter \&\& difer >= errorper}

\texttt{$>$ gx = g(x); ggx = g(gx);xnew
=(x*ggx-gx\^\,2)/(ggx-2*gx+x);}

\texttt{$>$ difer = abs(x-xnew);}

\texttt{$>$ x = xnew;}

 \texttt{$>$   itera = itera  + 1;}

\texttt{$>$ end}

\texttt{$>$ if}

\texttt{$>$ (itera==maxitera \&\& difer>errorper)}

 \texttt{$>$   fprintf(['El m\'etodo no converge en el n\'umero',...}

 \texttt{$>$            'm\'aximo de iteraciones fijado']);}

\texttt{$>$ end}

\texttt{$>$ sol = x;}

 \texttt{$>$ return}



\subsubsection{M\'etodo de Newton Raphson para una ecuaci\'on de una inc\'ognita}

\texttt{$>$ function[sol,itera]=metnewton1ec(fecu,dfecu,x0,errorper,maxitera)}\\

\textit{Este c\'odigo obtiene una soluci\'on aproximada de una
ra\'{\i}z de la ecuaci\'on fecu=0, con el m\'etodo de
Newton-Raphson, tomando como semilla x0. Los argumentos son:
\begin{itemize}
\item x0: semilla para arracar el esquema iterativo.
\item  fecu: es la funci\'on que define la ecuaci\'on de la cual
queremos obtener una ra\'{\i}z.
\item dfecu: derivada primera de fecu.
\item sol: la soluci\'on num\'erica obtenida.
\item errorper: el m\'aximo error permitido en la
aproximaci\'on.
\item  maxitera: el n\'umero m\'aximo de iteraciones permitidas.
\item itera: el n\'umero de iteracioes realizadas por el esquema
para   alcanzar un error por debajo del error cometido. El error
se mide con  el valor absoluto de fecu(x)/dfecu(x).
\end{itemize}}

\texttt{$>$ x = x0;}

\texttt{$>$  fx = fecu(x); dfx = dfecu(x); itera = 0;medidaerror =
errorper+1;}

\texttt{$>$  while medidaerror >= errorper \&\& itera < maxitera}

\texttt{$>$    itera = itera + 1;}

 \texttt{$>$   medidaerror = - fx/dfx;}

\texttt{$>$    x = x + medidaerror;}

\texttt{$>$    medidaerror= abs(medidaerror);}

   \texttt{$>$ fx = fecu(x);}

 \texttt{$>$   dfx = dfecu(x);}

\texttt{$>$ end}

\texttt{$>$  if (itera==maxitera \&\& medidaerror > errorper)}

\texttt{$>$   fprintf(['El m\'etodo no converge por alcanzar el
n\'umero m\'aximo',...}

 \texttt{$>$  'de iteraciones permitidas',...}

 \texttt{$>$  'sin llegar a un error por debajo del error m\'aximo
 permitido']);}

\texttt{$>$ end}

\texttt{$>$ sol = x;}

\texttt{$>$ return}

\vspace{1cm}

\subsubsection{M\'etodo de la secante}

\texttt{$>$ function
[sol,itera]=metsecante(fecu,a,b,errorper,maxitera)}

\textit{ Este c\'odigo encuentra una aproximaci\'on de una
ra\'{\i}z de la ecuaci\'on dada por fecu=0, utilizando el m\'etodo
de la secante. Los argumentos son:
\begin{itemize}
\item a y b: semillas.
\item  sol: la soluci\'on num\'erica obtenida.
\item errorper: el error m\'aximo permitido. El error se mide sobre el valor absoluto de la diferencia de dos valores
obtenidos en iteraciones consecutivas.
\item  maxitera: el n\'umero m\'aximo de iteraciones permitido.
\item  itera: el n\'umero de iteraciones realizadas para obtener
sol, cumpliendo con el error m\'aximo permitido.
\end{itemize}}

\texttt{$>$ fx0=feval(fecu,a); fx1=feval(fecu,b); x0=a; x1=b;}

\texttt{$>$  r=fx1*(x1-x0)/(fx1-fx0); x2=x1-r; c=1;}

\texttt{$>$ while (abs(r)>errorper \&\& c < maxitera)}

 \texttt{$>$   x0=x1;    x1=x2;}

\texttt{$>$    fx0=feval(fecu,x0);    fx1=feval(fecu,x1);}

\texttt{$>$    r=fx1*(x1-x0)/(fx1-fx0);}

\texttt{$>$    x2=x1-r;}

\texttt{$>$    c=c+1;}

\texttt{$>$ end}

\texttt{$>$ if c >= maxitera}

\texttt{$>$ fprintf('No converge en el m\'aximo',...}

\texttt{$>$          'n\'umero de iteraciones');}

\texttt{$>$ end}

\texttt{$>$ sol = x2;}

\texttt{$>$ itera = c;}

\texttt{$>$ end}



\vspace{1cm}

\subsubsection{M\'etodo de regulafalsi}

\texttt{$>$ function[sol,itera]=metregulafalsi(fecu,a,b,errorper,maxitera)}\\

\textit{Este c\'odigo encuentra una aproximaci\'on de una
ra\'{\i}z de la ecuaci\'on dada por fecu=0, utilizando el m\'etodo
de regulafalsi. Los argumentos son:
\begin{itemize}
\item fecu: funci\'on que define la ecuaci\'on.
\item a y b semillas.
\item sol: la soluci\'on num\'erica obtenida tomando como
semillas a y b.
\item errorper: el error m\'aximo permitido. El error se mide sobre el valor absoluto de la diferencia de dos
valores obtenidos en iteraciones consecutivas.
\item maxitera: el n\'umero m\'aximo de iteraciones permitido.
\item itera: el n\'umero de iteraciones realizadas para obtener sol,
cumsubpliendo con el error m\'aximo permitido.
\end{itemize}}

\texttt{$>$ fx0=feval(fecu,a); fx1=feval(fecu,b); x0=a; x1=b;}

\texttt{$>$ r=fx1*(x1-x0)/(fx1-fx0); x2=x1-r; c=1;}

\texttt{$>$ fx2=feval(fecu,x2);}

\texttt{$>$ while (abs(r)>errorper \&\& c < maxitera)}

\texttt{$>$    if fx2*fx1<0}

\texttt{$>$  x0=x2;   x1=x1;}

    \texttt{$>$ else}

 \texttt{$>$       x0=x0;       x1=x2;}

    \texttt{$>$ end}

 \texttt{$>$   fx0=feval(fecu,x0);    fx1=feval(fecu,x1);}

   \texttt{$>$ r=fx1*(x1-x0)/(fx1-fx0);}

 \texttt{$>$   x2=x1-r;    fx2=feval(fecu,x2);}

    \texttt{$>$c=c+1;}

\texttt{$>$ end}

\texttt{$>$ if c >= maxitera}

\texttt{$>$ fprintf(' No converge en el m\'aximo',...}

\texttt{$>$          'n\'umero de iteraciones ');}

\texttt{$>$ end}

\texttt{$>$ sol = x2; itera = c;}

\texttt{$>$ end}





\vspace{1cm}
\subsubsection{M\'etodo de Newton Raphson para un sistema de ecuaciones}

\texttt{$>$ function [vectorsol,itera] =
metnewtonsistema(fecusistema,...}\\

\texttt{$>$jacobiana,vectorx0,errorper,maxitera)}\\

\textit{Este c\'odigo obtiene una soluci\'on aproximada del
sistema no lineal dado por el vector de funciones fecusistema, F,
igualado a cero, F=0, partiendo de la semilla vectorx0. Los
argumentos son:
\begin{itemize}
\item fecusistema: funciones que definen el sistema de
ecuaciones.
\item  jacobiana: la funci\'on matriz jacobiana J.
\item fecusistema y jacobiana deben ser definidas como funciones en
scripts .m
\item vectorx0: vector semilla.
\item errorper: el m\'aximo error permitido, que se mide tomando el
m\'odulo del vector -[J(vectorx)]${}^{-1}$F(vectorx).
\item vectorsol: es la aproximaci\'on del vector ra\'{\i}z que hemos
obtenido.
\end{itemize}}

\texttt{$>$ itera = 0; medidaerror = errorper + 1; vectorx=
vectorx0;}

\texttt{$>$ while medidaerror >= errorper \&\& itera < maxitera}

\texttt{$>$    J = fecusistema(vectorx(1),vectorx(2));}

\texttt{$>$    F = jacobiana(vectorx(1),vectorx(2));}

\texttt{$>$    incremento = - inv(J)*F;}

\texttt{$>$ vectorx = vectorx + incremento;}

\texttt{$>$    medidaerror = norm(incremento);}

\texttt{$>$ itera = itera + 1;}

\texttt{$>$ end}

\texttt{$>$ vectorsol=vectorx;}

\texttt{$>$ if (itera==maxitera \&\& medidaerror> errorper)}

\texttt{$>$ fprintf(['El m\'etodo no converge en el n\'umero
m\'aximo',...}

\texttt{$>$   'de iteraciones',...}

\texttt{$>$       'por no alcanzarse un error inferior al
permitido'],F);}

\texttt{$>$ end}

\texttt{$>$ return}

\vspace{0.5cm}

 \textbf{Funci\'on fecusistema.m. Ejemplo de implementaci\'on}\\

\texttt{$>$ function F=fecusistema(x,y)}

 \texttt{$>$ F(1,1)=(x+y)\^\,2-2;}

 \texttt{$>$ F(2,1)=x\^\,2+1-y;}

 \texttt{$>$ return}


\vspace{0.5cm}

\textbf{Funci\'on jacobiana.m. Ejemplo de
implementaci\'on}\\

\texttt{$>$ function J=jacobiana(x,y)}

\texttt{$>$ J(1,1)=2*(x+y);}

\texttt{$>$ J(1,2)=2*(x+y);}

\texttt{$>$ J(2,1)=2*x;}

\texttt{$>$ J(2,2)=-1;}

\texttt{$>$ return}


\newpage


\subsection{ C\'odigos para resolver problemas de
valor inicial}

\subsubsection{M\'etodo de Euler expl\'{\i}cito}

\texttt{$>$ function[solt,soly]=eulerexplicito(f,intiempo,valorini,npasos)}\\

\textit{Este c\'odigo resuelve el problema de valor inicial dado
por la ecuaci\'on diferencial y' = f(t,y) y el valor inicial,
y(t0), denotado por valorini, en el intervalo temporal [t0,T)
(definido en intiempo), con el m\'etodo de Euler expl\'{\i}cito.
El resto de argumentos son:\begin{itemize}
\item npasos: el n\'umero de pasos que hay que dar, tomando un
tama\~no de discretizaci\'on constante h, para llegar a T, es
decir, (T-t0)/h.
\item solt: el vector que contiene los
nodos temporales, t0,t1, ...tn=T.
\item soly: el vector que
 contiene los valores num\'ericos aproximados de
y(t), obtenidos para los distintos nodos temporales.
\end{itemize}}

\texttt{$>$   h=(intiempo(2)-intiempo(1))/npasos;}
%inicializamos el vector que albergar� la soluci�n

\texttt{$>$   solt=ones(npasos+1,1);}

\texttt{$>$   soly=ones(npasos+1,1);}

\texttt{$>$   solt(1)=intiempo(1);}%t0

\texttt{$>$   soly(1)=valorini;}%y(t0)

\texttt{$>$   for i=2:npasos+1}

 \texttt{$>$   solt(i)=solt(i-1)+h;}

\texttt{$>$ soly(i)=soly(i-1)+h*f(solt(i-1),soly(i-1));}

\texttt{$>$   end}

\texttt{$>$   return}

\subsubsection{
M\'etodo de Euler impl\'{\i}cito}

\texttt{$>$ function
[solt,soly]=eulerimplicito(f,intiempo,valorini,npasos)}\\

\textit{Este c\'odigo resuelve un problema de valor inicial
y'=f(t,y), y(t0)=y0, uti\-lizando el m\'etodo de Euler
impl\'{\i}cito. Los argumentos son:
\begin{itemize}
\item f: la funci\'on que define y'.
\item intiempo: el intervalo de tiempo donde se
quiere resolver el problema.
\item valorini: el dato inicial.
\item npasos: el n\'umero de pasos que hay que dar para llegar al
tiempo final y depende del tama\~no de discretizaci\'on h.
\item solt y soly: los vectores soluci\'on, e indican los tiempos
intermedios considerados y los correspondientes valoresaproximados de y.
\item Obs: utilizamos el comando fsolve para
resolver la posible ecuaci\'on no lineal en y que pueda aparecer
al aplicar el esquema impl\'{i}cito.
\end{itemize}}

\texttt{$>$   h=(intiempo(2)-intiempo(1))/npasos;}

\texttt{$>$ solt=(intiempo(1):h:intiempo(2));}

\texttt{$>$  soly=ones(npasos+1:1); soly(1)=valorini;} % always create a vector column

\texttt{$>$ global $glob_h$ $glob_t$ $glob_y$ $glob_f$;}

\texttt{$>$ $glob_h=h$;$glob_y$=valorini; $glob_f$=f;}

\texttt{$>$   for i=2:npasos+1}

\texttt{$>$       $glob_t$=solt(i);}

\texttt{$>$    w = fsolve(@(w) beulerfun(w),$glob_y$);}

\texttt{$>$   soly(i)=w;}

\texttt{$>$     $glob_y$=w;}

\texttt{$>$    end}

\texttt{$>$   clear $glob_h$ $glob_t$ $glob_y$ $glob_f$;}

\texttt{$>$   end}

\texttt{$>$    function[z]=beulerfun(w)}

\texttt{$>$     global $glob_h$ $glob_t$ $glob_y$ $glob_f$;}

  \texttt{$>$   z=w-$glob_y$-$glob_h$*feval($glob_f$,$glob_t$,w);}

\texttt{$>$   end}






\subsubsection{M\'etodo de Crank-Nicolson}
\texttt{$>$
function[solt,soly]=cranknicolson(f,intiempo,valorini,npasos)}\\

\textit{Este c\'odigo resuelve un problema de valor inicial
y'=f(t,y), y(t0)=y0, utilizando el m\'etodo de Crank Nicolson. Los
argumentos son:
\begin{itemize}
\item f: funci\'on que define y'.
\item intiempo: el intervalo de tiempo donde se quiere resolver el problema.
\item valorini: el dato inicial.
\item npasos: el n\'umero de pasos que hay que dar para llegar al
tiempo final y depende del tama\~no de discretizaci\'on h.
\item solt y soly: los vectores soluci\'on, e indican los tiempos
intermedios considerados y los correspondientes valores
aproximados de y.
\item Obs: utilizamos el comando fsolve para resolver la posible ecuaci\'on no lineal
en y que pueda aparecer al aplicar el esquema impl\'{\i}cito.
\end{itemize}}

\texttt{$>$   h=(intiempo(2)-intiempo(1))/npasos;}

\texttt{$>$   solt=(intiempo(1):h:intiempo(2));}

\texttt{$>$   soly=ones(npasos+1:1); soly(1)=valorini;}

\texttt{$>$   global $glob_h$ $glob_t$ $glob_tt$ $glob_y$
$glob_f$;}\\

\texttt{$>$  $glob_h=h$; $glob_y$=valorini; $glob_f$=f;}\\

\texttt{$>$   for i=2:npasos+1}

    \texttt{$>$   $glob_t$=solt(i);}

    \texttt{$>$   $glob_tt$=solt(i-1);}

\texttt{$>$     w = fsolve(@(w) crank(w),$glob_y$);}

  \texttt{$>$   soly(i)=w;}

  \texttt{$>$   $glob_y$=w;}

\texttt{$>$   end}

\texttt{$>$   clear $glob_h$ $glob_t$ $glob_tt$ $glob_y$
$glob_f$;}

\texttt{$>$    end}

\texttt{$>$    function [z]=crank(w)}

\texttt{$>$   global $glob_h$ $glob_t$ $glob_tt$ $glob_y$
$glob_f$;}

\texttt{$>$ z=w-$glob_y$-0.5.*$glob_h$*feval($glob_f$,$glob_t$,w)-
...}

\texttt{$>$ 0.5.*$glob_h$*feval($glob_f$,$glob_tt$,$glob_y$);}

\texttt{$>$    end}





\subsubsection{M\'etodo de Heun}

\texttt{$>$   function
[solt,soly]=heun(f,intiempo,valorini,npasos)}\\

\textit{Este c\'odigo resuelve un problema de valor inicial
y'=f(t,y), y(t0)=y0, utilizando el m\'etodo de Heun visto como un
caso de m\'etodo Runge Kutta. Los argumentos son:
\begin{itemize}
\item f: la funci\'on que define y'.
\item  intiempo: el intervalo de tiempo donde se quiere resolver el problema.\item  valorini: el dato inicial.
\item npasos: el n\'umero de pasos que hay que dar para llegar al
tiempo final y depende del tama\~no de discretizaci\'on h.
\item solt y soly: son los vectores soluci\'on, e indican los
tiempos intermedios considerados y los correspondientes valores
aproximados de y.
\end{itemize}}

\noindent\texttt{$>$ h=(intiempo(2)-intiempo(1))/npasos;}\\
\texttt{$>$ soly=ones(npasos+1:1);}\\
\texttt{$>$ solt=linspace(intiempo(1),intiempo(2),npasos+1);}\\
\texttt{$>$ soly(1)=valorini;}\\
\texttt{$>$ for i=2:npasos+1}\\
 \texttt{$>$ tn1=solt(i-1);}\\
\texttt{$>$  tn2=solt(i);}\\
\texttt{$>$  yn1=soly(i-1);}\\
\texttt{$>$  yn2=soly(i-1)+h*f(tn1,yn1);}\\
 \texttt{$>$ soly(i) =soly(i-1)+
 0.5*h*(f(tn1,yn1)+f(tn2,yn2));}\\
\texttt{$>$ end}\\



\subsubsection{M\'etodo de Simpson}

\texttt{$>$   function
[solt,soly]=rungekuttao3(f,intiempo,valorini,npasos)}\\

\textit{Este c\'odigo resuelve un problema de valor inicial
y'=f(t,y), y(t0)=y0, utilizando el m\'etodo Runge Kutta de orden,
que denominamos Simpson por la f\'ormula de integraci\'on
num\'erica que usa. Los argumentos son:
\begin{itemize}
\item f: la funci\'on que define y'.
\item  intiempo: el intervalo de tiempo donde se
quiere resolver el problema.
\item valorini: el dato inicial.
\item npasos: el n\'umero de pasos que hay que dar para llegar al
tiempo final y depende del tama\~no de discretizaci\'on h.
\item solt y soly: los vectores soluci\'on, e indican los
tiempos intermedios considerados y los correspondientes valores
aproximados de y.
\end{itemize}}

\noindent\texttt{$>$ solt=linspace(intiempo(1),intiempo(2),npasos+1);}\\
\texttt{$>$   h=(intiempo(2)-intiempo(1))/npasos;}\\
\texttt{$>$ soly=ones(npasos+1:1); soly(1)=valorini;}\\
\texttt{$>$   for i=2:npasos+1}\\
\texttt{$>$ ti1=solt(i-1);}\\
\texttt{$>$ ti2=solt(i-1)+0.5.*h;}\\
\texttt{$>$ ti3=solt(i);}\\
\texttt{$>$    yi1=soly(i-1);}\\
\texttt{$>$   yi2=soly(i-1)+h.*f(ti1,yi1);}\\
\texttt{$>$    yi3=soly(i-1)+h.*(-f(ti1,yi1)+2.*f(ti2,yi2));}\\
\texttt{$>$soly(i)=soly(i-1)+(1/6).*h.*(f(ti1,yi1)+4.*f(ti2,yi2)+...}\\
\texttt{$>$         f(ti3,yi3));}\\
\texttt{$>$ end}\\


\subsubsection{Ejemplo de resoluci\'on de un sistema de edos con E. expl\'{\i}cito. Ejercicio 3.30 del libro de problemas}

\texttt{$>$  odefun1=@(t,x,y) y; odefun2=@(t,x,y)-6.54.*x-0.8.*y;}

\texttt{$>$   inicial1=0 inicial2=1; h=0.05; tspan=10;
n=1+(tspan/h);}

\texttt{$>$  t=ones(n); x=ones(n); y=ones(n);}

\texttt{$>$  x(1)=inicial1; y(1)=inicial2;}

\texttt{$>$  t(1)=0;}

\texttt{$>$   for i=2:n;}

\texttt{$>$   t(i)=h.*(i-1);}

\texttt{$>$  x(i)=x(i-1)+h*odefun1(t(i-1),x(i-1),y(i-1));}

\texttt{$>$   y(i)=y(i-1)+h*odefun2(t(i-1),x(i-1),y(i-1));}

\texttt{$>$  end}

\texttt{$>$   plot(t,x,'r'); legend('z');}

\newpage

\subsection{ C\'odigos para resolver problemas de
contorno }

\subsubsection{M\'etodo para resolver un problema de contorno 1D con condiciones Dirichlet}

\texttt{$>$ function
[xh,uh]=bvpdirichlet(a,b,numeronodos,D,V,Q,f,ua,ub)}


\textit{Este c\'odigo resuelve un problema de contorno
unidimensional con condiciones de contorno Dirichlet:}

     -D u''+Vu'+Qu=f en (a,b)

\textit{con condiciones de contorno u(a)=ua, u(b)=ub, utilizando
diferencias centradas. El resto de argumentos son:
\begin{itemize}
\item numeronodos: el n\'umero de nodos inclu\'{\i}dos los extremos del
intervalo.
\item xh: el vector que contiene los nodos.
\item uh: es la soluci\'on num\'erica.
\end{itemize}}


\texttt{$>$ h = (b-a)/(numeronodos-1);}
%n�mero de intervalos=n�mero de nodos-1

\texttt{$>$ xh = (linspace(a,b,numeronodos));}

\texttt{$>$ nodosinteriores=numeronodos-2;}
%construimos la matriz de coeficientes que resulta de aplicar el esquema
%centrado

\texttt{$>$ hD = D/ h\^\,2 ; hV = V/ (2*h);}

\texttt{$>$ e=ones(nodosinteriores,1);}

%matriz de coeficientes A

\texttt{$>$ A = spdiags([-hD*e-hV (2*hD+Q)*e
-hD*e+hV],...}\\
\texttt{$>$ -1:1,nodosinteriores,nodosinteriores);}\\

\texttt{$>$ xi = xh(2:end-1);}
%valores de nodos interiores
%lado derecho del sistema ff
\texttt{$>$ ff =feval(f,xi); ff(1) =  ff(1)+ua*(hD+hV);}

\texttt{$>$ ff(end) = ff(end)+ub*(hD-hV);}

\texttt{$>$ uh = A$\backslash$ ff; uh=[ua; uh; ub];}

\texttt{$>$ return}

\subsubsection{M\'etodo para resolver un problema de contorno 1D con condiciones mixtas}

\texttt{$>$function[xh,uh]=bvp2cvrobinup(a,b,N,D,V,q,bvpfun,c11,...}

\texttt{$>$c12,c21,c22,ua,ub,esquema,varargin)}

\textit{Este c\'odigo resuelve problemas de contorno
unidimensionales de la forma}
$$ -D(x)*(d^2U/dX^2)+V(x)*dU/dX+q(x)*U=BVPFUN$$
\textit{en el intervalo (A,B) con las condiciones de contorno}
$$c11*U'(A)+c12*U(A)=UA$$
$$c21*U'(B)+c22*U(B)=UB$$
\textit{mediante difirencias finitas en N nodos internos
equiespaciado en (A,B). Notaci\'on:
\begin{itemize}
\item BVPFUN, D, V, q:  funciones inline o an\'onimas.
\item  esquema = 'C', se usan f\'ormulas centradas.
\item  esquema = 'U', se usan formulas descentradas a
contracorriente ("upwind") para el t\'ermino convectivo.
\item varargin: permite pasar par\'ametros adicionales a la funci\'on
BVPFUN.
\item  xh: contiene los nodos de la discretizaci\'on.
\item uh: contiene la soluci\'on num\'erica.
\item N: es el n\'umero de nodos interiores.
\item c11,c12,c21,c22: coeficientes en las condiciones de
contorno.
\item ua, ub: valores de las condiciones de contorno.
\end{itemize}}


\texttt{$>$ h = (b-a)/(N+1);}\\

\texttt{$>$   xh = (linspace(a,b,N+2))'; xi = xh(2:N+1);}\\

\texttt{$>$ if (c11 == 0 \&\& c12 == 0)}\\

\texttt{$>$ disp('Condici\'on en la frontera izquierda
incorrecta')}\\

\texttt{$>$  return}\\

\texttt{$>$ end}\\

\texttt{$>$ if (c21 == 0 \&\& c22 == 0)}\\

\texttt{$>$  disp('Condici\'on de contorno en la frontera derecha
incorrecta')}\\

\texttt{$>$  return}\\

\texttt{$>$ end}\\

%% Ecuaciones en la frontera
\texttt{$>$ if c11 ~= 0}\\

\texttt{$>$ TrIzq=[-D(a)/(h\^\,2)\, 2*D(a)/(h\^\,2)-V(a)*c12/c11+q(a)\, -D(a)/(h\^\,2)];}\\

\texttt{$>$  FlIzq = [-c11/(2*h)\, c12\, c11/(2*h)];}\\

\texttt{$>$ end}\\

\texttt{$>$ if c21 ~= 0}\\

\texttt{$>$ TrDcha=-D(b)/(h\^\,2)\,
2*D(b)/(h\^\,2)-V(b)*c22/c21+q(b)\,-D(b)/(h\^\,2)];}\\


\texttt{$>$ FlDcha = [-c21/(2*h)\, c22\, c21/(2*h)];}\\

\texttt{$>$ end}\\

\texttt{$>$ if esquema == 'C'} \\   % Esquema centrado

\texttt{$>$  nW = -D(xi)/(h\^\,2)-V(xi)/(2*h);}\\

\texttt{$>$  nC = 2*D(xi)/(h\^\,2)+q(xi);}\\

\texttt{$>$  nE = -D(xi)/(h\^\,2)+V(xi)/(2*h);}\\

 \texttt{$>$ K = [nW nC nE];}\\

\texttt{$>$ elseif esquema == 'U'} \\   % Esquema descentrado
%  disp('Esquema upwind para la parte convectiva')

\texttt{$>$  for i=2:N+1 }\\ % parametro rho

\texttt{$>$ if V(a+(i-1)*h) ~= 0}\\

\texttt{$>$ rho(i-1) = 1/2+abs(V(a+(i-1)*h))./(2*V(a+(i-1)*h));}

\texttt{$>$ else}\\

\texttt{$>$ rho(i-1) = 1/2;}\\

\texttt{$>$ end}\\

 \texttt{$>$ end}\\

\texttt{$>$  nW = -D(xi)/(h\^\,2)-V(xi).*rho'/h;}\\

\texttt{$>$  nC = 2*D(xi)/(h\^\,2)+V(xi).*(2*rho'-1)/h+q(xi);}\\

\texttt{$>$  nE = -D(xi)/(h\^\,2)+V(xi).*(1-rho')/h;}\\

\texttt{$>$  K = [nW nC nE];}\\

\texttt{$>$ else}\\

 \texttt{$>$ disp('Error en el tipo de esquema. Las opciones
 son:')}\\

 \texttt{$>$ disp('C = centrado')}\\

\texttt{$>$  disp('U = upwind')}\\

\texttt{$>$  return}\\

\texttt{$>$ end}\\
% Construccion de las matrices para cada caso

\texttt{$>$ if c11 ~= 0}\\

%  disp('Cond. tipo Robin en Frontera izquierda')
\indent\texttt{$>$  if c21 ~= 0}\\
    %%% Caso Robin izq, Robin dcha
%    disp('Cond. tipo Robin en Frontera derecha')

\indent \texttt{$>$   A = full([spdiags(FlIzq,[0 1 2],1,N+4);...}\\

 \texttt{$>$ spdiags(TrIzq,[0 1 2],1,N+4);...}\\

 \texttt{$>$ spdiags(K,[1 2 3],N,N+4);...}\\

 \texttt{$>$ spdiags(TrDcha,[N+1 N+2 N+3],1,N+4);...}\\

 \texttt{$>$ spdiags(FlDcha,[N+1 N+2 N+3],1,N+4)]);}\\

    %B = [ua; bvpfun(a)-V(a)*ua/c11; feval(bvpfun,xi); bvpfun(b)-V(b)*ub/c21; ub];  %sin varargin

 \texttt{$>$ B = [ua; bvpfun(a)-V(a)*ua/c11;...}\\

 \texttt{$>$ feval(bvpfun,xi,varargin{:});...}\\

 \texttt{$>$ bvpfun(b)-V(b)*ub/c21; ub];}\\  %con varargin

 \texttt{$>$ uh = A$\backslash$ B;}\\

 \texttt{$>$ uh = [uh(2:N+3)];}\\

\texttt{$>$  else}\\
    %%% Caso Robin izq, Dirichlet dcha
%    disp('Cond. tipo Dirichlet en Frontera derecha')

 \texttt{$>$   cD = zeros(1,N+3); cD(1,end) = 1;}\\

 \texttt{$>$   A = full([spdiags(FlIzq,[0 1 2],1,N+3);...}\\

 \texttt{$>$ spdiags(TrIzq,[0 1 2],1,N+3);...}\\

 \texttt{$>$ spdiags(K,[1 2 3],N,N+3); cD]);}\\

    %B = [ua; bvpfun(a)-V(a)*ua/c11; feval(bvpfun,xi); ub/c22];  %sin varargin

 \texttt{$>$ B = [ua; bvpfun(a)-V(a)*ua/c11;feval(bvpfun,xi,varargin{:}); ub/c22];}\\  %con varargin

\texttt{$>$ uh = A$\backslash$ B;}\\

\texttt{$>$ uh = [uh(2:N+2); ub];}\\

\texttt{$>$  end}\\

\texttt{$>$ else}\\

%  disp('Cond. tipo Dirichlet en Frontera izquierda')
 \texttt{$>$ if c21 ~= 0}\\

    %%% Caso Diriclet izq, Robin dcha
%    disp('Cond. tipo Robin en Frontera derecha')

\texttt{$>$    cI = zeros(1,N+3); cI(1,1) = 1;}\\

 \texttt{$>$   A = full([cI; spdiags(K,[0 1 2],N,N+3);...}\\

 \texttt{$>$ spdiags(TrDcha,[N N+1 N+2],1,N+3);...}\\

 \texttt{$>$ spdiags(FlDcha,[N N+1 N+2],1,N+3)]);}\\

    %B = [ua/c12; feval(bvpfun,xi); bvpfun(b)-V(b)*ub/c21; ub];  %sin varargin
 \texttt{$>$   B = [ua/c12; feval(bvpfun,xi,varargin{:});...}\\

 \texttt{$>$ bvpfun(b)-V(b)*ub/c21; ub];}\\
   %con varargin

\texttt{$>$    uh = A$\backslash$ B;}\\

\texttt{$>$ uh = [ua; uh(2:N+2)];}\\

\texttt{$>$ else}\\

    %%% Caso Diriclet izq, Dirichlet dcha
%    disp('Cond. tipo Dirichlet en Frontera derecha')

 \texttt{$>$   cI = zeros(1,N+2); cI(1,1) = 1;}\\

 \texttt{$>$ cD = zeros(1,N+2); cD(1,end) = 1;}\\

\texttt{$>$    A = full([cI; spdiags(K,[0 1 2],N,N+2); cD]);}\\

    %B = [ua/c12; feval(bvpfun,xi); ub/c22];  %sin varargin

\texttt{$>$ B = [ua/c12; feval(bvpfun,xi,varargin{:}); ub/c22];} \\ %con varargin

\texttt{$>$ uh = A$\backslash$ B;}\\

\texttt{$>$    uh = [ua; uh(2:N+1); ub];}\\

 \texttt{$>$ end}\\

\texttt{$>$ end}\\


\texttt{$>$ return}\\



\subsubsection{M\'etodo para resolver la ecuaci\'on del calor 1D}

\texttt{$>$ function[x,uf]=ecucalor(C,intespacio,intiempo,...}\\

\texttt{$>$ pasosespacio,pasostiempo,theta,u0,cc,f)}\\

\textit{Este c\'odigo resuelve la ecuaci\'on del calor: en una
dimensi\'on espacial utilizando un theta-m\'etodo para la
discretizaci\'on temporal:}

$$\frac{du}{dt}=C \frac{d^2 u}{dx^2} +f, \quad (t,x)\in  (t_0,T)\times(a,b),$$
$$u(x,t)=cc(x,t),\quad x=a,\quad x=b,$$
$$u(x,t_0) = u_0(x),\quad  x \in [a,b],$$
\textit{en una dimensi\'on espacial utilizando un theta-m\'etodo
para la discretizaci\'on temporal.}

\textit{Los argumentos son: \begin{itemize}
\item intespacio: (a,b).
\item intiempo: (t0,T).
\item u0: condici\'on inicial.
\item cc: condici\'on de contorno tipo Dirichlet.
\item theta: el valor de theta empleado en el theta m\'etodo,
theta=0 es para Euler expl\'{\i}cito, theta=1 es para Euler
impl\'{\i}cito, theta=0.5 para Crank-Nicolson.
\item pasosespacio: el n\'umero de pasos en espacio.
\item pasostiempo: el n\'umero de pasos en tiempo.
\item x: el vector que contiene los nodos en espacio.
\item uf: contiene los valores num\'ericos obtenidos en t=T.
\end{itemize}}

\texttt{$>$ h  = (intespacio(2)-intespacio(1))/pasosespacio;}\\

\texttt{$>$ dt = (intiempo(2)-intiempo(1))/pasostiempo;}\\

\texttt{$>$ N = pasosespacio+1;}

%n�mero de nodos en espacio
%perparamos la matriz de coeficientes, que variar� dependiendo del
%theta-m�todo empleado

\texttt{$>$ e =ones(N,1);} \\

\texttt{$>$ D = spdiags([-e 2*e -e],[-1,0,1],N,N); I =
speye(N);}\\

\texttt{$>$ M = I+C*dt*theta*D/ h\^\,2; Mn =
I-C*dt*(1-theta)*D/h\^\,2;}

\texttt{$>$ M(1,:) = 0;
 M(1,1) = 1; M(N,:) = 0; M(N,N) = 1;}

\texttt{$>$ x =linspace(intespacio(1),intespacio(2),N);}

\texttt{$>$ x = x'; fn = feval(f,intiempo(1),x);}


\texttt{$>$ un =feval(u0,x);}

\texttt{$>$ [L,U]=lu(M);}

\texttt{$>$ for t = intiempo(1)+dt:dt:intiempo(2)}

\texttt{$>$    fn1 = feval(f,t,x);}

    %lado derecho
 \texttt{$>$   rhs = Mn*un+dt*(theta*fn1+(1-theta)*fn);}

\texttt{$>$    temp = feval(cc,t,[intespacio(1),intespacio(2)]);}

\texttt{$>$ rhs([1,N]) = temp;}
    %el sistema se resuleve utilizando una factorizaci�n LU
 \texttt{$>$   u = L$\backslash$ rhs;}

\texttt{$>$   u = U$\backslash$u;}

\texttt{$>$ fn = fn1;    un = u;}

\texttt{$>$ end}

\texttt{$>$ uf=u;}

\texttt{$>$ return}


\subsubsection{M\'etodo para resolver la ecuaci\'on de Poisson 2D}

\texttt{$>$function [u,x,y,error]=ecupoisson(a,b,c,d,dx,dy,f,...}

\texttt{$>$                       cc,solexacta)}


\textit{Este c\'odigo resuelve la ecuaci\'on de Poisson
bidimensional:}  $$-\Delta(u)=f (x,y),\,\,\mbox{en}\,\,
(a,b)\times(c,d)$$ \textit{con condiciones de frontera Dirichlet,
dadas por la funci\'on cc(x,y), utilizando el esquema de cinco
punto en cruz para aproximar el operador laplaciano $\Delta$.
Tambi\'en calcula el error cometido cuando se conoce la soluci\'on
exacta, definida en solexacta(x,y). En caso de no conocer la
soluci\'on exacta, no se mete como argumento. Los tama\~nos de
discretizaci\'on vienen dados por dx y dy.}

\noindent\texttt{$>$  if nargin == 8}\\
 \texttt{$>$ solexacta=@(x,y) 0.*x;}\\
\texttt{$>$  end}\\

\textit{N\'umero de intervalos en x, nx, en y, ny}

\noindent\texttt{$>$  nx=(b-a)/dx; ny=(d-c)/dy;}\\
\texttt{$>$ nx1=nx+1;}\\
\texttt{$>$ dx2=dx\^ \,2; dy2=dy\^ \,2;}\\
\texttt{$>$ kii=2/dx2+2/dy2;     kix=-1/dx2;  kiy=-1/dy2;}\\
\texttt{$>$ dim=(nx+1)*(ny+1);   K=speye(dim,dim);}\\
\texttt{$>$ rhs=zeros(dim,1);}\\
\texttt{$>$ rhs1=zeros(dim,1);}\\
\texttt{$>$ y = c;}\\

\textit{Calculamos la matriz de coeficientes del sistema}

\noindent\texttt{$>$ for m = 2:ny}\\
\texttt{$>$ x = a; y = y + dy;}\\
\texttt{$>$ for n = 2:nx}\\
\texttt{$>$ i = n+(m-1)*(nx+1);}\\
\texttt{$>$ x = x + dx;}\\
\texttt{$>$ rhs(i) = feval(f,x,y);}\\
\texttt{$>$ K(i,i) = kii;         K(i,i-1) = kix;}\\
\texttt{$>$ K(i,i+1) =kix; K(i,i+nx1)=kiy;}\\
\texttt{$>$ K(i,i-nx1)=kiy;}\\
\texttt{$>$ end}\\
\texttt{$>$ end}\\

\textit{Calculamos el lado derecho del sistema de ecuaciones}

\noindent\texttt{$>$ rhs1(1:nx1) = feval(cc,x,c);}\\
\texttt{$>$ rhs1(dim-nx:dim) = feval(cc,x,d);}\\
\texttt{$>$ y = [c:dy:d];}\\
\texttt{$>$ rhs1(1:nx1:dim-nx) = feval(cc,a,y);}\\
\texttt{$>$ rhs1(nx1:nx1:dim) = feval(cc,b,y);}\\
\texttt{$>$ rhs = rhs - K*rhs1;}\\
\texttt{$>$ nbound = [[1:nx1],[dim-nx:dim],...}\\
\texttt{$>$ [1:nx1:dim-nx],[nx1:nx1:dim]];}\\
\texttt{$>$ ninternal = setdiff([1:dim],nbound);}\\
\texttt{$>$ K = K(ninternal,ninternal);}\\

\textit{Resolvemos el sistema para los nodos interiores}

\noindent\texttt{$>$ rhs = rhs(ninternal);}\\
\texttt{$>$ utemp = K$\backslash$ rhs;}\\
\texttt{$>$ uh = rhs1;}\\

\textit{Imponemos las condiciones de contorno en los bordes}

\noindent\texttt{$>$ uh (ninternal) = utemp;}\\
\texttt{$>$ k = 1; y = c;}\\
\texttt{$>$ for j = 1:ny+1}\\
\texttt{$>$ x = a;}\\
\texttt{$>$ for i = 1:nx1}\\
\texttt{$>$ u(i,j) = uh(k);}\\
\texttt{$>$  k = k + 1;}\\
\texttt{$>$  ue(i,j) = feval(solexacta,x,y);}\\
\texttt{$>$  x = x + dx;}\\
\texttt{$>$ end}\\
\texttt{$>$ y = y + dy;}\\
\texttt{$>$ end}\\
\texttt{$>$ x = [a:dx:b];}\\
\texttt{$>$ y = [c:dy:d];}\\
\texttt{$>$ if nargout == 4}\\
\texttt{$>$ if nargin == 8}\\
\texttt{$>$ warning('La soluci\'on exacta no se conoce');}\\
\texttt{$>$ error = [ ];}\\
\texttt{$>$ else}\\
\texttt{$>$ error = max(max(abs(u-ue)))/max(max(abs(ue)));}\\
\texttt{$>$  end}\\
\texttt{$>$ end}\\
\texttt{$>$ return}\\



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
\newpage
\section{SEMINARIOS: EJERCICIOS RESUELTOS }
\vspace{1cm}


{\bf{Observaciones:}} Toda l\'{\i}nea que comience por \%
indicar\'a un comentario y lo que siga al s\'{\i}mbolo prompt del
sistema, \texttt{$>$}, es una l\'{\i}nea de comandos.

\vspace{0.5cm}


\begin{enumerate}
\item Sea dada la funci\'on $$f(x) = x^3 + 4x^2-10,$$
\begin{enumerate}
\item Aplicar el algoritmo de bisecci\'on para calcular todas las
soluciones de la ecuaci\'on $f(x) = 0$, trabajando con una
tolerancia de $10^{-3}$ y un n\'umero m\'aximo de 100 iteraciones.
\item Aplicar el m\'etodo de Newton para resolver la misma ecuaci\'on
con los mismos valores de tolerancia y n\'umero m\'aximo de
iteraciones. \end{enumerate}

\vspace{1cm} {\bf{Soluci\'on.}} Escribimos en un script el
siguiente texto:

\% Apartado a)\\

\texttt{$>$ fecu=@(x) x.\^ \, 3+4.*x.\^ \, 2-10;}\\
\texttt{$>$ a=0;b=2;maxitera=100;errorper=1e-3;}\\

\% Comprobamos que hemos elegido
bien los extremos del intervalo:\\

\texttt{$>$ fecu(a); fecu(b)}\\

\% y vemos que toman signos distintos, -10 y 14.\\

\% Hacemos la llamada al c\'odigo que resuelve el problema por
bisecci\'on, \\
\% metbiseccion.m: \\

\texttt{$>$
[sol,itera]=metbiseccion(fecu,a,b,errorper,maxitera)};\\

\% Apartado b).  A continuaci\'on resolvemos el problema por el
m\'etodo de
Newton-Raphson:\\

\% s\'olo queda por definir:\\

\texttt{$>$ dfecu=@(x) 3.*x.\^ \, 2+8.*x; x0=1;}\\

\% y llamamos al c\'odigo metnewton1ecu.m


\texttt{$>$ [sol2,itera2]=metnewton1ecu(fecu,dfecu,x0,errorper,maxitera)}\\

\%Fin del script.\\

Las soluciones obtenidas son: sol=1.3643, itera=10, con el
m\'etodo de bisecci\'on y sol2=1.3652, itera2=4, con el m\'etodo
de Newton.

\vspace{1cm}




\item Sea dada la funci\'on
$$f(x) = e^x-3x^2.$$
\begin{enumerate}
\item Considerando el esquema num\'erico asociado al m\'etodo de
Newton como un esquema de punto fijo, definir la funci\'on $g(x)$
asociada al m\'etodo y utilizar el algoritmo de punto fijo para
calcular las soluciones de la ecuaci\'on $f(x) = 0$, trabajando
con una tolerancia de $10^{-6}$ y un n\'umero m\'aximo de 1000
iteraciones.
\end{enumerate}

\vspace{1cm} {\bf{Soluci\'on.}} Escribimos en un script el
siguiente texto:

\texttt{$>$ f=@(x) exp(x)-3.*x.\^ \, 2;}\\

\% Comprobamos que en el intervalo [0,1] hay una soluci\'on


\texttt{$>$ a=0;b=1; f(a),f(b)}\\


\% f(a)=1, f(b)=-0.2817

\% En efecto, la funci\'on toma signos distintos en 0 y en 1: 1 y -0.28172, resp.\\

\%Definimos g(x)=x-f/f'\\

\texttt{$>$ g=@(x) x-(exp(x)-3.*x.\^ \, 2).*(exp(x)-6.*x).\^ \, (-1);}\\

\texttt{$>$ x0=0.5;errorper=1e-6;maxitera=1000;}\\

\% Hacemos la llamada al c\'odigo metpuntofijo.m\\

\texttt{$>$ [x,itera]=metpuntofijo(g,x0,errorper,maxitera)}\\

La soluci\'on obtenida es: x=0.91001 e itera=5. \vspace{1cm}




\item Resolver el sistema de ecuaciones no lineales:
\[
\left\{ \begin{array}{l} \mbox{ln}(x^2 + y^2)-\mbox{sin}(xy)- (\mbox{ln}( 2) + \mbox{ln}(\pi)) = 0\\
e^{x-y}+ \mbox{cos}(xy) = 0. \end{array} \right. \] utilizando el
m\'etodo de Newton y trabajando con una tolerancia de $10^{-6}$.
Utilizar como semilla $(x_0,y_0) = (2,2)$. \vspace{1cm}

{\bf{Soluci\'on}}.

Primero definimos los scripts con las funciones fecusistema.m y
jacobiana.m:


\texttt{$>$ function F=fecusistema(x,y)}\\

\texttt{$>$ F(1,1)=log(x.\^\, 2+y.\^\, 2)-sin(x.*y)-(log(2)+log(pi));}\\
\texttt{$>$ F(2,1)=exp(x-y)+cos(x.*y);}\\
\texttt{$>$ end}\\

\texttt{$>$ function J=jacobiana(x,y)}\\
\texttt{$>$ J(1,1)=2.*x.*(x.\^\, 2+y.\^\, 2).\^\, (-1)-y.*cos(x.*y);}\\
\texttt{$>$ J(1,2)=2.*y.*(x.\^\, 2+y.\^\, 2).\^\, (-1)-x.*cos(x.*y);}\\
\texttt{$>$ J(2,1)=exp(x-y)-y.*sin(x.*y);}\\
\texttt{$>$ J(2,2)=-exp(x-y)-x.*sin(x.*y);}\\
\texttt{$>$ end}\\

A continuaci\'on, llamamos al c\'odigo metnewtonsistemas.m, para
ellos escribimos un nuevo script con los comandos:

\texttt{$>$ vectorx0=[2;2]; errorper=1e-6;maxitera=100;}\\
\texttt{$>$ [X,itera]=metnewtonsistema(@fecusistema,...}\\
\texttt{$>$ @jacobiana,vectorx0,errorper,maxitera)}\\

La soluci\'on obtenida es: X=(x,y)=(1.77245,1.77245), itera=6.
\vspace{1cm}


\item Se considera el problema de valor inicial:
\[ \left\{
\begin{array}{l}
y^{'}=y-\mbox{sin}(t)+\mbox{cos(t)},\\
y(0)=1.
\end{array}\right.\]
En el intervalo temporal $[0,2]$, utilizar los esquemas
num\'ericos de Euler expl\'{\i}cito, Euler impl\'{\i}cito y
Crank-Nicolson para obtener una soluci\'on aproximada, tomando un
paso de discretizaci\'on $h=0.1$. Sabiendo que la soluci\'on
exacta es $y(t)=e^t+sin(t)$, dibujar las soluciones obtenidas con
los tres m\'etodos y la soluci\'on exacta en un mismo plot.
Comparar los resultados obtenidos. Escribir los valores de las
distintas soluciones en $t=0.5$ y $t=1.5$. \vspace{1cm}

{\bf{Soluci\'on.}} Abrimos un nuevo script para escribir el
siguiente c\'odigo y luego ejecutarlo:

\texttt{$>$ f=@(t,y) y-sin(t)+cos(t);}\\

\texttt{$>$ intiempo=[0 2]; valorini=1;npasos=20;}\\

\% Hacemos una llamada a los c\'odigos correspondientes:

\texttt{$>$ [solt,solye]=eulerexpicito(f,intiempo,valorini,npasos);}\\

\texttt{$>$ [solt,solyi]=eulerimplicito(f,intiempo,valorini,npasos);}\\

\texttt{$>$ [solt,solyc]=cranknicolson(f,intiempo,valorini,npasos);}\\

\% Observa que llamamos de forma distinta a los resultados para y.

\texttt{$>$ solexac=@(t) exp(t)+sin(t); figure;}\\

\% Dibujamos las soluciones junto con los valores de la soluci\'on
exacta para comparar

\texttt{$>$ plot(solt,solye,'ro',solt,solyi,'b+',...}\\

\texttt{$>$ solt,solyc,'g*',solt,solexac(solt),'k\^\,')}\\



\texttt{$>$ t1=min(find(solt>=0.5));}\\

\% tambi\'en t1=find(solt==0.5);\\

\texttt{$>$ t2=min(find(solt>=1.5));}\\

\% tambi\'en t2=find(solt==1.5);\\

\texttt{$>$ solye(t1),solyi(t1),solyc(t1), solye(t2),solyi(t2),solyc(t2)}\\


Las soluciones son: para t=0.5, con EE 2.096, con EI 2.1646 y con
CN 2.1283. Para t=1.5, con EE 5.2516, con EI 5.7585 y con CN,
5.4825.



\begin{figure}[h]
 \centering
       \includegraphics[width=0.6\textwidth]{figsem4.eps}
    \caption{Resultados obtenidos con los 3 m\'etodos junto con la exacta (en c\'{\i}rculos rojos EE, en + azules EI, en * verdes CN, y la exacta en tri\'angulos negros).}
%   % \label{figura101}
\end{figure}

\vspace{1cm}

\item Se considera el problema de valor inicial:
\[ \left\{
\begin{array}{l}
y^{'}=-2ty^2,\\
y(0)=1.
\end{array}\right.\]
Aplicar los m\'etodos tipo Runge-Kutta dados por los algoritmos de
Heun (heun.m) y Simpson (rungekuttao3.m), para resolver el PVI
anterior en el intervalo temporal [0,1], tomando como paso de
discretizaci\'on $h=0.01$. Escribir los valores obtenidos con los
dos m\'etodos para los tiempos $t=0.1$, $t=0.5$ y $t=1$.

\vspace{1cm} {\bf{Soluci\'on}}. Abrimos un script con el siguiente
c\'odigo :

\texttt{$>$ f=@(t,y) -2.*t.*y.\^\, 2; intiempo=[0 1]; }\\

\texttt{$>$ npasos=100;valorini=1;}\\

\texttt{$>$ [solt,soly1]=heun(f,intiempo,valorini,npasos);}\\

\texttt{$>$ [solt,soly2]=rungekuttao3(f,intiempo,valorini,npasos);}\\

\% Dibujamos los resultados obtenidos con los dos m\'etodos:

\texttt{$>$ plot(solt,soly1,'ro',solt,soly2,'b*')}\\

\% Calculamos los valores num\'ericos en los tiempos pedidos:

\texttt{$>$ t1=min(find(solt>=0.1));}\\ % tambi�n t1=find(tt>=0.1)
\texttt{$>$ t2= min(find(solt>=0.5));}\\ %t2=find(tt>=0.5)
\texttt{$>$ t3=min(find(solt>=1));}\\ %t3=find(tt>=1)
\texttt{$>$ soly1(t1),soly1(t2),soly1(t3),soly2(t1),soly2(t2),soly2(t3)}\\

Los resultados obtenidos con Heun son: en t=0.1, 0.9901, en t=0.5,
0.8000, en t=1, 0.5174. Los resultados obtenidos con Simpson son:
en t=0.1, 0.9901, en t=0.5, 0.8006, en t=1, 0.5014.

\begin{figure}[h]
 \centering
       \includegraphics[width=0.6\textwidth]{sem5.eps}
    \caption{Resultados obtenidos con el m\'etodo de Heun en rojo y con Simpson en azul.}
%   % \label{figura101}
\end{figure}


\item Sea dado el Problema de Transporte Difusivo y Convectivo
definido por el PVC:

\[ \left\{
\begin{array}{l}
-u^{''}-4u^{'}=-16x^3+34x-1,\quad \forall x\in (0,2),\\
u(0)=4,\quad u(2)=2,\end{array}\right.\] cuya soluci\'on
anal\'{\i}tica es: $$sol(x) = x^4-x^3-3.5x^2 + 2x + 4,$$ siendo u
la concentraci\'on de un contaminante. \begin{enumerate}
\item Aplicar el algoritmo bvpdirichlet.m para calcular la soluci\'on en el
intervalo [0, 2] con paso de discretizaci\'on $h = 0.125$. \item
Dibujar la soluci\'on anal\'{\i}tica junto con la soluci\'on
num\'erica. \item Determinar los valores m\'aximos y m\'{\i}nimos
de contaminaci\'on as\'{\i} como su localizaci\'on utilizando la
soluci\'on num\'erica y la anal\'{\i}tica. 
\end{enumerate}

\vspace{1cm} {\bf{Soluci\'on}} Creamos un script con el siguiente
c\'odigo para luego ejecutarlo:

\texttt{$>$ a=0; b=2; D=1; V=-4; Q=0;}\\

\texttt{$>$ f=@(x) -16.*x.\^\, 3+34.*x-1; ua=4; ub=2; NumeroNodos=17;}\\

\texttt{$>$ [xh,uh]=bvpdirichlet(a,b,NumeroNodos,D,V,Q,f,ua,ub);}\\

\texttt{$>$ solexac=xh.\^\, 4-xh.\^\, 3-3.5.*xh.\^\, 2+2.*xh+4;}\\

\texttt{$>$ figure; plot(xh,uh,'r',xh,solexac,'g')}\\

\begin{figure}[h]
 \centering
       \includegraphics[width=0.6\textwidth]{fig1sem6.eps}
    \caption{Resultados obtenidos con el m\'etodo num\'erico en rojo, y la soluci\'on exacta en verde.}
%   % \label{figura101}
\end{figure}


\texttt{$>$ nmax=max(uh),nmin=min(uh)}\\

\texttt{$>$ emin=min(solexac),emax=max(solexac)}\\

Los resultados son: nmax=4.3235, nmin=0.7233, emin=0.6897,
emax=4.2695.

\% Dibujamos el error cometido:

\texttt{$>$ figure; err=uh-solexac; plot(xh,err)}\\

\begin{figure}[h]
 \centering
       \includegraphics[width=0.6\textwidth]{figsem6.eps}
    \caption{Error cometido.}
%   % \label{figura101}
\end{figure}


\% Vemos d\'onde se alcanzan el valor m\'aximo y el m\'{\i}nimo


\texttt{$>$ xm=find(uh>=nmax), xmin=find(uh<=nmin),
xh(xm),xh(xmin)}\\

Los valores obtenidos (componentes) son: xm=3 y xmin=14. Por
tanto, el valor m\'aximo se alcanza en x=0.25, y el valor
m\'{\i}nimo en x=1.6250.

\vspace{1cm}



\item Sea dado el PVIC:
\[ \left\{
\begin{array}{l}
\frac{\partial u}{\partial t} = \frac{\partial ^2 u}{\partial x^2}
+ f(x,t),\\
u(x,t)=cc(x,t),\quad x=0,\quad x=1,\\
u(x,0) = u_0(x),\quad  x \in(0,1)\end{array}\right.\] donde
$$f(x,t) = 2t-x^2 + 10xt,$$
$$ cc(x,t) = x^2(x-1),$$
$$u_0(x) = 0.$$
Aplicar el algoritmo definido en ecucalor.m para calcular la soluci\'on en
el intervalo temporal [0,1], con paso de discretizaci\'on
espacial, $\Delta x = 0.05$ y de discretizaci\'on temporal,
$\Delta t = 0.02$. Dibujar la soluci\'on obtenida para $\theta= 0
$ (Euler Expl\'{\i}cito). ?`Puedes justificar la gr\'afica
obtenida ? Determinar el paso de discretizaci\'on temporal
m\'aximo que asegura estabilidad y calcular la soluci\'on para ese
paso temporal.  ?`Cu\'anto vale la soluci\'on en el punto $x$ =
0.85 ?

\vspace{1cm} {\bf{Soluci\'on.}} Abrimos un nuevo script con el
siguiente c\'odigo:

\texttt{$>$ intespacio=[0 1]; intiempo=[0 1];}\\

\texttt{$>$ pasosespacio=20,pasostiempo=100;}\\

\texttt{$>$ theta=0; c=1; u0=@(x) 0.*x;}\\

\texttt{$>$ cc=@(t,x) (x.\^\,2).*(x-1);}\\

\texttt{$>$ f=@(t,x) 2.*t-x.\^\, 2+10.*x.*t;}\\

\texttt{$>$ [x,u]=ecucalor(c,intespacio,intiempo,pasosespacio,...}\\

\texttt{$>$ pasostiempo,theta,u0,cc,f);}\\

\texttt{$>$ figure; plot(x,u)}\\

\% Se puede observar que el plot no aparece puesto que hay
valores\\
\% que son del orden de 1e+11


\texttt{$>$ x1=min(find(x>=0.85)); u(x1)}\\

\% Se tiene que x1=18 y u(18)=3.6504e+11

\% Esto se debe a la inestabilidad asociada al m\'etodo
expl\'{\i}cito

\% Para que sea estable, debemos reducir el paso temporal

\texttt{$>$ dx=0.05; dt=(0.05)\^\, 2; pasost=1/(0.05)\^\, 2 ;}\\

\texttt{$>$ pasosespacio=20; pasosespacio=800;}


\texttt{$>$ [x,u1]=ecucalor(c,intespacio,intiempo,pasosespacio,pasostiempo,...}\\

\texttt{$>$ theta,u0,cc,f);}\\

\texttt{$>$ figure;}\\

\texttt{$>$ plot(x,u1), x1=min(find(x>=0.85)), u1(x1)}\\
\begin{figure}
 \centering
       \includegraphics[width=0.6\textwidth]{figsem7.eps}
    \caption{Resultados obtenidos cumpliendo el criterio de estabilidad.}
%   % \label{figura101}
\end{figure}

El resultado obtenido es u1(18)=0.4496.




\end{enumerate}



\newpage 



\section{EJERCICIOS DE PR\'ACTICAS}


\vspace{1cm}

\begin{enumerate}
\item Sea dada la funci\'on $$f(x) = e^x-15-arctg(x).$$
\begin{enumerate}
\item Aplicar el algoritmo de bisecci\'on para calcular la
soluci\'on de la ecuaci\'on $f(x) = 0$, trabajando con una
tolerancia de $10^{-3}$ y un n\'umero m\'aximo de 100 iteraciones.
\item Aplicar el m\'etodo de Newton para resolver la misma ecuaci\'on
con los mismos valores de tolerancia y n\'umero m\'aximo de
iteraciones. \end{enumerate}

\item Sea dada la funci\'on
$$f(x) =\sqrt{x}sin(x)-x^3+2.$$
\begin{enumerate}
\item Considerando el esquema num\'erico asociado al m\'etodo de
Newton como un esquema de punto fijo, definir la funci\'on $g(x)$
asociada al m\'etodo y utilizar el algoritmo de punto fijo para
calcular la soluci\'on de la ecuaci\'on $f(x) = 0$, trabajando con
una tolerancia de $10^{-6}$ y un n\'umero m\'aximo de 1000
iteraciones.
\item Define un esquema de punto fijo diferente al considerado en
el apartado anterior justificando la elecci\'on, para calcular la
soluci\'on de la ecuaci\'on $f(x) = 0$ trabajando con una
tolerancia de $10^{-6}$ y un n\'umero m\'aximo de 1000
iteraciones.
\end{enumerate}
\item Resolver el sistema de ecuaciones no lineales:
\[
\left\{ \begin{array}{l} -x-y^3+x^2-2 = 0,\\
1-2x+y-y^2 = 0. \end{array} \right. \] utilizando el m\'etodo de
Newton y trabajando con una tolerancia de $10^{-6}$. Utilizar como
semilla $(x_0,y_0) = (-0.8,-0.4)$.

\item Se considera el problema de valor inicial:
\[ \left\{
\begin{array}{l}
y^{'}=-kt^{\beta}y,
\\ y(0)=1,
\end{array}\right.\] donde $k=4$ y $\beta=0.25$. Resolver el PVI en el intervalo $[0,4]$, utilizando los
m\'etodos de Euler expl\'{\i}cito y Euler impl\'{\i}cito. Analizar
la estabilidad de los dos m\'etodos utilizando para ello distintos
pasos de discretizaci\'on, por ejemplo $h_1=0.5$, $h_2=0.4$,
$h_3=0.2$ y $h_4=0.1$.
\item Se considera el problema de valor inicial:
\[ \left\{
\begin{array}{l}
y^{'}=\frac{2y}{t+1},\\
y(0)=1.
\end{array}\right.\]
En el intervalo temporal $[0,1]$, utilizar los esquemas
num\'ericos de Euler expl\'{\i}cito, Euler impl\'{\i}cito y
Crank-Nicholson para obtener el valor aproximado de la soluci\'on
en $t=1$, tomando un paso de discretizaci\'on $h=0.1$. Sabiendo
que la soluci\'on exacta es $y(t)=(t+1)^2$, dibujar las soluciones
obtenidas con los tres m\'etodos y la soluci\'on exacta en un
mismo plot. Comparar los resultados obtenidos. Escribir los
valores de las distintas soluciones en $t=0.5$ y $t=0.75$.
Encontrar los intervalos temporales para los cuales las distintas
soluciones num\'ericas alcanzan valores superiores a tres.

\item Se considera el problema de valor inicial:
\[ \left\{
\begin{array}{l}
y^{'}=-y(y-1),\\
y(0)=2.
\end{array}\right.\]
Aplicar los m\'etodos tipo Runge-Kutta dados por los algoritmos
heun.m y rungekuttao3.m, para resolver el PVI anterior en el
intervalo temporal [0,1], tomando como paso de discretizaci\'on
$h=0.01$ Escribir los valores obtenidos con los dos m\'etodos para
los tiempos $t=0.1$, $t=0.5$ y $t=1$.
\item Sea dado el Problema de Transporte Difusivo y Convectivo
definido por el PVC:
\[ \left\{
\begin{array}{l}
-u^{''}+u^{'}=-2e^{-x},\quad \forall x\in (0,1),\\
u(0)=2,\quad u(1)=2cosh(1),\end{array}\right.\] cuya soluci\'on
anal\'{\i}tica es: $$sol(x) =e^x+e^{-x}$$.
\begin{enumerate}
\item Aplicar el algoritmo bvpdirichlet.m para calcular la soluci\'on en el
intervalo [0, 1] con paso de discretizaci\'on $h = 0.01$.
\item
Dibujar la soluci\'on anal\'{\i}tica junto con la soluci\'on
num\'erica.
\item Determinar los valores m\'aximos y m\'{\i}nimos
de $u$ y su localizaci\'on, para la soluci\'on num\'erica y la
soluci\'on anal\'{\i}tica. Calcular los errores cometidos.
\end{enumerate}

%\item Sea dado el Problema de Transporte Difusivo y Convectivo
%definido por el PVC:

%\[ \left\{
%\begin{array}{l}
%-u^{''}-4u^{'}=-16x^3+34x-1,\quad \forall x\in (0,2),\\
%u(0)=4,\quad u(2)=2,\end{array}\right.\] cuya soluci\'on
%anal\'{\i}tica es: $$sol(x) = x^4-x^3-3.5x^2 + 2x + 4,$$ siendo u
%la concentraci\'on de un contaminante. \begin{enumerate}
%\item Aplicar el algoritmo bvp.m para calcular la soluci\'on en el
%intervalo [0, 2] con paso de discretizaci\'on $h = 0.125$. \item
%Dibujar la soluci\'on anal\'{\i}tica junto con la soluci\'on
%num\'erica. \item Determinar los valores m\'aximos y m\'{\i}nimos
%de contaminaci\'on as\'{\i} como su localizaci\'on utilizando la
%soluci\'on num\'erica y la anal\'{\i}tica. Calcular los errores
%cometidos. \item Determinar la regi\'on m\'axima de seguridad (
%donde la concentraci\'on del contaminante es menor que uno) es
%decir $u(x) < 1$. Utilizar la soluci\'on anal\'{\i}tica.
%\end{enumerate}
\item Sea dado el PVIC:
\[ \left\{
\begin{array}{l}
\frac{\partial u}{\partial t} = \frac{\partial ^2 u}{\partial x^2}
+ f(x,t),\\
u(x,t)=cc(x,t),\quad x=0,\quad x=2,\\
u(x,0) = u_0(x),\quad  x \in(0,2),\end{array}\right.\] donde
$$f(x,t) =tsin(x),$$
$$ cc(x,t) =0,$$
$$u_0(x) = x, \, \mbox{para}\, 0\le x \le 1 \quad \mbox{y}\quad u_0(x) = 2-x, \, \mbox{para}\, 1< x \le 2 .$$
Aplicar el algoritmo de ecucalor.m para calcular la soluci\'on en
el intervalo temporal [0,1], con paso de discretizaci\'on espacial
$\Delta x = 0.05$ y discretizaci\'on temporal, $\Delta t = 0.02$.
Dibujar la soluci\'on obtenida para $\theta= 0.5 $ en $t=1$.
?`Cu\'anto vale la soluci\'on en el punto $x$ = 1 ?

\end{enumerate}
\end{document}










\end{document}