...
 
Commits (3)
......@@ -110,14 +110,14 @@
\end{equation*}
De plus, à l'aide des méthodes de {\it Jacobi}
et {\it Gauss-Seidel} (voir Série 2) essayez d'obtenir la figure suivante.
\begin{figure}[H]
\centering
\includegraphics[scale=0.6]{convergence.pdf}
\caption{La ligne bleue correspond à la méthode de Jacobi,
la ligne rouge correspond à la méthode de Gauss-Seidel
et la ligne violette correspond à la méthode SOR.
Les lignes noires pointillées représentent les estimations de la convergence des trois méthodes.}
\end{figure}
% \begin{figure}[H]
% \centering
% \includegraphics[scale=0.6]{convergence.pdf}
% \caption{La ligne bleue correspond à la méthode de Jacobi,
%la ligne rouge correspond à la méthode de Gauss-Seidel
%et la ligne violette correspond à la méthode SOR.
%Les lignes noires pointillées représentent les estimations de la convergence des trois méthodes.}
% \end{figure}
\end{enumerate}
......
AAltermatt
Anne-CatherineBerisV
CanosaAlvesL
CorviniV
FolL
GeimerA
GrosjeanL
HorvathR
MaizCalleP
Perrin-GanierN
PistoresiL
RivollatQ
RohrbasserL
StraubC
WahlN
......@@ -5,22 +5,32 @@ function A=FDLaplacian(n,d)
h=1/(n+1);
D=1/h^2*spdiags(ones(n,1)*[1 -2*d 1],[-1 0 1],n,n);
D=1/h^2*spdiags(ones(n,1)*[1 -2 1],[-1 0 1],n,n);
if d==1
A=D;
A=D;
elseif d==2
X=spdiags(ones(n,1),0,n,n);
I2=1/h^2*spdiags(ones(n^2,1)*[1 1],[-n n],n^2,n^2);
A=kron(X,D)+I2;
I = eye(n);
A=kron(I,D) + kron(D,I);
elseif d==3
X=spdiags(ones(n^2,1),0,n^2,n^2);
I2=1/h^2*spdiags(ones(n^3,1)*[1 1],[-n n],n^3,n^3);
I3=1/h^2*spdiags(ones(n^3,1)*[1 1],[-n^2 n^2],n^3,n^3);
A=kron(X,D)+I2+I3;
I = eye(n);
D1=kron(I,D); D2=kron(D,I); II=kron(I,I);
A=kron(D1,I) + kron(D2,I) + kron(II,D);
end
% if d==1
% A=D;
% elseif d==2
% X=spdiags(ones(n,1),0,n,n);
% I2=1/h^2*spdiags(ones(n^2,1)*[1 1],[-n n],n^2,n^2);
% A=kron(X,D)+I2;
% elseif d==3
% X=spdiags(ones(n^2,1),0,n^2,n^2);
% I2=1/h^2*spdiags(ones(n^3,1)*[1 1],[-n n],n^3,n^3);
% I3=1/h^2*spdiags(ones(n^3,1)*[1 1],[-n^2 n^2],n^3,n^3);
% A=kron(X,D)+I2+I3;
% end
% avec blkdiag
......
% Series 1 Test Suite
n=101; d=1;
T=FDLaplacianTrue(n-2,d);
A=FDLaplacien(n-2,d);
disp(normest(A-T,2));
disp(normest((n-1)^2 * A-T,2));
% x=linspace(0,1,n)'; f=sin(x(2:end-1));
% error=abs(A*f + f);
% disp(norm(error(2:end-1),2))
n=11; d=2;
T=FDLaplacianTrue(n-2,d);
A=FDLaplacien(n-2,d);
disp(normest(A-T,2));
disp(normest((n-1)^2 * A - T,2));
% x=linspace(0,1,n);
% y=linspace(0,1,n)';
% f=sin(y(2:end-1)) * sin(x(2:end-1));
% error=abs(A*f(:) + f(:));
% disp(norm(error(2:end-1),2));
n=7; d=3;
T=FDLaplacianTrue(n-2,d);
A=FDLaplacien(n-2,d);
disp(normest(A-T,2));
disp(normest((n-1)^2 * A - T,2));
\ No newline at end of file
function u=Heat2D(f,g,p)
function [u,F]=Heat2D(f,g,p)
% HEAT2D solves the stationary heat equation over a square
% u=Heat2D(f,g,p) solves the heat equation over a square of size one
% by one. The source function is given by f, the boundary condition is
% given by g, the number of grid points for each direction is a quarter
% of the length(g) and one chooses the iterative solver with p,
% i.e. p=0 uses Jacobi, p=1 uses Gauss-Seidel.
n = length(g); n = n/4; ind=1:n;
I = eye(n);
G = kron(g(ind),I(:,1)) + kron(g(ind+n),I(:,end)) + ...
kron(I(:,1),g(ind+2*n)) + kron(I(:,end),g(ind+3*n));
G = (n+1)^2 * G;
F = f - G;
if p==0
LU = (n+1)^2*spdiags(ones(n,1)*[1 0 1],[-1 0 1],n,n);
LU = kron(I,LU) + kron(LU,I);
D = 1/( 4 * (n+1)^2 );
elseif p==1
end
u0 = kron(zeros(n,1),zeros(n,1));
err= 1;
tol= 1e-8;
itermax = 1000;
iter = 1;
while err > tol && iter < itermax
u = D*LU*u0 - D*F;
err = norm(u-u0);
u0= u;
iter = iter+1;
end
\ No newline at end of file
function D=Laplacian1D(n)
% LAPLACIAN1D compute a finite difference Laplacian
% A=FDLaplacian(n,d) computes a finite difference Laplacian on the
% unit interval using n interior points: 0 .. h .. 2h .. .. nh .. 1
h=1/(n+1);
D=1/h^2*spdiags(ones(n,1)*[1 -2 1],[-1 0 1],n,n);
\ No newline at end of file
% test suite for serie 2
n = 10;
g = zeros(4*n,1);
f = -ones(n^2,1);
[u,F] = Heat2D(f,g,0);
u = reshape(u,n,n);
surf(u)
L = FDLaplacianTrue(n,2);
L*u(:) - F
\ No newline at end of file
......@@ -10,6 +10,7 @@
\addtolength{\voffset}{-2cm}
\usepackage{hyperref}
\usepackage{float}
%% commandes utiles %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newcommand{\R}{\mathbb{R}}
......@@ -33,7 +34,7 @@
\begin{enumerate}
\item (2 pts) Ecrire les méthodes itératives de {\it Jacobi} et {\it Gauss-Seidel}
sous les deux formes standard, de correction, de résidu et de différences.
sous les formes standard, de correction, de résidu et de différences.
\item (4 pts) Soit $\vec{B}\in\C^{n\times n}$ une matrice carrée dont le
rayon spectral $\rho(\vec{B})$ vérifie $\rho(\vec{B})<1$. Démontrer que
......@@ -74,7 +75,7 @@
En Matlab, la fonction \texttt{kron(A,B)} calcule le produit de Kronecker.
\begin{enumerate}
\item Construire le Laplacien discret en dimension deux en utilisant \texttt{kron}.
Montrer qu'il est le m\^eme du r\'esultat de la fonction \texttt{FDLaplacian(n,2)} d'exercice 5 de s\'erie 1.
Montrer qu'il est le m\^eme du r\'esultat de la fonction \texttt{FDLaplacian(n,2)} de l'exercice 5 de la s\'erie 1.
\item Soit $\left \{ \lambda_i \right \}$ l'ensemble des valeurs propres du Laplacien discret en dimension une.
D\'emontrer que $\left \{ \lambda_i + \lambda_j \right \}$ est l'ensemble des valeurs propres du Laplacien discret en dimension deux.
Quelle est la relation entre les vecteurs propres?
......
\documentclass[a4paper]{article}
\input{style/head.tex}
%-------------------------------
% TITLE VARIABLES (identify your work!)
%-------------------------------
\newcommand{\yourname}{} % replace YOURNAME with your name
\newcommand{\yournetid}{} % replace YOURNETID with your NetID
\newcommand{\youremail}{} % replace YOUREMAIL with your email
\newcommand{\assignmentnumber}{2} % replace X with assignment number
\usepackage{preamble}
\usepackage{pdfpages}
\theoremstyle{definition}
\newtheorem{ex}{Exercice}
\begin{document}
%-------------------------------
% TITLE SECTION (do not modify unless you really need to)
%-------------------------------
\input{style/header.tex}
%-------------------------------
% ASSIGNMENT CONTENT (add your responses)
%-------------------------------
\begin{ex}
Puisque $\rho(B) < 1$, $1 \notin \sigma(B)$, o\`u $\sigma(B)$ est le spectre de $B$.
Alors, $\det(I-B)\neq 0$, qui veut dire que $(I-B)$ est inversible.
To show the Neumann series, we first consider the partial sum $\sum_{j=0}^n B^j$
and multiply it by $(I-B)$ to get
\begin{equation*}
(I-B) \sum_{j=0}^n B^j = \sum_{j=0}^n B^j - B^{j+1} = B^0 - B^{n+1} = I - B^{n+1},
\end{equation*}
where we used the fact that the second term is a telescoping sum.
Now, we multiply both terms by $(I-B)^{-1}$ to obtain
\begin{equation*}
\sum_{j=0}^n B^j = (I-B)^{-1} ( I - B^{n+1} ) = (I-B)^{-1} - (I-B)^{-1} B^{n+1}.
\end{equation*}
By subtracting $(I-B)^{-1}$, we get
\begin{equation*}
\sum_{j=0}^n B^j - (I-B)^{-1} = - (I-B)^{-1} B^{n+1}.
\end{equation*}
For any arbitrary $\epsilon>0$ and $k$ large enough
we have that $\| B^k \| \leq ( \rho(B) + \epsilon )^k$.
Now, we write
\begin{equation*}
\| \sum_{j=0}^n B^j - (I-B)^{-1} \| = \| (I-B)^{-1} B^{n+1} \| \leq \| (I-B)^{-1} \| \| B^{n+1} \| \leq \| (I-B)^{-1} \| ( \rho(B) + \epsilon )^{n+1}.
\end{equation*}
Since $\rho(B)<1$, we can choose $\epsilon$ sufficiently small such
that $( \rho(B) + \epsilon ) < 1$ holds.\footnote{Notice also that $\rho(B)<1$, we can find a matrix norm $\| \cdot \|$
such that $\| B \| < 1$ (see Problem \ref{prob:2.3}), therefore we have $\|B^k\|\leq\|B\|^k$
and $\lim_{k\rightarrow\infty}\|B\|^k = 0$.}
Taking the limit for $n \rightarrow \infty$ we get
\begin{equation*}
\lim_{n \rightarrow \infty} \| \sum_{j=0}^n B^j - (I-B)^{-1} \| \leq
\lim_{n \rightarrow \infty} \| (I-B)^{-1} \| ( \rho(B) + \epsilon )^{n+1} = 0.
\end{equation*}
This means that
\begin{equation*}
\lim_{n \rightarrow \infty} \sum_{j=0}^n B^j = (I-B)^{-1}.
\end{equation*}
\end{ex}
\begin{ex}
J: M=D, N=A-D.
GS: M=D+L, N=-U.
For a splitting $A = M - N$ we have that
\begin{eqnarray*}
\vec{u}_{k+1} &=& M^{-1}N \vec{u}_k + M^{-1}\vec{f} \qquad \quad \; \; \: \text{(standard form)}, \nonumber \\
\vec{u}_{k+1} &=& \vec{u}_k + M^{-1}\vec{r}_k \qquad \qquad \qquad \; \, \text{(correction form)}, \nonumber \\
\vec{e}_{k+1} &=& M^{-1}N \vec{e}_k \qquad \qquad \qquad \qquad \text{(error form)}, \nonumber \\
\vec{d}_{k+1} &=& M^{-1}N \vec{d}_{k} \qquad \qquad \qquad \qquad \text{(difference form)}, \nonumber \\
\vec{r}_{k+1} &=& (I-AM^{-1})\vec{r}_{k} = N M^{-1} \vec{r}_{k} \; \text{(residual form)}. \nonumber \\
\end{eqnarray*}
Now, it is sufficient to replace into these formulas the splitting matrices $M$ and $N$ corresponding to the methods
of Jacobi, Gauss-Seidel, and SOR. For the method of Jacobi, we have
\begin{eqnarray*}
\vec{u}_{k+1} &=& -D^{-1}(L+U) \vec{u}_k + D^{-1}\vec{f} , \nonumber \\
\vec{u}_{k+1} &=& \vec{u}_k + D^{-1}\vec{r}_k , \nonumber \\
\vec{e}_{k+1} &=& -D^{-1}(L+U) \vec{e}_k , \nonumber \\
\vec{d}_{k+1} &=& -D^{-1}(L+U) \vec{d}_{k} , \nonumber \\
\vec{r}_{k+1} &=& (I-AD^{-1})\vec{r}_{k} = -(L+U)D^{-1} \vec{r}_{k} . \nonumber \\
\end{eqnarray*}
For the method of Gauss-Seidel, we have
\begin{eqnarray*}
\vec{u}_{k+1} &=& -(D+L)^{-1}U \vec{u}_k + (D+L)^{-1}\vec{f} , \nonumber \\
\vec{u}_{k+1} &=& \vec{u}_k + (D+L)^{-1}\vec{r}_k , \nonumber \\
\vec{e}_{k+1} &=& -(D+L)^{-1}U \vec{e}_k , \nonumber \\
\vec{d}_{k+1} &=& -(D+L)^{-1}U \vec{d}_{k} , \nonumber \\
\vec{r}_{k+1} &=& (I-A(D+L)^{-1})\vec{r}_{k} = -U(D+L)^{-1} \vec{r}_{k} . \nonumber \\
\end{eqnarray*}
Finally, for the SOR method, we have
\begin{eqnarray*}
\vec{u}_{k+1} &=& (\frac{1}{\omega}D+L)^{-1}(-U+(\frac{1}{\omega}-1)D) \vec{u}_k + (\frac{1}{\omega}D+L)^{-1}\vec{f}, \nonumber \\
\vec{u}_{k+1} &=& \vec{u}_k + (\frac{1}{\omega}D+L)^{-1}\vec{r}_k , \nonumber \\
\vec{e}_{k+1} &=& (\frac{1}{\omega}D+L)^{-1}(-U+(\frac{1}{\omega}-1)D) \vec{e}_k , \nonumber \\
\vec{d}_{k+1} &=& (\frac{1}{\omega}D+L)^{-1}(-U+(\frac{1}{\omega}-1)D) \vec{d}_{k} , \nonumber \\
\vec{r}_{k+1} &=& (I-A(\frac{1}{\omega}D+L)^{-1})\vec{r}_{k} = (-U+(\frac{1}{\omega}-1)D)(\frac{1}{\omega}D+L)^{-1} \vec{r}_{k} , \nonumber \\
\end{eqnarray*}
that are equivalent to
\begin{eqnarray*}
\vec{u}_{k+1} &=& (D+\omega L)^{-1}(-\omega U+(1-\omega)D) \vec{u}_k + \omega(D+\omega L)^{-1}\vec{f}, \nonumber \\
\vec{u}_{k+1} &=& \vec{u}_k + \omega (D+\omega L)^{-1}\vec{r}_k , \nonumber \\
\vec{e}_{k+1} &=& (D+\omega L)^{-1}(-\omega U+(1-\omega)D) \vec{e}_k , \nonumber \\
\vec{d}_{k+1} &=& (D+\omega L)^{-1}(-\omega U+(1-\omega)D) \vec{d}_{k} , \nonumber \\
\vec{r}_{k+1} &=& (I-\omega A(D+\omega L)^{-1})\vec{r}_{k} = (-\omega U+(1-\omega)D)(D+\omega L)^{-1} \vec{r}_{k} . \nonumber \\
\end{eqnarray*}
\end{ex}
\begin{ex}
La convergence arrive si $\norm{M^{-1}N}<1$.
Toutes les normes de matrices sont \'equivalentes.
\begin{align*}
\norm{M^{-1}N}_\infty & = \max_i \sum_{j \neq i} \frac{\abs{a_{ij}}}{a_{ii}} < 1
\end{align*}
\end{ex}
\newcommand{\G}{\begin{bmatrix} 5/2 & -1 \\ 1 & 0 \end{bmatrix}}
\begin{ex}
\begin{align*}
v_1 & = \G \begin{bmatrix} 1 \\ 2 \end{bmatrix} \\
& = \begin{bmatrix} 1/2 \\ 1 \end{bmatrix} \\
v_2 & = \G \begin{bmatrix} 1/2 \\ 1 \end{bmatrix} \\
& = \begin{bmatrix} 1/4 \\ 1/2 \end{bmatrix} \\
\det \left (\G - \lambda I \right ) & = \lambda^2 - \frac{5}{2}\lambda + 1 = 0 \\
\implies \lambda & = \frac{5/2 \pm \sqrt{25/4 - 16/4}}{2} \\
& = \frac{5/2 \pm 3/2}{2} \\
& = \begin{cases} 2 \\ 1/2 \end{cases} \\
\left (\G - 2I \right ) v = 0 & \implies v = \alpha \begin{bmatrix} 2 \\ 1 \end{bmatrix} \\
\left (\G - 1/2 I \right ) v = 0 & \implies v = \alpha \begin{bmatrix} 1/2 \\ 1 \end{bmatrix}
\end{align*}
\end{ex}
%------------------------------------------------
\end{document}
% Master Preamble
\NeedsTeXFormat{LaTeX2e}
\ProvidesPackage{preamble}
\RequirePackage{amsmath}
\RequirePackage{amssymb}
\RequirePackage{amsthm}
\RequirePackage{graphicx}
\RequirePackage{hyperref}
\RequirePackage{subcaption}
\newcommand{\dxdy}[2]{\frac{d #1}{d #2}}
\newcommand{\dxdyk}[3]{\frac{d^{#3} #1}{d {#2}^{#3}}}
\newcommand{\pdxdy}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\liminfty}[1]{\lim_{#1 \to \infty}}
\newcommand{\limab}[2]{\lim_{#1 \to #2}}
\newcommand{\abs}[1]{\left \vert #1 \right \vert}
\newcommand{\norm}[1]{\left \Vert #1 \right \Vert}
\newcommand{\order}[1]{\mathcal{O} \left ( #1 \right )}
\newcommand{\set}[1]{\left \{ #1 \right \}}
\newcommand{\Set}[2]{\left \{ #1 \ \middle \vert \ #2 \right \}}
\newcommand{\vmat}[1]{\begin{vmatrix} #1 \end{vmatrix}}
\DeclareMathOperator{\sign}{sign}
\newcommand{\bbn}{\mathbb{N}}
\newcommand{\bbz}{\mathbb{Z}}
\newcommand{\bbq}{\mathbb{Q}}
\newcommand{\bbr}{\mathbb{R}}
\newcommand{\bbc}{\mathbb{C}}
\newcommand{\bbf}{\mathbb{F}}
\newtheorem{lemma}{Lemma}
\newtheorem{thm}{Theorem}
\newtheorem{cor}{Corollary}
\newtheorem{prop}{Proposition}
\newtheorem{conj}{Conjecture}
\newtheorem{defn}{Definition}
\endinput
\ No newline at end of file
\addtolength{\hoffset}{-2.25cm}
\addtolength{\textwidth}{4.5cm}
\addtolength{\voffset}{-3.25cm}
\addtolength{\textheight}{5cm}
\setlength{\parskip}{0pt}
\setlength{\parindent}{0in}
\usepackage[square,sort,comma,numbers]{natbib}
\usepackage{blindtext} % Package to generate dummy text
\usepackage{charter} % Use the Charter font
\usepackage[utf8]{inputenc} % Use UTF-8 encoding
\usepackage{microtype} % Slightly tweak font spacing for aesthetics
\usepackage{amsthm, amsmath, amssymb} % Mathematical typesetting
\usepackage{float} % Improved interface for floating objects
\usepackage{hyperref} % For hyperlinks in the PDF
\usepackage{graphicx, multicol} % Enhanced support for graphics
\usepackage{xcolor} % Driver-independent color extensions
\usepackage{pseudocode} % Environment for specifying algorithms in a natural way
\usepackage[yyyymmdd]{datetime} % Uses YEAR-MONTH-DAY format for dates
\usepackage{fancyhdr} % Headers and footers
\pagestyle{fancy} % All pages have headers and footers
\fancyhead{}\renewcommand{\headrulewidth}{0pt} % Blank out the default header
\fancyfoot[L]{} % Custom footer text
\fancyfoot[C]{} % Custom footer text
\fancyfoot[R]{\thepage} % Custom footer text
\newcommand{\note}[1]{\marginpar{\scriptsize \textcolor{red}{#1}}} % Enables comments in red on margin
%----------------------------------------------------------------------------------------
\fancyhead[C]{}
\hrule \medskip
\begin{minipage}{0.295\textwidth}
\raggedright
\footnotesize
\yourname \hfill\\
\yournetid \hfill\\
\youremail
\end{minipage}
\begin{minipage}{0.4\textwidth}
\centering
\large
S\'{e}rie \assignmentnumber\\
\end{minipage}
\begin{minipage}{0.295\textwidth}
\raggedleft
\hfill\\
\end{minipage}
\medskip\hrule
\bigskip
% Source: ss17_wissschreib (Eva)
\lstset{
basicstyle=\ttfamily\scriptsize\mdseries,
keywordstyle=\bfseries\color[rgb]{0.171875, 0.242188, 0.3125},
identifierstyle=,
commentstyle=\color[rgb]{0.257813, 0.15625, 0},
stringstyle=\itshape\color[rgb]{0.0195313, 0.195313, 0.0117188},
numbers=left,
numberstyle=\tiny,
stepnumber=1,
breaklines=true,
frame=none,
showstringspaces=false,
tabsize=4,
backgroundcolor=\color[rgb]{0.98,0.98,0.98},
captionpos=b,
float=htbp,
language=Python,
xleftmargin=15pt,
xrightmargin=15pt
}
%(deutsche) Sonderzeichen
\lstset{literate=%
{Ä}{{\"A}}1
{Ö}{{\"O}}1
{Ü}{{\"U}}1
{ä}{{\"a}}1
{ö}{{\"o}}1
{ü}{{\"u}}1
{ß}{{\ss}}1
}
%Verwendung im Text:
%-> \begin{lstlisting}[language=Python,firstnumber=27] ... \end{lstlisting}
%-> \begin{lstlisting}[language=Python,numbers=none] ... \end{lstlisting}
%-> \lstinline[language=JAVA]{...}
\ No newline at end of file
\documentclass[12pt,a4paper]{article}
\usepackage{amsmath,amssymb,amsthm,epsfig,amscd,verbatim}
\usepackage[french]{babel}
\usepackage{fancyvrb}
%\usepackage[T1]{fontenc}
\usepackage[utf8]{inputenc}
\addtolength{\textwidth}{3cm}
\addtolength{\textheight}{4cm}
\addtolength{\hoffset}{-2cm}
\addtolength{\voffset}{-2cm}
\usepackage{hyperref}
\usepackage{float}
%% commandes utiles %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newcommand{\R}{\mathbb{R}}
\newcommand{\C}{\mathbb{C}}
\renewcommand{\vec}[1]{\mathbf{#1}}
\newcommand{\tr}{\mbox{tr}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}
\noindent UNIVERSIT\'E DE GEN\`EVE \hfill Section de Math\'ematiques \\
\noindent Facult\'e des Sciences \hfill \`A rendre le mercredi 1 avril \\[-3mm] \hrule
\bigskip
\begin{center}
\textbf{\Large{Méthodes Itératives}} \\[1mm]\bigskip
\textbf{Série 3}
\end{center}
\bigskip
\begin{enumerate}
\item \textbf{M\'ethodes it\'eratives stationnaires}
\begin{enumerate}
\item (3 pts) Considérer le Laplacien discret $A$ de $-\Delta$ en dimension 2.
Démontrer que la méthode de Jacobi et la méthode de Richardson (avec le paramètre optimal $\alpha_{opt}$) ont le même comportement de convergence.
\textit{Indice:} Pour trouver $\alpha_{opt}$, rappeler exercice 5(b) de la s\'erie 2 et que les valeurs propres d'une matrice tridiagonale et Toeplitz de taille $n \times n$ sont
\begin{equation*}
\lambda_n = a - 2 \sqrt{bc} \cos \left ( \frac{k \pi}{n+1} \right ), \ k = 1,...,n
\end{equation*}
o\`u $a$, $b$ et $c$ sont les valeurs sur les diagonales (primaire, supérieure et inférieure, respectivement).
\item (3 pts) Considérer une matrice $A \in \R^{2 \times 2}$ et le splitting $A = D + L + U$ avec $D$ inversible.
Démontrer que si la méthode de Jacobi converge, alors la méthode de Gauss-Seidel converge aussi et plus rapidement.
\item (4 pts) Considérer le Laplacien discret $A\in \R^{n^2\times n^2}$ de $-\Delta$ en dimension 2
et soit $h=\frac{1}{n+1}$ la taille de la grille.
Dénoter avec $\rho(G_{\rm J};h)$ et $\rho(G_{\rm GS};h)$ les rayons spectraux
des méthodes de Jacobi et Gauss-Seidel en fonction de $h$.
Démontrer que :
\begin{itemize}%\itemsep0em
\item $\rho(G_{\rm J};h) = 1 - \frac{\pi^2}{2} h^2 + O(h^4)$.
\item $\rho(G_{\rm GS};h) = 1 - \pi^2 h^2 + O(h^4)$.
\end{itemize}
\end{enumerate}
\item \textbf{M\'ethode de descente maximale}
\begin{enumerate}
\item (2 pts) Consid\'erer la fonction $F: \R^n \rightarrow \R$ d\'efinit par $F(u) = \frac{1}{2} u^\top A u - u^\top f$ pour une $A \in \R^{n \times n}$ sym\'etrique et $f \in \R^n$.
D\'emontrer que $\nabla F(u) = Au - f$.
\item (4 pts) Considérer une fonction $f : \R^n \rightarrow \R$ et on dénote avec $\nabla f(\vec{v})$ le
gradient de $f$ en $\vec{v} \in \R^n$ obtenu au moyen du produit scalaire habituel pour $\R^n$.
Démontrer que $\nabla f(\vec{v})$ est orthogonal à la ligne de niveau de $f$ en $\vec{v}$ et que
$-\nabla f(\vec{v})$ est la direction de descente maximale de $f$ en $\vec{v}$.
\item (4 pts) Implémenter une fonction Matlab qui résout un système linéaire $A \vec{u} = \vec{f}$ à l'aide de la méthode
du gradient.
(On peut utiliser le code \texttt{SteepestDescent} dans le polycopi\'e.)
R\'esoudre un syst\`eme $A \vec{u} = \vec{f}$ o\`u $A$ est le Laplacien en dimension deux et utiliser les r\'esultats \`a estimer le valeur de $\frac{\kappa(A) - 1}{\kappa(A) + 1}$.
\end{enumerate}
\end{enumerate}
$\;$\\
\\
{\bf \'Evaluation:}
\begin{itemize}
\item[$\bullet$] Les exercices
\item[$\bullet$] Un examen oral durant la session d'examens sur le cours.
\end{itemize}
La note finale est de : $30 \%$ exercices et $70\%$ examen oral.\\
Les exercices sont notés $15\%$ pour la présentation de séries et
$15\%$ pour la présentation au tableau. Il faut passer au tableau au
moins 2 fois dans le semestre.
\end{document}
function [u,uk]=SteepestDescent(A,f,u0,tol,maxiter)
% STEEPESTDESCENT solves Au=f using the steepest descent method
% [u,uk]=SteepestDescent(A,f,u0,tol,maxiter); solves Au=f using steepest
% descent starting at the initial guess u0 up to a tolerance tol using
% at most maxiter iterations. A has to be symmetric positive definite.
% uk contains all the iterates.
if nargin<5, maxiter = 100; end
if nargin<4, tol =0.0001; end
r=f-A*u0;
uk(:,1)=u0;
k=0;
while norm(r)/norm(f)>tol && k<maxiter
k=k+1;
al=r'*r/(r'*A*r);
uk(:,k+1)=uk(:,k)+al*r;
r=f-A*uk(:,k+1);
end
u=uk(:,k+1);
\ No newline at end of file
% Test suite for Serie 3
n = 10;
A = rand(n); A = -(A+A')/2;
A = diag(-sum(A,2)+0.1) + A;
u = 2*rand(n,1)-1; f = A*u;
[~,vk]= SteepestDescent(A,f,rand(n,1));
[~,vcol] = size(vk);
error = zeros(1,vcol);
for i=1:vcol
error(i) = norm(vk(:,i) - u);
end
semilogy(1:vcol,error)
%% Laplacian convergence
n = 7;
A = FDLaplacianTrue(n,2); A = -A;
x = linspace(0,1,n);
u = cos(x') * sin(x); u = u(:);
f = A*u;
C = condest(A); C = (C-1)/(C+1);
[v,vk] = SteepestDescent(A,f,ones(n^2,1));
[vrow,vcol] = size(vk);
error = zeros(vcol,1);
for i=1:vcol
w = vk(:,i) - u;
error(i) = sqrt(w'*A*w);
end
P = fit((1:vcol)',log(error),'poly1');
plot(1:vcol,log(error),'.-',1:vcol,P(1:vcol),'--')
xlabel('Iteration')
ylabel('log(Error)')
Cest = exp(P.p1); %Cest = (1+Cest)/(1-Cest);
\ No newline at end of file
\documentclass[a4paper]{article}
\input{style/head.tex}
%-------------------------------
% TITLE VARIABLES (identify your work!)
%-------------------------------
\newcommand{\yourname}{} % replace YOURNAME with your name
\newcommand{\yournetid}{} % replace YOURNETID with your NetID
\newcommand{\youremail}{} % replace YOUREMAIL with your email
\newcommand{\assignmentnumber}{3} % replace X with assignment number
\usepackage{preamble}
\usepackage{pdfpages}
\usepackage{enumitem}
\theoremstyle{definition}
\newtheorem{ex}{Exercice}
\begin{document}
%-------------------------------
% TITLE SECTION (do not modify unless you really need to)
%-------------------------------
\input{style/header.tex}
%-------------------------------
% ASSIGNMENT CONTENT (add your responses)
%-------------------------------
\begin{ex}
\begin{enumerate}[label={\alph*)}]
\item Noter que les valeurs propres de Laplacien en dimension deux sont:
\begin{equation*}
\lambda_{i,j} = -\frac{2}{h^2} \left ( 2 + \cos \left ( \frac{i \pi}{n+1} \right ) + \cos \left ( \frac{j \pi}{n+1} \right ) \right ).
\end{equation*}
Le plus grand et le plus petit sont $\lambda_{1,1}$ et $\lambda_{n,n}$, respectivement.
Rappeler que $\cos(\pi - \theta) = -\cos(\theta)$, alors $\lambda_{1,1} + \lambda_{n,n} = -8/h^2$.
Donc, $\alpha_{opt} = -h^2/4$.
Consid\'erer le splitting de Laplacien $A=D + L + U$.
Dans ce cas, $D^{-1} = \alpha_{opt} I$.
Alors, l’itération de Jacobi est
\begin{align}
u_{k+1} & = -D^{-1} (L+U) u_k + D^{-1} f \\
& = -D^{-1} (A - D) u_k + \alpha_{opt} f \\
& = (I - \alpha_{opt} A) a_k + \alpha_{opt} f
\end{align}
qui est exactement l'it\'eration de Richardson.
\item D'abord, \'ecrire $A=\begin{bmatrix}
a & b \\ c & d
\end{bmatrix}$.
Pour Jacobi:
\begin{align}
M^{-1} N & = -D^{-1}(U+L) = \begin{bmatrix}
0 & b/a \\ c/d & 0
\end{bmatrix} \\
\implies \rho(M^{-1}N) & = \abs{\sqrt{\frac{bc}{ad}}}.
\end{align}
On suppose que Jacobi converge, donc $\abs{bc/ad}<1$.
Pour GS:
\begin{align*}
M^{-1} N & = -(D+L)^{-1} U = \begin{bmatrix}
0 & b/a \\ 0 & -bc/ad
\end{bmatrix} \\
\implies \rho(M^{-1}N) & = \abs{\frac{bc}{ad}} < 1.
\end{align*}
\item
\end{enumerate}
\end{ex}
\begin{ex}
\begin{enumerate}[label={\alph*)}]
\item Consid\'erer la d\'eriv\'ee partielle de $F(u)$ en $u_i$:
\begin{align*}
\pdxdy{F(u)}{u_i} & = \frac{1}{2} \pdxdy{}{u_i} \begin{bmatrix}
u_1 & \dots & u_n
\end{bmatrix} A \begin{bmatrix}
u_1 \\ \vdots \\ u_n
\end{bmatrix} - \pdxdy{}{u_i} \begin{bmatrix}
u_1 & \dots & u_n
\end{bmatrix} f \\
& = \frac{1}{2} e_i^\top A u + \frac{1}{2} u A e_i^\top - e_i^\top f \\
& = \frac{1}{2} e_i^\top A u + \frac{1}{2} e_i^\top A^\top u - e_i^\top f \\
& = e_i^\top \left( \frac{1}{2} A u + \frac{1}{2} A^\top u - f \right) \\
& = e_i^\top \left( A u - f \right),
\end{align*}
o\'u $e_i$ est le $i$--i\`eme colonne de la matrice idendit\'e.
Quand on construit $\nabla F(u)$ avec $\pdxdy{F(u)}{u_i}$ on a $\nabla F(u) = I (Au - f) = Au - f$.
\item Soit $\vec{w}$ le vecteur tangent \`a la ligne de niveau de $f(\vec{v})$.
La d\'eriv\'ee directionnelle de $f$ dans la direction de $\vec{w}$ est $\nabla f(\vec{v}) \cdot \vec{w}$.
Dans cette direction le valeur de $f$ ne change pas, alors $\nabla f(\vec{v}) \cdot \vec{w} = 0$.
(Une d\'emonstration irait bien)
Soit $\vec{u}$ un vecteur unitaire dans une direction arbitraire qui fait un angle avec $\nabla f(\vec{v})$ de $\theta$.
Le produit scalaire entre $\vec{u}$ et $\nabla f(\vec{v})$ est:
\begin{align*}
\langle \nabla f(\vec{v}), \vec{u} \rangle & = \norm{\nabla f(\vec{v})} \norm{\vec{u}} \cos \theta \\
& = \norm{\nabla f(\vec{v})} \cos \theta
\end{align*}
qui est minimis\'e par $\theta = -\pi$.
Alors, $\vec{u}$ point dans la direction opposite de $\nabla f(\vec{v})$.
\end{enumerate}
\end{ex}
%------------------------------------------------
\end{document}
% Master Preamble
\NeedsTeXFormat{LaTeX2e}
\ProvidesPackage{preamble}
\RequirePackage{amsmath}
\RequirePackage{amssymb}
\RequirePackage{amsthm}
\RequirePackage{graphicx}
\RequirePackage{hyperref}
\RequirePackage{subcaption}
\newcommand{\dxdy}[2]{\frac{d #1}{d #2}}
\newcommand{\dxdyk}[3]{\frac{d^{#3} #1}{d {#2}^{#3}}}
\newcommand{\pdxdy}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\liminfty}[1]{\lim_{#1 \to \infty}}
\newcommand{\limab}[2]{\lim_{#1 \to #2}}
\newcommand{\abs}[1]{\left \vert #1 \right \vert}
\newcommand{\norm}[1]{\left \Vert #1 \right \Vert}
\newcommand{\order}[1]{\mathcal{O} \left ( #1 \right )}
\newcommand{\set}[1]{\left \{ #1 \right \}}
\newcommand{\Set}[2]{\left \{ #1 \ \middle \vert \ #2 \right \}}
\newcommand{\vmat}[1]{\begin{vmatrix} #1 \end{vmatrix}}
\DeclareMathOperator{\sign}{sign}
\newcommand{\bbn}{\mathbb{N}}
\newcommand{\bbz}{\mathbb{Z}}
\newcommand{\bbq}{\mathbb{Q}}
\newcommand{\bbr}{\mathbb{R}}
\newcommand{\bbc}{\mathbb{C}}
\newcommand{\bbf}{\mathbb{F}}
\newtheorem{lemma}{Lemma}
\newtheorem{thm}{Theorem}
\newtheorem{cor}{Corollary}
\newtheorem{prop}{Proposition}
\newtheorem{conj}{Conjecture}
\newtheorem{defn}{Definition}
\endinput
\ No newline at end of file
\addtolength{\hoffset}{-2.25cm}
\addtolength{\textwidth}{4.5cm}
\addtolength{\voffset}{-3.25cm}
\addtolength{\textheight}{5cm}
\setlength{\parskip}{0pt}
\setlength{\parindent}{0in}
\usepackage[square,sort,comma,numbers]{natbib}
\usepackage{blindtext} % Package to generate dummy text
\usepackage{charter} % Use the Charter font
\usepackage[utf8]{inputenc} % Use UTF-8 encoding
\usepackage{microtype} % Slightly tweak font spacing for aesthetics
\usepackage{amsthm, amsmath, amssymb} % Mathematical typesetting
\usepackage{float} % Improved interface for floating objects
\usepackage{hyperref} % For hyperlinks in the PDF
\usepackage{graphicx, multicol} % Enhanced support for graphics
\usepackage{xcolor} % Driver-independent color extensions
\usepackage{pseudocode} % Environment for specifying algorithms in a natural way
\usepackage[yyyymmdd]{datetime} % Uses YEAR-MONTH-DAY format for dates
\usepackage{fancyhdr} % Headers and footers
\pagestyle{fancy} % All pages have headers and footers
\fancyhead{}\renewcommand{\headrulewidth}{0pt} % Blank out the default header
\fancyfoot[L]{} % Custom footer text
\fancyfoot[C]{} % Custom footer text
\fancyfoot[R]{\thepage} % Custom footer text
\newcommand{\note}[1]{\marginpar{\scriptsize \textcolor{red}{#1}}} % Enables comments in red on margin
%----------------------------------------------------------------------------------------
\fancyhead[C]{}
\hrule \medskip
\begin{minipage}{0.295\textwidth}
\raggedright
\footnotesize
\yourname \hfill\\
\yournetid \hfill\\
\youremail
\end{minipage}
\begin{minipage}{0.4\textwidth}
\centering
\large
S\'{e}rie \assignmentnumber\\
\end{minipage}
\begin{minipage}{0.295\textwidth}
\raggedleft
\hfill\\
\end{minipage}
\medskip\hrule
\bigskip
% Source: ss17_wissschreib (Eva)
\lstset{
basicstyle=\ttfamily\scriptsize\mdseries,
keywordstyle=\bfseries\color[rgb]{0.171875, 0.242188, 0.3125},
identifierstyle=,
commentstyle=\color[rgb]{0.257813, 0.15625, 0},
stringstyle=\itshape\color[rgb]{0.0195313, 0.195313, 0.0117188},
numbers=left,
numberstyle=\tiny,
stepnumber=1,
breaklines=true,
frame=none,
showstringspaces=false,
tabsize=4,
backgroundcolor=\color[rgb]{0.98,0.98,0.98},
captionpos=b,
float=htbp,
language=Python,
xleftmargin=15pt,
xrightmargin=15pt
}
%(deutsche) Sonderzeichen
\lstset{literate=%
{Ä}{{\"A}}1
{Ö}{{\"O}}1
{Ü}{{\"U}}1
{ä}{{\"a}}1
{ö}{{\"o}}1
{ü}{{\"u}}1
{ß}{{\ss}}1
}
%Verwendung im Text:
%-> \begin{lstlisting}[language=Python,firstnumber=27] ... \end{lstlisting}
%-> \begin{lstlisting}[language=Python,numbers=none] ... \end{lstlisting}
%-> \lstinline[language=JAVA]{...}
\ No newline at end of file
......@@ -24,6 +24,7 @@ while ~isempty(bl)
ad(ta)=1; % flag new entries of al
cand(ta,bc)=1; % indicate candidate (superfluous?)
tb=Tb(bc,4:6); % nhbrs of bc
tb=tb(nc==1); % nhbrs of bc which intersect ac
bl=[bl tb(bd(tb)==0)]; % add those nhbrs not treated to bl
bd(tb)=1; % flag new entries of bl
cand(ac,tb)=1; % ac is candidate for new bl entries
......
......@@ -47,127 +47,21 @@ function [P,n,M]=Intersect(X,Y)
% element Y. The numerical challenges are handled by including
% points on the boundary and removing duplicates at the end.
[P,Q,R,n] = Geometry(X,Y);
if size(P,2)>1 % if two or more interior points
n=[1,1,1]; % the triangle is candidate for all
end % neighbors
P=[P Q R];
P=SortAndRemoveDoubles(P); % sort counter clock wise
[P1,Q,R,n] = Geometry(X,Y);
if size(P1,2)>1 % if two or more interior points
n=[1 1 1]; % the triangle is candidate for all
end % neighbors
P=[P1 Q R];
P=SortAndRemoveDoubles(P); % sort counter clock wise
M=zeros(3,3);
if size(P,2)>0
for j=2:size(P,2)-1 % compute interface matrix
M=M+MortarInt(P(:,[1 j j+1]),X,Y);
end
figure(1)
patch(P(1,:),P(2,:),'m') % draw intersection for illustration
for j=2:size(P,2)-1 % compute interface matrix
M=M+MortarInt(P(:,[1 j j+1]),X,Y);
end
patch(P(1,:),P(2,:),'m') % draw intersection for illustration
% H=line([P(1,:) P(1,1)],[P(2,:),P(2,1)]);
% set(H,'LineWidth',3,'Color','m');
pause(0)
end
% figure(2)
% triX = [ X , X(:,1)];
% triY = [ Y , Y(:,1)];
% plot(triX(1,:),triX(2,:),'-',triY(1,:),triY(2,:),'-')
% pause(0)
end
function [P,Q,R,n] = Geometry(U,V)
% GEOMETRY intersection of two triangles
% [P,Q,R,n] = Geometry(U,V) gives the vertices of U inside V (P),
% intersections between the edges of U and those of V (Q), the vertices
% of V inside U (R), and the neighbours of U that also intersect V (n).
v1=V(:,2)-V(:,1); % vector between V1 and V2
v2=V(:,3)-V(:,1); % vector between V1 and V3
d00=v1'*v1; % norm of v1
d01=v1'*v2;
d11=v2'*v2; % norm of v2
id =d00*d11 - d01*d01; % this value is used to normalize quantities
%---Change of Coordinates (sec.2)---%
X = zeros(2,3);
for i=1:3
r=U(:,i)-V(:,1); % vector between ith vertex of triangle U and V1
d02=v1'*r; % projection of r in v1 direction
d12=v2'*r; % projection of r in v2 direction
X(1,i)=(d11*d02-d01*d12)/id; % amount of r in v1 direction
X(2,i)=(d00*d12-d01*d02)/id; % amount of r in v2 direction
% end result: r = X(1,i)*v1 + X(2,i)*v2
end
%---Intersections (sec.3)---%
if nargout>1
Q = zeros(2,9);
indQ = zeros(1,9);
n = [0,0,0];
indR = n; % 1-(0,0); 2-(1,0); 3-(0,1)
indP = ones(1,3);
for i=1:3
[Qtemp,ntemp,indRtemp,sp] = EdgeIntersect(X,i);
ind = (1:3) + 3*(i-1);
Q(:,ind) = Qtemp;
indQ(ind)= ntemp;
n = max(n,ntemp);
indR = max(indR,indRtemp);
indP = indP.*sp;
end
P = U(:,indP==1);
Q = V(:,1) + [v1,v2]*Q(:,indQ==1);
R = V(:,indR==1);
end
end
function [Q,ind,indR,sp]=EdgeIntersect(X,edge)
% EDGEINTERSECT intersection of edges of X with reference line of Y
% [Q,ind,indR,sp] = EdgeIntersect(X,edge) calculates the intersections
% between all edges of X and the reference line of Y indicated by edge
% and stores the values in Q.
% Boolean indicators for neighbours of U intersecting V, vertices of V
% inside U, and vertices of U inside V are stored in ind, indR and sp,
% respectively.
% Particulars for each edge of reference triangle
if edge==1 % v1
Param1 = @(x,y) y; Param2 = @(x,y) x; % parametrization of edge
Retrieve = @(q) [q;0]; % inversion of parametrization
R1 = 1; R2 = 2; % vertices of V on edge
elseif edge==2 % v2
Param1 = @(x,y) x; Param2 = @(x,y) y;
Retrieve = @(q) [0;q];
R1 = 1; R2 = 3;
elseif edge==3 % line connecting v1 and v2
Param1 = @(x,y) 1-x-y; Param2 = @(x,y) 0.5*(1-x+y);
Retrieve = @(q) [1-q;q];
R1 = 2; R2 = 3;
end
Q = zeros(2,3); ind = zeros(1,3); indR = ind;
p = Param1(X(1,:),X(2,:)); q = NaN*ones(1,3);
sp= sign(p)>=0;
for i = 1:3
j = mod(i,3) + 1;
if sp(i) ~= sp(j) % test if vertices of X lie on opposite sides of edge
if isnan(q(i))% test if q-coord. of vertex already calculated
q(i) = Param2(X(1,i),X(2,i));
end
if isnan(q(j))
q(j) = Param2(X(1,j),X(2,j));
end
q0=(q(i)*p(j) - q(j)*p(i))/(p(j)-p(i)); % calculate intersection
if q0>=0 && q0<=1 % test if intersection lies on edge of Y
Q(:,i) = Retrieve(q0);
ind(i) = 1;
end
if exist('qold','var') % test if 1st intersect. already calculated
indR(R1) = indR(R1) | sign(qold)~=sign(q0);
indR(R2) = indR(R2) | sign(qold-1)~=sign(q0-1);
end
qold = q0;
end
pause(0)
end
end
......
......@@ -24,6 +24,7 @@ while ~isempty(bl)
ad(ta)=1; % flag new entries of al
cand(ta,bc)=1; % indicate candidate (superfluous?)
tb=Tb(bc,4:6); % nhbrs of bc
tb=tb(nc==1); % nhbrs of bc which intersect ac
bl=[bl tb(bd(tb)==0)]; % add those nhbrs not treated to bl
bd(tb)=1; % flag new entries of bl
cand(ac,tb)=1; % ac is candidate for new bl entries
......@@ -32,42 +33,6 @@ while ~isempty(bl)
end
end
% bl =1; % bl: list of triangles of Tb to treat
% bil=1; % bil: list of triangles Ta to start with
% bd=zeros(size(Tb,1)+1,1); % bd: flag for triangles in Tb treated
% bd(end)=1; % guard, to treat boundaries
% bd(1)=1; % mark first triangle in b list.
% M=sparse(size(Nb,2),size(Na,2));
% while ~isempty(bl)
% bc=bl(1); bl=bl(2:end); % bc: current triangle of Tb
% al=bil(1); bil=bil(2:end); % triangle of Ta to start with
% ad=zeros(size(Ta,1)+1,1); % same as for bd
% ad(end)=1;
% ad(al)=1;
% n=[0 0 0]; % triangles intersecting with neighbors
% while ~isempty(al)
% ac=al(1); al=al(2:end); % take next candidate
% [P,nc,Mc]=Intersect(Nb(:,Tb(bc,1:3)),Na(:,Ta(ac,1:3)));
% if ~isempty(P) % intersection found
% M(Tb(bc,1:3),Ta(ac,1:3))=M(Tb(bc,1:3),Ta(ac,1:3))+Mc;
% t=Ta(ac,3+find(ad(Ta(ac,4:6))==0));
% al=[al t]; % add neighbors
% ad(t)=1;
% n(nc>0)=ac; % ac is starting candidate for neighbor
% end
% end
% tmp=find(bd(Tb(bc,4:6))==0); % find non-treated neighbors
% idx=find(n(tmp)>0); % take those which intersect
% t=Tb(bc,3+tmp(idx));
% bl=[bl t]; % and add them
% bil=[bil n(tmp(idx))]; % with starting candidates Ta
% bd(t)=1;
% end
% end
% try same without advancing front - is the problem the intersection or the
% advancing front?
function [P,n,M]=Intersect(X,Y)
% INTERSECT intersection of two triangles and mortar contribution
% [P,n,M]=Intersect(X,Y); computes for the two given triangles X and
......@@ -79,15 +44,12 @@ function [P,n,M]=Intersect(X,Y)
% points on the boundary and removing duplicates at the end.
[P,Q,R,n] = RefFree(Y,X);
if size(P,2)>1 %|| size(R,2)>1 % if two or more interior points
if size(P,2)>1 % if two or more interior points
n=[1,1,1]; % the triangle is candidate for all
end % neighbors
P=[P Q R];
P=SortAndRemoveDoubles(P); % sort counter clock wise
% [Pn,Qn,Rn,nn] = Geometry(X,Y);
% Pn = SortAndRemoveDoubles([Pn,Qn,Rn]);
M=zeros(3,3);
if size(P,2)>0
for j=2:size(P,2)-1 % compute interface matrix
......@@ -99,147 +61,6 @@ if size(P,2)>0
% set(H,'LineWidth',3,'Color','m');
pause(0)
end
% if size(P,2)>0 && size(Pn,2)>0
% polyP = [ P, P(:,1)];
% polyPn= [Pn,Pn(:,1)];
% plot(polyP(1,:),polyP(2,:),'r*',polyPn(1,:),polyPn(2,:),'bo','markersize',10)
% % triX = [ X , X(:,1)];
% % triY = [ Y , Y(:,1)];
% % plot(triX(1,:),triX(2,:),'-',triY(1,:),triY(2,:),'-')
% pause(0.5)
% end
end
function [P,Q,R,n]=RefFree(U,V)
% REFFREE intersection of two triangles
W=V(:,[2,3,1])-V;
n=zeros(1,3); indR=n; indP=ones(1,3);
Q=zeros(2,9); indQ=zeros(1,9);
for i=1:3
j=mod(i,3)+1; w=W(:,i); wt=[-w(2);w(1)];
if wt'*W(:,j)<0
wt=-wt;
end
qp=[w,wt]\(U-V(:,i)); sp=sign(qp(2,:))>=0; qold=NaN;
for k=1:3
l=mod(k,3)+1; m=k+3*(i-1);
if sp(k)~=sp(l) % vertices on opposite sides of edge?
q0=(qp(1,k)*qp(2,l)-qp(1,l)*qp(2,k))/(qp(2,l)-qp(2,k));
if q0>=0 && q0<=1 % intersection lies on edge of Y?
Q(:,m)=V(:,i)+q0*w; indQ(m)=1;
n(i)=1;
end
if isnan(qold) % 1st intersection not calculated?
qold=q0;
else
indR(i)=indR(i)|sign(qold)~=sign(q0); % 0 in [q0,qold]?
indR(j)=indR(j)|sign(qold-1)~=sign(q0-1); % 1 in [q0,qold]?
end
end
end
indP=indP.*sp;
end
P=U(:,indP==1); Q=Q(:,indQ==1); R=V(:,indR==1);
end
function [P,Q,R,n] = Geometry(U,V)
% GEOMETRY intersection of two triangles
% [P,Q,R,n] = Geometry(U,V) gives the vertices of U inside V (P),
% intersections between the edges of U and those of V (Q), the vertices
% of V inside U (R), and the neighbours of U that also intersect V (n).
v1=V(:,2)-V(:,1); % vector between V1 and V2
v2=V(:,3)-V(:,1); % vector between V1 and V3
d00=v1'*v1; % norm of v1
d01=v1'*v2;
d11=v2'*v2; % norm of v2
id =d00*d11 - d01*d01; % this value is used to normalize quantities
%---Change of Coordinates (sec.2)---%
X = zeros(2,3);
for i=1:3
r=U(:,i)-V(:,1); % vector between ith vertex of triangle U and V1
d02=v1'*r; % projection of r in v1 direction
d12=v2'*r; % projection of r in v2 direction
X(1,i)=(d11*d02-d01*d12)/id; % amount of r in v1 direction
X(2,i)=(d00*d12-d01*d02)/id; % amount of r in v2 direction
% end result: r = X(1,i)*v1 + X(2,i)*v2
end
%---Intersections (sec.3)---%
if nargout>1
Q = zeros(2,9);
indQ = zeros(1,9);
n = [0,0,0];
indR = n; % 1-(0,0); 2-(1,0); 3-(0,1)
indP = ones(1,3);
for i=1:3
[Qtemp,ntemp,indRtemp,sp] = EdgeIntersect(X,i);
ind = (1:3) + 3*(i-1);
Q(:,ind) = Qtemp;
indQ(ind)= ntemp;
n = max(n,ntemp);
indR = max(indR,indRtemp);
indP = indP.*sp;
end
P = U(:,indP==1);
Q = V(:,1) + [v1,v2]*Q(:,indQ==1);
R = V(:,indR==1);
end
end
function [Q,ind,indR,sp]=EdgeIntersect(X,edge)
% EDGEINTERSECT intersection of edges of X with reference line of Y
% [Q,ind,indR,sp] = EdgeIntersect(X,edge) calculates the intersections
% between all edges of X and the reference line of Y indicated by edge
% and stores the values in Q.
% Boolean indicators for neighbours of U intersecting V, vertices of V
% inside U, and vertices of U inside V are stored in ind, indR and sp,
% respectively.
% Particulars for each edge of reference triangle
if edge==1 % v1
Param1 = @(x,y) y; Param2 = @(x,y) x; % parametrization of edge
Retrieve = @(q) [q;0]; % inversion of parametrization
R1 = 1; R2 = 2; % vertices of V on edge
elseif edge==2 % v2
Param1 = @(x,y) x; Param2 = @(x,y) y;
Retrieve = @(q) [0;q];
R1 = 1; R2 = 3;
elseif edge==3 % line connecting v1 and v2
Param1 = @(x,y) 1-x-y; Param2 = @(x,y) 0.5*(1-x+y);
Retrieve = @(q) [1-q;q];
R1 = 2; R2 = 3;
end
Q = zeros(2,3); ind = zeros(1,3); indR = ind;
p = Param1(X(1,:),X(2,:)); q = NaN*ones(1,3);
sp= sign(p)>=0;
for i = 1:3
j = mod(i,3) + 1;
if sp(i) ~= sp(j) % test if vertices of X lie on opposite sides of edge
if isnan(q(i))% test if q-coord. of vertex already calculated
q(i) = Param2(X(1,i),X(2,i));
end
if isnan(q(j))
q(j) = Param2(X(1,j),X(2,j));
end
q0=(q(i)*p(j) - q(j)*p(i))/(p(j)-p(i)); % calculate intersection
if q0>=0 && q0<=1 % test if intersection lies on edge of Y
Q(:,i) = Retrieve(q0);
ind(i) = 1;
end
if exist('qold','var') % test if 1st intersect. already calculated
indR(R1) = indR(R1) | sign(qold)~=sign(q0);
indR(R2) = indR(R2) | sign(qold-1)~=sign(q0-1);
end
qold = q0;
end
end
end
function P=SortAndRemoveDoubles(P)
......
@@ -8,13 +8,11 @@
% element Y. The numerical challenges are handled by including
% points on the boundary and removing duplicates at the end.
-[P,n]=EdgeIntersections(X,Y);
-P1=PointsOfXInY(X,Y);
+[P1,Q,R,n] = Geometry(X,Y);
if size(P1,2)>1 % if two or more interior points
n=[1 1 1]; % the triangle is candidate for all
end % neighbors
-P=[P P1];
-P=[P PointsOfXInY(Y,X)];
+P=[P1 Q R];
P=SortAndRemoveDoubles(P); % sort counter clock wise
M=zeros(3,3);
if size(P,2)>0
function [P,Q,R,n]=RefFree(U,V)
% REFFREE intersection of two triangles
% [P,Q,R,n] = RefFree(U,V) produces the vertices of U inside V (P), the
% intersections between the edges of U and V (Q), the vertices of V
% inside U (R) and the neighbours of U that also intersect V (n).
% The result is the same as Geometry(U,V) up to error on the order of
% machine epsilon. The calculations are performed without change of
% coordinates to a reference triangle.
W=V(:,[2,3,1])-V;
W=V(:,[2,3,1])-V; % vectors between vertices of V
n=zeros(1,3); indR=n; indP=ones(1,3);
Q=zeros(2,9); indQ=zeros(1,9);
for i=1:3
j=mod(i,3)+1; w=W(:,i); wt=[-w(2);w(1)];
if wt'*W(:,j)<0
wt=-wt;
j=mod(i,3)+1; w=W(:,i); wt=[-w(2);w(1)]; % find w and vector perp. to w
if wt'*W(:,j)<0 % if V is not on + side of w
wt=-wt; % change sign of wt
end
qp=[w,wt]\(U-V(:,i)); sp=sign(qp(2,:))>=0; qold=NaN;
for k=1:3
......@@ -16,7 +22,7 @@ for i=1:3
q0=(qp(1,k)*qp(2,l)-qp(1,l)*qp(2,k))/(qp(2,l)-qp(2,k));
if q0>=0 && q0<=1 % intersection lies on edge of Y?
Q(:,m)=V(:,i)+q0*w; indQ(m)=1;
n(i)=1; % nhbrs of V intersect U?
n(k)=1; % nhbrs of U intersect V
end
if isnan(qold) % 1st intersection not calculated?
qold=q0;
......
N=1000;
error_geo=zeros(N,1);
error_free=error_geo;
aa=pi*linspace(0,1,N+2); aa=aa(2:end-1);
eps_array = eps * 2.^(-4);
eps_contours = zeros(N,length(eps_array));
for i = 1:N
a=aa(i);
% d=0.5*sqrt(1-0.5*(1-cos(a)));
d=0.5;
V=[0,1,cos(a);0,0,sin(a)]; % triangle to make reference
U=d*[cos(a/2),cos(a/2 + 2*pi/3),cos(a/2 + 4*pi/3);
sin(a/2),sin(a/2 + 2*pi/3),sin(a/2 + 4*pi/3)];
[P,Q,R,~] = Geometry(U,V);
P_geo = SortAndRemoveDoubles([P Q R]);
area_geo = polyarea(P_geo(1,:)',P_geo(2,:)');
[P,Q,R,~] = RefFree(U,V);
P_free = SortAndRemoveDoubles([P Q R]);
area_free = polyarea(P_free(1,:)',P_free(2,:)');
area_exact= d^2 * sin(a/2) * sin(pi/6) / sin(pi - pi/6 - a/2);
if cos(a/2)<d
area_exact= area_exact - (d-cos(a/2))^2*tan(pi/6);
end
error_geo(i) = abs(area_geo-area_exact)/abs(area_exact);
error_free(i)= abs(area_free-area_exact)/abs(area_exact);
eps_contours(i,:) = eps_array/abs(area_exact);
triU = [U, U(:,1)];
triV = [V, V(:,1)];
plot(triU(1,:),triU(2,:),'b-',triV(1,:),triV(2,:),'r-')
patch(P_geo(1,:),P_geo(2,:),'m')
axis([-1,1,-1,1])
axis square
pause(0)
end
figure(1)
semilogy(aa/pi,error_geo,'r*',aa/pi,error_free,'bo',aa/pi,eps_contours,'k-')
xlabel('\alpha/\pi')
ylabel('Rel. error')
% ylabel('Abs. error')
%
% figure(2)
% plot(aa/pi,error_geo>error_free,'k^')
% xlabel('\alpha/\pi')
% ylabel('Geometry<RefFree')
\ No newline at end of file
[N,T,Np]=CircleMesh(20,eps);
[N,T,Np]=CircleMesh(20,5*eps);
%%
figure(1)
......