\documentclass[12pt,twoside,a4paper]{article}
\usepackage{amsmath, amstext, amssymb}
\usepackage{datetime}
\usepackage[most]{tcolorbox}

\newcommand{\bash}[1]{\begin{tcolorbox}\texttt{\$ #1}\end{tcolorbox}}
\newcommand{\bashio}[2]{\begin{tcolorbox}\texttt{\$ #1}\tcblower\texttt{#2}\end{tcolorbox}}
\newcommand{\fn}[2]{\begin{tcolorbox}[enhanced,colback=white,colframe=gray!50!white,attach boxed title to top center={yshift=-3mm,yshifttext=-1mm},title={Funció \texttt{#1}},boxed title style={size=small,colframe=black},colbacktitle=white,coltitle=black]\texttt{#2}\end{tcolorbox}}
\renewcommand{\Bbb}[1]{\mathbb{#1}}
\newcommand{\set}[1]{\left\{\,#1\,\right\}}

\title{Memòria Pràctica 3 - Mètodes Numèrics}
\author{Arnau del Río}
\date{\today}

\begin{document}
\maketitle

\section{Problema}
L'objectiu d'aquesta pràctica és utilitzar mètodes numèrics per calcular la trajectòria de flyby lunar de la missió Artemis II.

Ho farem a través del mètode de Newton per trobar zeros de funcions i el mètode de Nevill d'interpolació numèrica.

A més a més, utilitzarem certes funcions donades per a propagar la trajectòria de la missió.

En particular, utilitzarem el model del problema restringit de tres cossos, en el qual considerem que la missió es propaga en un pla, on dos cossos, la terra i la lluna, amb massa significativa es troben a l'espai i un cos amb massa insignificativa es troba en una òrbita generada per l'atracció gravitatòria de la terra i la lluna.

Els mètodes de Newton i de Neville estan implementats en un fitxer de biblioteca anomenat \texttt{nwtnev.c}.

Les utilitats que desenvoluparem són les següents:
\begin{itemize}
    \item \texttt{detcsig}: Propaga la trajectòria i troba les posicions i velocitats en cada moment.
    \item \texttt{trobalf}: Resol el problema trobant el moment del flyby i la distància a la lluna.
    \item \texttt{interp}: Interpola un conjunt de punts per a obtenir una funció que aproxima la trajectòria.
\end{itemize}

\newpage

\section{Implementació}

\subsection{Mètode de Newton}
El mètode de Newton permet aproximar arrels d'equacions de la forma $f(x) = 0$.

Suposem que $f: [a, b] \to \Bbb{R}$ és de classe $\mathcal{C}^2$ i que existeix $\alpha \in (a, b)$ amb $f(\alpha) = 0$. Aleshores, el mètode és el següent algorisme iteratiu:

\begin{tcolorbox}[title={Mètode de Newton}, colback=red!5!white,colframe=red!75!black]
    Escollim $x_0 \in [a, b]$ aproximació inicial, $\epsilon > 0$ tolerància i $n_{max} > 0$ el nombre màxim d'iteracions.

    $\forall i = 0, 1, \dots, n_{max} - 1$

    \qquad $x_{i+1} = x_i - \frac{f(x_i)}{f'(x_i)}$

    \qquad $\text{si } |x_{i+1} - x_i| < \epsilon, \text{ STOP}$
\end{tcolorbox}

\newpage

\subsection{Algoritme de Neville}
L'algoritme de Neville ens permet trobar el polinomi interpolador de Lagrange a partir d'un conjunt de punts de suport donats.

Considerem que tenim un conjunt de $n+1$ punts de suport $\set{(x_i, f_i)}_{i=0}^n$. Aleshores denotem per $P_{i_0, i_1, \dots, i_k}$ el polinomi interpolador de Lagrange resol $P_{i_0, i_1, \dots, i_k}(x_{i_j}) = f_{i_j}$, per $j = 0, 1, \dots, k$.

\begin{tcolorbox}[title={Algoritme de Neville}, colback=blue!5!white,colframe=blue!75!black, breakable, enhanced]
Siguin $\set{(x_i, f_i)}_{i=0}^n$ el conjunt de punts de suport. Aleshores se satisfà:

\begin{itemize}
    \item $P_i(x) = f_i$
    \item $P_{i_0, i_1, \dots, i_k}(x) = \frac{(x - x_{i_0})P_{i_1, \dots, i_k}(x) - (x - x_{i_k})P_{i_0, \dots, i_{k-1}}(x)}{x_{i_k} - x_{i_0}}$
\end{itemize}

Podem organitzar els càlculs de l'algoritme de Neville en una taula:
\\ \\
\begin{tabular}{c|cccccc}
    & $k = 0$ & $k = 1$ & $k = 2$ & \dots \\
    \hline
    $x_0$ & $f_0 = P_0(x)$ \\
    & & $P_{0, 1}(x)$ & & \phantom{\vdots} \\
    $x_1$ & $f_1 = P_1(x)$ & & $P_{0, 1, 2}(x)$ \\
    & & $P_{1, 2}(x)$ & & \vdots \\
    $x_2$ & $f_2 = P_2(x)$ & & \vdots \\
    \vdots & \vdots & \vdots \\
\end{tabular}
\\ \\

Finalment, el polinomi interpolador de Lagrange que busquem serà $P_{0, 1, \dots, n}(x)$.

\end{tcolorbox}

\newpage

\section{Manual de software}

\subsection{Funció \texttt{newton}}
La funció \texttt{newton} té el prototipus següent:

\fn{newton}{int newton (double *x, double tolf, double tolx, int maxit, void (*fdf)(double,double*,double*,void*), void *prm, int ixrr);}

Aquesta funció troba el valor de $x$ que, dins la tolerància donada, troba el zero de una funció. Els arguments són:
\begin{itemize}
    \item \texttt{x}: Punter al valor inicial de $x$.
    \item \texttt{tolf}: Tolerància per a la imatge de $x$.
    \item \texttt{tolx}: Tolerància per a la diferència entre valors consecutius de $x$.
    \item \texttt{maxit}: Nombre màxim d'iteracions.
    \item \texttt{fdf}: Punter a la funció que calcula la funció i la seva derivada.
    \item \texttt{prm}: Punter a paràmetres addicionals a la funció \texttt{fdf}.
    \item \texttt{ixrr}: Índex de xarrera, indica si imprimeix valors sobre les iteracions per \texttt{stderr}.
\end{itemize}

\newpage

\subsection{Funció \texttt{neville}}
La funció \texttt{neville} té el prototipus següent:

\fn{neville}{double neville (int n, double xi[], double yi[], double x);}

Aquesta funció calcula l'interpolació de Newton per un conjunt de punts donats. Els arguments són:
\begin{itemize}
    \item \texttt{n}: Hi ha $n + 1$ punts.
    \item \texttt{xi}: Vector de valors de $x$.
    \item \texttt{yi}: Vector de valors de $y$.
    \item \texttt{x}: Valor de $x$ per el qual calcular l'interpolació.
\end{itemize}

\newpage

\subsection{Utilitat \texttt{detcsig}}
La utilitat \texttt{detcsig} reb unes condicions inicials i propaga la trajectòria, imprimint certs valors a cada pas. No reb cap valor per \texttt{stdin}.

\bash{./detcsig r0 v0 alf0 alff nalf nt pmax}

on \texttt{r0} és el valor inicial de $r$, \texttt{v0} és el valor inicial de $v$, \texttt{alf0} és el valor inicial de $\alpha$, \texttt{alff} és el valor final de $\alpha$ i \texttt{nalf} és el nombre de valors de $\alpha$ a considerar. Els valors de \texttt{nt} i de \texttt{pmax} serveixen per passar-los a \texttt{proptraj}.

Exemple d'ús:

\bashio{./detcsig 0.01701377045867686=r0 10.67764174790536=v0 0.95=alf0 1.25=alff 100=nalf 2=nt .05=pmax}{0.95 0.07049004991490235 7385.357327340641\\0.953030303030303 0.05298628767590897 7803.10624193974\\...\\1.246969696969697 0.05226682892839966 35245.90328074239\\1.25 0.05497146275142845 35339.34197112153}

En particular, si podem redirigir aquesta sortida al programa gnuplot per obtenir:

\begin{figure}[!h]
    \centering
    \includegraphics[width=0.6\textwidth]{detcsig_plot.png}
    \caption{Grafica de la funció \texttt{detcsig}}
    \label{fig:detcsig_plot}
\end{figure}

Observem que en aquest cas, el flyby es produeix en una alçada on la velocitat s'anula, és a dir, en una alçada aproximada de $6047.70$ km (ja que el valor al voltant de $30000$ no és el què ens interessa).

\newpage

\subsection{Utilitat \texttt{trobalf}}
La utilitat \texttt{trobalf} reb unes condicions inicials de forma iterada, i per a cada iteració, propaga la trajectòria i calcula el punt allà on la trajectòria passa a prop de la lluna i talla l'eix d'ordenades.

\bash{./trobalf tolu maxit pmax ixrr}

on \texttt{tolu} és la tolerància, \texttt{maxit} el nombre màxim d'iteracions, i \texttt{ixrr} el valor de \texttt{ixrr}. El valor de \texttt{pmax} serveixe per passar-lo a \texttt{proptraj}.

A més a més, per \texttt{stdin} reb quatre valors per a cada iteració: \texttt{r0, v0}, $\alpha_0$, \texttt{nt}, on \texttt{r0} és la distància inicial, \texttt{v0} és la velocitat inicial, $\alpha_0$ és l'angle inicial, i \texttt{nt} és el nombre de talls amb l'eix d'ordenades.

Exemple d'ús (seguint l'exemple de la utilitat anterior):

\bashio{echo 0.01701377045867686 10.67764174790536 0.9651515151515151 2 | ./trobalf 1e-10 10 .05 1}{newton: it=0 x=0.9651515151515151 f(x)=0.0007233922573357029 f'(x)=-3.387472267594433\\newton: it=0 x=0.9653650644296172 |dx|=0.000213549278102132\\newton: it=1 x=0.9653650644296172 f(x)=2.880323103687877e-06 f'(x)=-3.360536680934217\\newton: it=1 x=0.9653659215317343 |dx|=8.571021170578064e-07\\newton: it=2 x=0.9653659215317343 f(x)=4.648310961316571e-11 f'(x)=-3.360429056389705\\0.01701377045867686 10.67764174790536 0.9653659215317343 6046.557984050551}

Obtenim un resultat coherent que indica que l'angle que busquem és $\alpha = 0.9653659215317343$ per a aconseguir una velocitat horitzontal nul·la al creuar l'eix d'ordenades.

\newpage
Amb la utilitat donada \texttt{shtmn} (shoot the moon) podem obtenir un conjunt de punts, que, després de graficar-los, ens mostren la òrbita que busquem. En el nostre cas, si executem (noteu que hem utilitzat els valors obtinguts en executar \texttt{trobalf}):

\bash{echo 0.01701377045867686 10.67764174790536 0.9653659215317343 3 | ./shtmn .05 trajsex.txt}

\bash{gnuplot\\gnuplot> plot "trajsex.txt" u 2:3 w l}

Obtindrem el gràfic que volíem, amb els tres talls a l'eix horitzontal:
\begin{figure}[!h]
    \centering
    \includegraphics[width=0.8\textwidth]{trobalf_plot.png}
    \caption{Órbita de la nau espacial amb els valors trobats}
    \label{fig:trobalf}
\end{figure}

\newpage

\subsection{Utilitat \texttt{interp}}

La utilitat \texttt{interp} calcula el polinomi interpolador a partir d'uns punts de suport i avalua el polinomi en uns altres punts donats.

\bash{./interp n}

on \texttt{n} és el nombre de punts de suport. Per altra banda, reb per \texttt{stdin} $n+1$ punts de suport, i a continuació un nombre indeterminat de punts a avaluar.

Per exemple, suposem que tenim els punts de suport $(0, 1), (1, 3), (3, 2)$ i volem avaluar el polinomi en el punt $x = 2$. Aleshores, executem:

\bashio{echo 0 1 1 3 3 2 2 | ./interp n}{2 3.333333333333333}

Podem comprovar que això és correcte fent la taula de l'algoritme de Neville:

\begin{tabular}{c|ccc}
    & $k = 0$ & $k = 1$ & $k = 2$ \\
    \hline
    $x_0 = 0$ & $f_0 = P_0(2) = 1$ \\
    & & $P_{0,1}(2) = 5$ \\
    $x_1 = 1$ & $f_1 = P_1(2) = 3$ & & $P_{0,1,2}(2) = 10/3$ \\
    & & $P_{1,2}(2) = 5/2$ \\
    $x_2 = 3$ & $f_2 = P_2(2) = 2$ \\
\end{tabular}

Per tant el resultat és correcte.

\newpage

\section{Resolució del problema}

Ara, ens volem plantejar aconseguir una alçada de persileni de $6545$ km.

Per fer-ho, modificarem els valors de $v_0$ per a obtenir una alçada propera a $6545$ km, i quan la tinguem, aplicarem interpolació polinòmica amb l'algoritme de Neville utilitzant la utilitat \texttt{interp} per a interpolar la alçada que busquem i trobar $v_0$.

Provem alguns valors i creem la següent taula:

\begin{tabular}{c|c}
    $v_0$ & $h$ \\
    \hline
    $10.67764174790536$ & $6046.557984050551$ \\
    $10.67754174790536$ & $6201.149491899968$ \\
    $10.67744174790536$ & $6364.005404339322$ \\
    $10.67735174790536$ & $6518.393064334263$ \\
    $10.67734174790536$ & $6536.038703405597$ \\
    $10.67733174790536$ & $6553.786761598189$ \\
    $10.67732174790536$ & $6571.638426123778$ \\
    $10.67724174790536$ & $6718.331514385751$ \\
    $10.67714174790536$ & $6912.181634825405$ \\
\end{tabular}

Apliquem \texttt{interp} a la taula anterior per a interpolar la alçada que busquem i trobar $v_0$. Guardem la taula en un fitxer que anomenarem \texttt{trobalf.txt} i executem \texttt{interp} amb el següent comandament:

\bashio{(cat trobalf.txt | awk '{print \$2,\$1}' ; echo 6545) | ./interp 8}{6545 10.67733669148103}

Comprovem aquest resultant executant de nou la utilitat \texttt{trobalf} amb el valor interpolat $v_0$ per veure si obtenim l'alçada que busquem:

\bashio{echo 0.01701377045867686 10.67733669148103 0.9653659215317343 2 | ./trobalf 1e-10 10 .05 0}{0.01701377045867686 10.67733669148103 0.9731833801332418 6545.000000012091}

Efectivament, hem comprovat que, amb la velocitat inicial $v_0=10.67733669148103$, i amb l'angle inicial $\alpha_0 = 0.9653659215317343$, l'alçada que obtenim és $6545$ km.

\end{document}
