\documentclass[10pt]{beamer}
%\documentclass[ignorenoframetext,handout]{beamer}  % per stampare senza pause
\usetheme{Padova}
\usepackage{amsfonts}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsopn}
\usepackage{amsthm}
\usepackage{mathrsfs}
\usepackage{soul}
\usepackage{enumerate}

\usepackage{marvosym}
\usepackage{yhmath}
\usepackage{xcolor}
\usepackage[T1]{fontenc}
\usepackage[utf8]{inputenc}
\usepackage{subfigure}
\usepackage{multirow}
\usepackage{graphicx}
\usepackage{esvect}
\usepackage{booktabs}
\usepackage{adjustbox}
\usepackage{caption}
\usepackage{tikz}
\usepackage{wrapfig}
\usepackage{lipsum}
\usepackage{wasysym}

\captionsetup{labelformat=empty}

\newtheorem{lem}{Lemma}
\newtheorem{rem}{Remark}
\newtheorem{thm}{Teorema}
\newtheorem{proposition}{Proposition}
\theoremstyle{definition}
% \newtheorem{definition}{Definizione}

\definecolor{blu}{RGB}{0,155,135}

\newenvironment{system}%
	{\left\lbrace \begin{array}{@{} l @{} }}%
	{ \end{array}\right.}

\newcommand{\bbR}{\mathbb{R}}
\newcommand{\bbN}{\mathbb{N}}
\newcommand{\bbP}{\mathbb{P}}
\newcommand{\bbZ}{\mathbb{Z}}
\newcommand{\calX}{\mathcal{X}}
\newcommand{\calV}{\mathcal{V}}
\newcommand{\calR}{\mathcal{R}}

\newcommand*{\vimagea}{\vcenter{\hbox{\includegraphics[width=0.4\textwidth]{Esp01_Nodi Full.png}}}}
\newcommand*{\vimageb}{\vcenter{\hbox{\includegraphics[width=0.4\textwidth]{Esp01_Nodi compressi.png}}}}
\newcommand*{\vpointer}
{\vcenter{\hbox{\scalebox{2}{\Huge\pointer}}}}

\newcommand*{\vimagebb}{\vcenter{\hbox{\includegraphics[width=0.25\textwidth]{Esp01_Nodi compressi.png}}}}
\newcommand*{\vimageaa}{\vcenter{\hbox{\includegraphics[width=0.25\textwidth]{Esp01_Nodi Full.png}}}}

\newcommand{\de}{\,\mathrm{d}}
\newcommand{\dist}{\,\mathrm{dist}}
\newcommand{\dx}{\,\mathrm{d} x_1 \ldots \mathrm{d} x_n}

\newcommand{\px}{p(x_1, \ldots, x_n)}
\newcommand{\bs}{\boldsymbol}

\DeclareMathOperator{\ctg}{ctg}
\DeclareMathOperator{\cst}{cst}

\usepackage[tikz]{bclogo}
\presetkeys{bclogo}{
ombre=true,
epBord=2,
couleur = white,
couleurBord = rossoPantano,
arrondi = 0.1,
logo= 
}{}

\definecolor{azzurro}{RGB}{0,155,135}
\newcommand{\gel}[1]{\textcolor{azzurro}{#1}}


\title[{\small{Quadratura di elementi poligonali con un lato curvo in Python}}\hspace{2em}\insertframenumber/
\inserttotalframenumber]{Quadratura di elementi poligonali con un lato curvo in Python}

\author{GIOVANNI TRAVERSIN}
\institute{Universit\`a degli Studi di Padova}
\date{Laurea Triennale in Matematica \\ 17/09/2025}

\begin{document}

\frame{\titlepage}







\begin{frame}{L'obiettivo}
L'obiettivo di questo lavoro \`e fornire una fedele traduzione in {\sc Python} di routines {\sc Matlab} implementate in 
\begin{center}
{\em{Algebraic cubature on polygonal elements with a circular edge}}
\end{center}
da E. Artioli, A. Sommariva e M. Vianello.

{\vspace{0.3cm}}

Tali codici calcolano una formula di quadratura algebrica per funzioni aventi come dominio un poligono circolare, sfruttando tecniche quali
\begin{itemize}
\item il {\it linear blending di archi ellittici},
\item la {\it compressione della formula di quadratura}.
\end{itemize}
{\vspace{0.3cm}}

{\bf{Nota:}} Queste formule hanno impiego nei {\it{virtual elements}} (risoluzione di certi problemi alle derivate parziali).
\end{frame}%







\begin{frame}{I domini considerati}

Con {\it{elemento poligonale con un lato curvo}} considereremo un qualsiasi poligono convesso $\mathcal{P}$ a cui \`e stato sostituito uno dei lati con un arco di circonferenza. Tale arco pu\`o essere sia convesso che concavo.

\begin{figure}[h!]
    \centering
    \includegraphics[width=4cm]{PC1.png}
    \qquad
    \includegraphics[width=4.5cm]{PC2.PNG.png} 
    \caption{}%
    \label{fig:P_C_1_2 }%
\end{figure}
\end{frame}







\begin{frame}{Le formule di integrazione}

Nello specifico, i programmi costruiti restituiscono due formule di quadratura di tipo algebrico ({\it{integrano esattamente polinomi algebrici fino ad un certo grado}}), entrambe a nodi interni e pesi positivi.

\begin{itemize}
    \item La prima ha grado di precisione $n$ e cardinalit\`a dell'insieme dei nodi pari a $M_n \geq \frac{(n+1)(n+2)}{2}$:
    \[ \int \int_\Omega p(x,y) \ dxdy = \sum_{i=1}^{M_n} \ \omega_i \ p(x_i,y_i), \ \forall \ p \in \mathbb{P}_n ^2 \  .\]

    \item La seconda ha nuovamente grado di precisione $n$, cardinalit\`a $M^{(c)}_n=\frac{(n+1)(n+2)}{2}$ ed \`e una compressione della prima tramite l'algoritmo di Lawson-Hanson. 
\end{itemize}


\end{frame}



\begin{frame}{Compressione delle formule}

\centering
$\vimagea{}\ \ \ \ \ \vpointer \ \ \ \ \
\vimageb{}$


\end{frame}



\begin{frame}{Il fulcro della quadratura}

{\bf{Punto chiave}}: sfruttando la propriet\`a additiva dell'integrale, suddividiamo i poligoni circolari in sottodomini pi\`u semplici 
per definire progressivamente su ogni suddivisione una formula di quadratura del tipo richiesto.

{\vspace{0.5cm}}

Sui sottodomini disponiamo di formule di quadratura algebriche di grado prefissato avente nodi interni e pesi positivi e quindi, per l'additivit\'a dell'operatore di integrazione, abbiamo simili formule anche sull'elemento poligonale con un lato curvo.

\end{frame}



\begin{frame}{Elementi poligonali con lato curvo convesso}

Per elementi poligonali con lato curvo convesso, quindi, consideriamo:

{\vspace{0.3cm}}

\begin{itemize}

    
    \item {\bf{suddivisione}}: un poligono e un settore circolare; il poligono viene quindi triangolato in modo minimale;

    \item {\bf{quadratura}}: formule algebriche, a bassa cardinalit\`a e del grado richiesto, su ogni dominio della suddivisione:

{\vspace{0.3cm}}

    \begin{enumerate}
\item {\bf{sul triangolo}}: formule quasi minimali;

{\vspace{0.3cm}}

\item {\bf{sul settore circolare}}: formule di tipo {\it{linear blending di archi ellittici}}.
\end{enumerate}

\end{itemize}

\end{frame}












\begin{frame}{Elementi poligonali con lato curvo concavo}

Per elementi poligonali con lato curvo concavo, quindi, consideriamo:
{\vspace{0.3cm}}

{\bf{Suddivisione:}} due poligoni e alcuni quadrilateri circolari.

       \begin{figure}
\includegraphics[width=0.4\textwidth]{dhdhfh.PNG}
\end{figure}

\begin{itemize}
\item $\{P_1,..., P_{\ell_1} \}$, che contiene i vertici consecutivi raccolti a partire da $P_1$ esterni all'angolo $(\theta_1, \theta_\ell)$;
\item $\{ P_{\ell_2},..., P_\ell \}$, che contiene i vertici consecutivi raccolti a partire da $P_{\ell}$ esterni all'angolo $(\theta_1, \theta_\ell)$;
\item $\{ P_{\ell_1 + 1} ,..., P_{\ell_2 - 1} \}$ che contiene tutti gli eventuali vertici restanti.
\end{itemize}




\end{frame}






\begin{frame}{Elementi poligonali con lato curvo concavo}

{\bf{Formule}}: di tipo algebrico, a bassa cardinalit\`a e del grado richiesto su ogni dominio della suddivisione:

{\vspace{0.3cm}}

    \begin{enumerate}
\item {\bf{triangolo}}: formule quasi-minimali ({\it{cardinalit\`a prossima alla minima possibile per il grado di precisione fissato}});

{\vspace{0.3cm}}

\item {\bf{quadrilateri circolari}}: algebrico di tipo {\it{linear blending di archi ellittici}}.
\end{enumerate}

\end{frame}















\begin{frame}{Esperimenti numerici}
Di seguito facciamo degli esperimenti su un elemento poligonale con lato curvo di tipo concavo per mostrare che
\begin{itemize}
\item le formule non compresse dei due sistemi sono numericamente identiche,
\item le formule compresse hanno le propriet\`a richieste ma possono essere diverse viste le differenti routines che implementano l'algoritmo di Lawson-Hanson nei due sistemi di calcolo.
\end{itemize}

Completeremo con alcuni esperimenti numerici su una certa integranda, ottenendo risultati comparabili offerti dai codici {\sc{Matlab}} e {\sc{Python}}.

\end{frame}



\begin{frame}{Nodi e pesi delle formule}
Mettendo a confronto, tramite la norma infinito, i nodi e i pesi trovati in {\sc Matlab} con quelli raggiunti in {\sc Python}, si nota che:
\begin{itemize}
\item i nodi e i pesi trovati per le formule originali sono gli stessi in entrambi i linguaggi;
\item i nodi e i pesi trovati per le formule compresse possono non essere gli stessi.
\end{itemize}

{\vspace{0.3cm}}
\centering
\scalebox{0.78}{



\begin{tabular}{|c||c|c|c|c|}
    \hline
    ADE & $E_{nodes}$ & $E_{weights}$ & $E^{(c)}_{nodes}$ & 
    $E^{(c)}_{weights}$ \\
    \hline
    \hline
    2 & 5.0 $\cdot 10^{-16}$  & 4.7 $\cdot 10^{-16}$ & 4.4 $\cdot 10^{-16}$ & 4.9 $\cdot 10^{-16}$ \\
    \hline
    4 & 5.6 $\cdot 10^{-16}$ & 5.0 $\cdot 10^{-16}$ & 4.4 $\cdot 10^{-16}$ & 4.8 $\cdot 10^{-16}$ \\
    \hline
    6 & 5.6 $\cdot 10^{-16}$ & 5.0 $\cdot 10^{-16}$ & {\color{red}{3.3 $\cdot 10^{-1}$}} & {\color{red}{1.7 $\cdot 10^{-2}$}} \\
    \hline
    8 & 5.3 $\cdot 10^{-16}$ & 5.0 $\cdot 10^{-16}$ & {\color{red}{1.9 $\cdot 10^{-1}$}} & {\color{red}{5.0 $\cdot 10^{-3}$}} \\
    \hline
    10 & 5.6 $\cdot 10^{-16}$ & 5.0 $\cdot 10^{-16}$ & {\color{red}{3.0 $\cdot 10^{-1}$}} & {\color{red}{6.7 $\cdot 10^{-3}$}} \\
    \hline
\end{tabular}
}
\end{frame}




\begin{frame}{Formule non completamente equivalenti}

Nonostante le formule compresse siano differenti sono entrambe accettabili, in quanto soddisfano la condizione dei momenti:

{\vspace{0.3cm}}

\begin{center}
\scalebox{0.78}{

\begin{tabular}{|c||c|c|}
    \hline
    ADE & ${\tt{momerr}}_M$ & ${\tt{momerr}}_P$ \\
    \hline
    \hline
    2 & 7.4 $\cdot 10^{-18}$ & 5.9 $\cdot 10^{-18}$ \\
    \hline
    4 & 9.0 $\cdot 10^{-18}$ & 7.1 $\cdot 10^{-18}$ \\
    \hline
    6 & 9.9 $\cdot 10^{-18}$ & 6.8 $\cdot 10^{-18}$ \\
    \hline
    8 & 9.3 $\cdot 10^{-18}$ & 6.3 $\cdot 10^{-18}$ \\
    \hline
    10 & 1.4 $\cdot 10^{-17}$ & 5.5 $\cdot 10^{-18}$ \\
    \hline
\end{tabular}
}
\end{center}

{\vspace{0.3cm}}

Il parametro {\tt{momerr}} valuta l'errore sui momenti in norma 2 di una certa base polinomiale a grado $n$ e, qualora sia vicino alla precisione di macchina, indica che la formula \`e numericamente esatta. 


\end{frame}





\begin{frame}{Equivalenza degli integrali calcolati}

Calcoliamo l'integrale su tale elemento poligonale di tipo concavo $\Omega$ della funzione
\begin{equation}
\int_{\Omega} e^{-(x-0.2)^2-(y-0.2)^2} dx dy \approx 0.145174464220773263845.
\end{equation}

{\vspace{0.3cm}}
\begin{center}
\begin{tabular}{|c||c|c|c|c|}
    \hline
    ADE & ER\_M & ER\_P & ERC\_M & ERC\_P \\
    \hline
    \hline
    2 & 6.3 $\cdot {10}^{-5}$ & 6.3 $\cdot {10}^{-5}$ &  1.2 $\cdot 10^{-4}$ & 1.2 $\cdot {10}^{-4}$ \\
    \hline
    4 & 2.0 $\cdot {10}^{-7}$ & 2.0 $\cdot {10}^{-7}$ & 3.3 $\cdot {10}^{-7}$ & 3.3 $\cdot {10}^{-7}$ \\
    \hline
    6 & 4.7 $\cdot {10}^{-10}$ & 4.7 $\cdot {10}^{-10}$ & 1.0 $\cdot {10}^{-9}$ & 1.0 $\cdot {10}^{-9}$ \\
    \hline
    8 & 8.6 $\cdot {10}^{-13}$ & 8.6 $\cdot {10}^{-13}$ & 1.9 $\cdot {10}^{-12}$ & 1.9 $\cdot {10}^{-12}$ \\
    \hline
    10 & 2.5 $\cdot {10}^{-15}$ & 2.3 $\cdot {10}^{-15}$ & 9.6 $\cdot {10}^{-16}$ & 1.1 $\cdot {10}^{-15}$ \\
    \hline
\end{tabular}
\end{center}
{\vspace{0.3cm}}

La tabella mostra che i risultati sono numericamente equivalenti.

\end{frame}






\begin{frame}{Tempi di calcoli}

Concludiamo coi tempi di calcolo che indicano che i codici in {\sc{Matlab}} sono leggermente pi\`u rapidi di quelli in {\sc{Python}}.

\begin{center}
{\it{Formule originali}}
\end{center}

{\vspace{-0.1cm}}

\centering
\scalebox{0.9}{
\begin{tabular}{|c||c|c|}
    \hline
    ADE & {\sc{Matlab}}& {\sc{Phyton}}\\
    \hline
    \hline
    2 & 2.545 $\cdot {10}^{-3}$ & 5.410 $\cdot {10}^{-3}$ \\
    \hline
    4 & 1.546 $\cdot {10}^{-3}$ & 5.432 $\cdot {10}^{-3}$ \\
    \hline
    6 & 1.845 $\cdot {10}^{-3}$ & 6.905 $\cdot {10}^{-3}$ \\
    \hline
    8 & 1.991 $\cdot {10}^{-3}$ & 1.698 $\cdot {10}^{-2}$ \\
    \hline
    10 & 2.786 $\cdot {10}^{-3}$ & 2.169 $\cdot {10}^{-2}$ \\
    \hline
\end{tabular}
}

\begin{center}
{\it{Formule compresse}}
\end{center}

{\vspace{-0.1cm}}
\scalebox{0.9}{
\begin{tabular}{|c||c|c|}
    \hline
    ADE & {\sc{Matlab}}& {\sc{Phyton}}\\
    \hline
    \hline
    2 & 2.406 $\cdot {10}^{-3}$ & 5.375 $\cdot {10}^{-3}$ \\
    \hline
    4 & 2.337 $\cdot {10}^{-3}$ & 6.218 $\cdot {10}^{-3}$ \\
    \hline
    6 & 3.376 $\cdot {10}^{-3}$ & 8.459 $\cdot {10}^{-3}$ \\
    \hline
    8 & 4.514 $\cdot {10}^{-3}$ & 2.845 $\cdot {10}^{-2}$ \\
    \hline
    10 & 9.514 $\cdot {10}^{-3}$ & 5.807 $\cdot {10}^{-2}$ \\
    \hline
\end{tabular}
}
\end{frame}















\begin{frame}
  \centering
  \color{white}
  \Huge{Grazie per l'attenzione}
\end{frame}



\begin{frame} {Appendice: Quadratura su quadrilateri circolari}

Il processo di quadratura per un quadrilatero circolare $\mathcal{Q}$ \`e pi\`u complesso rispetto a quello per un poligono.
Innanzitutto, si riscrive $\mathcal{Q}$ in forma parametrica mediante la tecnica di arc blending:
\[
\mathcal{Q} = \{ (x,y) = \sigma (t,\theta) = tP(\theta) + (1 - t)Q(\theta) \ \ ,  \ (t,\theta) \in [0,1] \times [\alpha,\beta] \}.
\]


con \[
P(\theta) = A_1\cos (\theta ) + B_1\sin (\theta )+ C_1, \quad Q(\theta) = A_2\cos (\theta ) + B_2\sin (\theta )+ C_2 . 
\]

Definiti i vertici del quadrilatero circolare come $V_i = (\xi_i,\eta_i)$ per $i \in \{1,2,3,4\}$, con $\wideparen{V_1 V_4}$ come lato curvo e $\alpha$ e $\beta$ le rispettive coordinate angolari di $V_1$ e $V_4$, si definiscono 

\[s = \sin \left(  \frac{\beta \ - \  \alpha}{2} \right) \ \ \ \ \ \ \ \ \ \ \ \ \ \
\phi = \frac{\beta \ + \ \alpha}{2}.\]

\end{frame}







\begin{frame} {Appendice: Quadratura su quadrilateri circolari}

\`E ora possibile quindi definire ogni poligono circolare come funzione bilineare ponendo
\[A_1 = (r,0), \qquad B_1 = (0,r) ,\qquad C_1 = (x_0,y_0),\]
\[
 A_2 = \frac{\sin ( \phi )}{2s} (V_3 - V_2), \qquad
 B_2 = \frac{\cos ( \phi )}{2s} (V_3 - V_2), \qquad
 C_2 = \frac{1}{2} (V_3 + V_2).
\]

{\vspace{0.3cm}}


Da qui \`e immediata la formula di cubatura con le caratteristiche richieste:

\[
\int \int_\Omega p(x,y) \ dxdy = \int \int_{[0,1] \times [\alpha, \beta]} p(\sigma(t,\theta))[\pm \det(J\sigma(t,\theta))] \ dtd\theta \ .
\]
\end{frame}




\begin{frame} {Appendice: Corrispondenze con i codici}


\begin{itemize}
    \item {\bf Formule di quadratura}
    {\vspace{0.3cm}}
    \begin{itemize}
        \item {\bf Polygauss}: produce nodi e pesi per la formula di quadratura gaussiana sui triangoli che compongono le parti poligonali dell'elemento circolare.
        {\vspace{0.3cm}}
        \item {\bf Circtrap}: produce nodi e pesi per la formula di quadratura gaussiana sulle componenti di lato curvo dell'elemento circolare. La tecnica di {\it arc blending} \`e implementata direttamente all'interno della funzione, mentre il calcolo dei valori \`e affidato alla funzione {\tt gqellblend}.
    \end{itemize}
\end{itemize}






\end{frame}















\begin{frame} {Appendice: Corrispondenze con i codici}



\begin{itemize}
    \item {\bf Compressione}
    {\vspace{0.3cm}}
    \begin{itemize}
        \item {\bf Chebvand}: costruisce la matrice di Chebyshev-Vandermonde per il programma lineare da risolvere.
        {\vspace{0.3cm}}
        \item {\bf Comprexcub}: risolve il sistema lineare a partire da Chebvand.
    \end{itemize}
\end{itemize}

{\vspace{0.3cm}}

\begin{itemize}
    \item {\bf Visualizzazione}
    {\vspace{0.3cm}}
    \begin{itemize}
        \item {\bf Demo\_polygcirc}: permette di visualizzare graficamente i gli elementi circolari usati per l'integrazione, e di stampare i valori riportati all'interno delle tabelle, oltre alle specifiche sulle formule utilizzate (grado di precisione, numero di nodi, algoritmo di compressione, etc.)
    \end{itemize}
\end{itemize}


\end{frame}





\begin{frame} {Appendice: Note sulla traduzione}
Per mantenere in {\sc Python} uno stile quanto pi\`u simile possibile al codice sorgente, all'interno della traduzione sono state prese le seguenti accortezze:
\begin{itemize}
    \item tutti i valori definiti in {\sc{Matlab}} tramite il comando {\tt nargin} sono definiti a priori come input facoltativi durante la costruzione delle varie funzioni in {\sc Python};
    \item i valori non definiti ma non calcolati in {\sc Matlab} tramite il comando {\tt nargout} sono stati riproposti in modo facoltativo anche nella traduzione tramite un comando {\tt if} controllato da una {\tt flag} facoltativa in input;
    \item all'interno dei cicli {\tt for} che vagliano gli indici delle variabili, il valore analizzato di volta in volta \`e gi\`a diminuito di $1$ rispetto agli indici {\sc{Matlab}}, salvo rare eccezioni in cui tali indici hanno effettivo valore numerico.
\end{itemize}



\end{frame}


\begin{frame} {Appendice: Compressione delle formule}
Per un teorema di Tchakaloff, da formule a pesi positivi con grado di precisione $n$ e un numero di nodi maggiore della dimensione ${\tt{dim}}_n=\frac{(n+1)(n+2)}{2}$ dello spazio polinomiale ${\mathbb{P}}_n$, possiamo ricavare una formula con cardinalit\`a al pi\`u ${\tt{dim}}_n$, pesi positivi e grado di esattezza nuovamente $n$.

{\vspace{0.3cm}}

Tale compressione avviene mediante l'applicazione dell'algoritmo di Lawson-Hanson
\[
 V^t u = b \ ,\  \ u\geq 0 ,
\]
dove $V$ \'e la matrice di Chebyshev-Vandermonde sui nodi della formula iniziale.

{\vspace{0.3cm}}

\centering
$\vimageaa{}\ \ \ \ \ \vpointer \ \ \ \ \
\vimagebb{}$

\end{frame}

\end{document}