...
 
Commits (8)
\documentclass[a4paper, 12pt]{article}
\usepackage{preamble}
\usepackage[french]{babel}
\usepackage[utf8]{inputenc}
\usepackage{tikz}
\usepackage[margin=2cm]{geometry}
\title{S\'erie 1, Exercices 4 et 5}
\date{}
\begin{document}
\maketitle
D'abord, qu'est-ce que c'est l'ordre d'un schéma des différences finies (SDF)?
Un SDF est une approximation num\'erique de la dérivée.
L'ordre définit le comportement de l'erreur quand la distance entre les points devient plus petite.
Un SDF est défini par les points utilisés à construire l'approximation, en anglais appelé le 'stencil'.
Pour les deuxièmes dérivées, vous devriez conna\^itre deux en dimension un.
\begin{figure}[h!]
\begin{subfigure}{0.4\textwidth}
\centering
\begin{tikzpicture}
\filldraw[black] (-1,0) circle (3pt) node[above] {$u_{i-1}$};
\filldraw[black] (1,0) circle (3pt) node[above] {$u_{i+1}$};
\draw[black, thick] (-1,0) -- (1,0);
\draw[black,thick] (0,0) circle (3pt) node[above] {$u_i$};
\end{tikzpicture}
\caption{Centr\'ee}
\end{subfigure}
\begin{subfigure}{0.4\textwidth}
\centering
\begin{tikzpicture}
\filldraw[black] (-2,0) circle (3pt) node[above] {$u_{i-2}$};
\filldraw[black] (-1,0) circle (3pt) node[above] {$u_{i-1}$};
\draw[black,thick] (0,0) circle (3pt) node[above] {$u_i$};
\draw[black,thick] (-2,0) -- (0,0);
\end{tikzpicture}
\caption{D\'ecentr\'ee}
\end{subfigure}
\caption{Deux SDF en dimension une.
Les deux approximent la deuxi\`eme d\'eriv\'ee de $u(x)$ au point $x_i$ ($u(x_i) = u_i$).}
\end{figure}
En dimension plus \'el\'ev\'ee, les 'stencils' deviennent plus compliqu\'es.
\begin{figure}[h!]
\begin{subfigure}{0.4\textwidth}
\centering
\begin{tikzpicture}
\filldraw[black] (-1,0) circle (3pt) node[above] {$u_{i-1,j}$};
\filldraw[black] (0,-1) circle (3pt) node[right] {$u_{i,j-1}$};
\filldraw[black] (1,0) circle (3pt) node[above right] {$u_{i+1,j}$};
\filldraw[black] (0,1) circle (3pt) node[right] {$u_{i,j+1}$};
\draw[black,thick] (0,0) circle (3pt) node[above right] {$u_{i,j}$};
\draw[black,thick] (-1,0) -- (1,0);
\draw[black,thick] (0,-1) -- (0,1);
\end{tikzpicture}
\caption{Centr\'ee, en dimension deux}
\end{subfigure}
\begin{subfigure}{0.4\textwidth}
\centering
\begin{tikzpicture}
\filldraw[black] (-1,0,0) circle (3pt);
\filldraw[black] (0,-1,0) circle (3pt);
\filldraw[black] (0,0,-1) circle (3pt);
\filldraw[black] (1,0,0) circle (3pt);
\filldraw[black] (0,1,0) circle (3pt);
\filldraw[black] (0,0,1) circle (3pt);
\draw[black,thick] (0,0,0) circle (3pt);
\draw[black,thick] (-1,0,0) -- (1,0,0);
\draw[black,thick] (0,-1,0) -- (0,1,0);
\draw[black,thick] (0,0,-1) -- (0,0,1);
\end{tikzpicture}
\caption{Centr\'ee, en dimension trois}
\end{subfigure}
\caption{Les SDF centr\'ees en dimensions deux et trois.
Celui en dimension deux est appel\'e l'\'etoile \`a cinque points et est le SDF du Laplacien discret en dimension deux (des diff\'erences finies; il y a d'autres representations du Laplacien discret, qui ne nous concernent pas).}
\end{figure}
Rappeler en dimension une avec un maillage uniforme:
\begin{align*}
u(x_{i+1}) && = u(x_i) & + u'(x_i) (x_{i+1} - x_i) &&+ u''(x_i) \frac{(x_{i+1}-x_i)^2}{2} && + u^{(3)}(x_i) \frac{(x_{i+1}-x_i)^3}{6} && + \order{h^4}\\
&& = u(x_i) & + u'(x_i) h && + u''(x_i) \frac{h^2}{2} && + u^{(3)}(x_i) \frac{h^3}{6} && + \order{h^4}, \\
u(x_{i-1}) && = u(x_i) & + u'(x_i) (x_{i-1} - x_i) && + u''(x_i) \frac{(x_{i-1}-x_i)^2}{2} && + u^{(3)}(x_i) \frac{(x_{i-1}-x_i)^3}{6} && + \order{h^4}\\
&& = u(x_i) & - u'(x_i) h && + u''(x_i) \frac{h^2}{2} && - u^{(3)}(x_i) \frac{h^3}{6} && + \order{h^4}, \\ \hline
u(x_{i+1}) + u(x_{i-1}) && = 2u(x_i) & && + u''(x_i) h^2 && && + \order{h^4},
\end{align*}
parce que dans en maillage uniforme $x_{i+1} - x_i = -(x_{i-1} - x_i) = h$ pour tous les $i$.
Donc, on peut \'ecrire la deuxi\`eme d\'eriv\'ee de $u$ sur $x_i$ comme:
\begin{equation}
u''(x_i) = \frac{u(x_{i+1}) - 2 u(x_i) + u(x_{i-1})}{h^2} + \frac{\order{h^4}}{h^2}
\end{equation}
qui est un SDF d'ordre deux.
En dimension deux, le Laplacien est le somme de deux deuxi\`emes d\'eriv\'ees: $\Delta u = u_{xx} + u_{yy}$.
On peut utiliser le SDF pr\'ec\'edent dans les deux directions, $x$ et $y$:
\begin{equation}
\begin{split}
\Delta u & = \frac{u_{i+1,j} - 2u_{i,j} + u_{i-1,j}}{h^2} + \order{h^2} + \frac{u_{i,j+1} - 2 u_{i,j} + u_{i,j-1}}{h^2} + \order{h^2} \\
& = \frac{u_{i+1,j} + u_{i,j+1} - 4 u_{i,j} + u_{i-1,j} + u_{i,j-1}}{h^2} + \order{h^2}.
\end{split}
\end{equation}
En dimension trois:
\begin{equation}
\Delta u = \frac{u_{i+1,j,k} + u_{i,j+1,k} + u_{i,j,k+1} - 6u_{i,j,k} + u_{i-1,j,k} + u_{i,j-1,k} + u_{i,j,k-1}}{h^2} + \order{h^2}.
\end{equation}
Le Laplacien discret en dimension une ($D^2$) est juste une matrice.
On consid\'ere un vecteur $u$ qui represent une fonction $u(x)$:
\begin{equation*}
D^2 u = \frac{1}{h^2} \begin{bmatrix} -2 & 1 & & \\ 1 & -2 & 1 & \\ & & \ddots & \end{bmatrix} \begin{bmatrix} u_1 \\ u_2 \\ \vdots \\ u_n \end{bmatrix} = \begin{bmatrix} \frac{-2u_1 + u_2}{h^2} \\ \frac{u_1 - 2u_2 + u_3}{h^2} \\ \vdots \\ \frac{u_{n-1} - 2u_n}{h^2} \end{bmatrix}.
\end{equation*}
Pour cr\'eer le Laplacien discret en dimension deux en Matlab, on veut premi\`erement \'ecrire une fonction $u(x,y)$ comme un vecteur.
\begin{equation*}
u = \begin{bmatrix} u_{1,1} & \dots & u_{1,n} \\ \vdots & \ddots & \vdots \\ u_{n,1} & \dots & u_{n,n} \end{bmatrix} \to
\begin{bmatrix} u_{1,1} \\ \vdots \\ u_{n,1} \\ \hline u_{1,2} \\ \vdots \\ u_{n,2} \\ \hline \vdots \\ \hline u_{1,n} \\ \vdots \\ u_{n,n} \end{bmatrix} = u(:) .
\end{equation*}
Consid\'erer le Laplacien discret comme somme de deux matrices qui representent les operateurs de la deux\`ieme d\'eriv\'ee dans les deux directions ($x$ et $y$): $L = D_{xx} + D_{yy}$.
La matrice $D_{xx}$ agit juste sur les colonnes de $u$, $u(:,i)$ pour $i=1,...,n$.
Ca veut dire que:
\begin{equation*}
D_{xx} u(:) = \begin{bmatrix} \frac{-2u_{1,1} + u_{2,1}}{h^2} \\ \frac{u_{1,1} - 2u_{2,1} + u_{3,1}}{h^2} \\ \vdots \\ \frac{u_{n-1,1} - 2u_{n,1}}{h^2} \\ \hline \\
\frac{-2u_{1,2} + u_{2,2}}{h^2} \\ \vdots \\ \frac{u_{n-1,2} - 2u_{n,2}}{h^2} \\ \hline \vdots \\ \hline
\frac{-2u_{1,n} + u_{2,n}}{h^2} \\ \vdots \\ \frac{u_{n-1,n} - 2u_{n,n}}{h^2} \end{bmatrix}
= \begin{bmatrix} D^2 & & \\ & \ddots & \\ & & D^2 \end{bmatrix} u(:) = (I \otimes D^2) u(:).
\end{equation*}
o\`u $I$ est la matrice identit\'e et $\otimes$ est la multiplication de Kronecker.
Pour la matrice $D_{yy}$, noter que la deuxi\`eme matrice dans le produit de Kronecker agit sur les colonnes de $u$ et la premi\`ere sur les lignes.
Alors, $D_{yy} = D^2 \otimes I$.
On peut confirmer par un calcul direct que $D_{yy} u(:)$ donne les deuxi\`emes d\'eriv\'ees correctes.
En fin, $L = (I \otimes D^2) + (D^2 \otimes I)$.
En dimension trois, on peut \'etendre la formule:
$L = (I \otimes I \otimes D^2) + (I \otimes D^2 \otimes I) + (D^2 \otimes I \otimes I)$.
Chaque endroit dans les produits de Kronecker signifie une direction: $x$, $y$ et $z$.
%On peut aussi consid\'erer $D^2 u$:
%\begin{equation*}
%D^2 u = \begin{bmatrix} \frac{-2u_{1,1} + u_{2,1}}{h^2} & \frac{-2u_{1,2} + u_{2,2}}{h^2} & \dots & \frac{-2u_{1,n} + u_{2,n}}{h^2} \\
% \frac{u_{1,1} - 2u_{2,1} + u_{3,1}}{h^2} & \vdots & & \vdots \\
% \vdots & \vdots & & \vdots \\
% \frac{u_{n-1,1} - 2u_{n,1}}{h^2} & \frac{u_{n-1,2} - 2u_{n,2}}{h^2} & \dots & \frac{u_{n-1,n} - 2u_{n,n}}{h^2} \end{bmatrix}.
%\end{equation*}
%Chaque \'el\'ement est la deuxi\`eme d\'eriv\'ee (appoximation) en direction $x$.
%
%Pour la matrice $D_{yy}$, noter qu'on peut consid\'erer $D^2 u^\top$, et chaque \'el\'ement de ce produit est la deuxi\`eme d\'eriv\'ee (appoximation) en direction $y$.
\end{document}
\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}{1} % 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}
\begin{enumerate}
\item $$\norm{A}_1 = \sup \frac{\norm{Ax}_1}{\norm{x}_2} = \sup \frac{\sum_i \abs{(Ax)_i}}{\sum_j \abs{x_j}}$$
On peut représenter $A$ par ses éléments $a_{ij}$ ($i$ pour les ranges, $j$ pour les colonnes).
Alors, $(Ax)_i$, le $i$-ième range du produit $Ax$, est $\sum_j a_{ij} x_j$.
\begin{align*}
\sum_i \abs{\sum_j a_{ij} x_j} & \leq \sum_i \sum_j \abs{a_{ij}} \abs{x_j} \\
& = \sum_j \abs{x_j} \sum_i \abs{a_{ij}} \\
& \leq \max_k \sum_i \abs{a_{ik}} \sum_j \abs{x_j} \\
\implies \norm{A}_1 & \leq \max_k \sum_i \abs{a_{ik}}
\end{align*}
Il suffit de prouver qu'il existe un $x \in \bbc^n$ tel que $\norm{A}_1 = \max_k \sum_i \abs{a_{ik}}$.
Soit $x = e_k = (0 ... 1 ... 0)^\top$.
Donc, $Ax$ est le $k$-ième colonne de $A$ et $\norm{Ax}_1/\norm{x}_1 = \sum_i \abs{a_{ik}}$.
C'est vrai pour tous les $k$, alors $\max_k \norm{Ae_k}_1 \leq \sup \norm{Ax}_1/\norm{x}_1 = \norm{A}_1$.
\item $$\norm{A}_\infty = \sup \frac{\norm{Ax}_\infty}{\norm{x}_\infty} = \sup \frac{\max_i \abs{\sum_j a_{ij} x_j}}{\max_j \abs{x_j}}$$
\begin{align*}
\abs{\sum_j a_{ij} x_j} & \leq \sum_j \abs{a_{ij}} \abs{x_j} \\
& \leq \max_k \abs{x_k} \sum_j \abs{a_{ij}} \\
& = \norm{x}_\infty \sum_j \abs{a_{ij}} \\
\implies \max_i \abs{\sum_j a_{ij} x_j} & \leq \norm{x}_\infty \max_i \sum_j \abs{a_{ij}} \\
\implies \norm{A}_\infty & \leq \max_i \sum_j \abs{a_{ij}}
\end{align*}
Soit $x = (1 ... 1)^\top$, alors $Ax = ( \sum_j a_{1j} ... \sum_j a_{nj} )^\top$ et $\norm{Ax}_\infty/\norm{x}_\infty = \max_i \abs{\sum_j a_{ij}} \leq \max_i \sum_j \abs{a_{ij}}$, qui n'est pas l'inégalité qu'on veux.
Choisissons $1 \leq k \leq m$ et soit $x = ( \sign(a_{k1}) ... \sign(a_{km}) )^\top$.
Donc, $Ax = ( ? ... \sum_j \abs{a_{kj}} ... ? )^\top$ et $\norm{Ax}_\infty \geq \sum_j \abs{a_{kj}}$ pour $k$ quelconque.
Alors, $\sup \norm{Ax}_\infty \geq \sum_j \abs{a_{ij}}$ pour tout les $i$, et certainement pour le maximum sur les mêmes.
\item $$\norm{Ax}_2 = \sqrt{(Ax)^*(Ax)} = \sqrt{x^* A^* A x}$$
La matrice $A^* A$ est carr\'ee et, en fait, Hermitian et (positive semi-definite).
Alors, il existe une décomposition de celle-l\`a en éléments propres telle que $x = \sum_i \alpha_i v_i$, o\`u $\set{v_i}$ est l'ensemble des vecteurs propres qui sont (orthonormal).
\begin{align*}
\norm{Ax}_2 & = \sqrt{\sum_i \sum_j \alpha_i v_i^* A^* A \alpha_j v_j} \\
& = \sqrt{\sum_i \sum_j \lambda_i \alpha_i \alpha_j v_i^* v_j} \\
& = \sqrt{\sum_i \lambda_i \alpha_i^2} \\
& \leq \max_k \sqrt{\lambda_k} \sqrt{\sum_i \alpha_i^2} \\
& = \sqrt{\rho(A^* A)} \norm{x}_2
\end{align*}
Certainement, $\norm{A}_2 \leq \sqrt{\rho(A^* A)}$.
De plus, si $x = v_i$ tel que $\lambda_i = \rho(A^* A)$, alors $\norm{Ax}_2 = \sqrt{v_i^* A^* A v_i} = \sqrt{\rho(A^* A) v_i^* v_i} = \sqrt{\rho(A^* A)} \norm{x}_2$ et $\norm{A}_2 \geq \sqrt{\rho(A^* A)}$.
\item $$trace(A^*A) = \sum_i \begin{pmatrix} a_{i1}^* & \dots & a_{im}^* \end{pmatrix} \begin{pmatrix} a_{i1} \\ \vdots \\ a_{im} \end{pmatrix} = \sum_i \sum_j \abs{a_{ij}}^2 = \norm{A}_F$$
\end{enumerate}
\end{ex}
\begin{ex}
Rappeler l'inégalité de Cauchy-Schwarz:
\begin{equation}
\abs{\sum_i u_i^* v_i}^2 \leq \sum_i \abs{u_i}^2 \sum_j \abs{v_j}^2 .
\end{equation}
\begin{align*}
\norm{A B}_F & = \sum_i \sum_j \abs{(AB)_{ij}}^2 \\
& = \sum_i \sum_j \abs{\sum_k a_{ik} b_{kj}}^2 \\
& \leq \sum_i \sum_j \sum_k \abs{a_{ik}}^2 \sum_l \abs{b_{lj}}^2 \\
& = \sum_i \sum_k \abs{a_{ik}}^2 \sum_j \sum_l \abs{b_{lj}}^2 \\
& = \norm{A}_F \norm{B}_F
\end{align*}
\end{ex}
\begin{ex}
Premièrement, si $Av = \rho(A) v$, alors $\rho(A) \norm{[v ... v]} = \norm{\rho(A) [v ... v]} = \norm{A [v ... v]} \leq \norm{A} \norm{[v ... v]} \implies \rho(A) \leq \norm{A}$.
Rappeler le forme canonique de Jordan:
\begin{equation*}
A = P \begin{bmatrix} J(\lambda_1) & & \\ & \ddots & \\ & & J(\lambda_k) \end{bmatrix} P^{-1}
\end{equation*}
o\`u $J(\lambda_i)$ est une matrice triangulaire supérieure de taille $n_i \times n_i$ telle que
\begin{equation*}
J(\lambda_i) = \begin{bmatrix} \lambda_i & 1 & & \\ & \ddots & \ddots & \\ & & \lambda_i & 1 \\ & & & \lambda_i \end{bmatrix}.
\end{equation*}
Noter que $\sum_i n_i = n$.
Si on prend la 1-norme de $P^{-1} A P$ on a $\norm{P^{-1} A P}_1 \leq \rho(A)+1$, parce que chaque colonne contient seulement un valeur propre de $A$ et possiblement 1.
On veut des matrices \`a changer les éléments sur la diagonale de $J(\lambda_i)$ \`a $\epsilon$, mais pas changer la diagonale.
La solution est \`a multiplier $J(\lambda_i)$ sur la gauche et la droit:
\begin{equation*}
\begin{bmatrix} \lambda_i & \epsilon & & \\ & \ddots & \ddots & \\ & & \lambda_i & \epsilon \\ & & & \lambda_i \end{bmatrix} = \begin{bmatrix} 1/\epsilon & & \\ & \ddots & \\ & & 1/\epsilon^{n_i} \end{bmatrix} J(\lambda_i) \begin{bmatrix} \epsilon & & \\ & \ddots & \\ & & \epsilon^{n_i} \end{bmatrix} = D_i(1/\epsilon) J(\lambda_i) D_i(\epsilon).
\end{equation*}
Soit $D(\epsilon)$ la matrice par les blocs $D_i(\epsilon)$ dans le même ordre des $J(\lambda_i)$ dans $P^{-1} A P$.
Définir la norme suivante:
\begin{equation*}
\norm{M} = \norm{D(1/\epsilon) P^{-1} M P D(\epsilon)}_1.
\end{equation*}
Clairement $\norm{A} \leq \rho(A) + \epsilon$.
Il faut prouver qu'elle est une norme.
\begin{lemma}
Si $B$ est une matrice inversible, alors la norme définir par $\norm{A}_* = \norm{B^{-1}AB}$ est vraiment une norme de matrices.
\end{lemma}
\begin{proof}
La matrice $B$ est inversible, et donc pas zéro.
Alors, $\norm{A}_* \geq 0$.
\begin{align*}
\norm{B^{-1}(A+C)B} & = \norm{B^{-1}AB + B^{-1}CB} \\
& \leq \norm{B^{-1}AB} + \norm{B^{-1}CB} \\
\norm{B^{-1}ACB} & = \norm{B^{-1}ABB^{-1}CB} \\
& \leq \norm{B^{-1}AB} \norm{B^{-1}CB} \\
\norm{B^{-1}\lambda A B} & = \abs{\lambda} \norm{B^{-1} A B}
\end{align*}
C'est tous les conditions, donc $\norm{A}_*$ est vraiment une norme.
\end{proof}
Pour le deuxième partie de l'exercice, rappeler que toutes les normes des matrices sont équivalentes.
Ça veut dire qu'il existe $C_p$ tel que $\norm{A}_p \leq C_p \norm{A} \leq C_p (\rho(A) + \epsilon/2)$.
(Soit $\epsilon>0$ quelconque, alors ils existent une norme et un $C_p>0$ qui satisfont cette inégalité.)
Alors,
\begin{align*}
\norm{A^k}_p & \leq C_p \norm{A^k} \\
& \leq C_p \norm{A}^k \\
& \leq C_p (\rho(A) + \epsilon/2)^k \\
\implies \norm{A^k}_p^{1/k} & \leq C_p^{1/k} (\rho(A) + \epsilon/2)
\end{align*}
qui approche $\rho(A) + \epsilon$ quand $k$ tend vers infinité.
Aussi, $\rho(A) = \rho(A^k)^{1/k} \leq \norm{A^k}_p^{1/k}$.
Donc, pour tout les $\epsilon >0$ il existe un $N \in \bbn$ tel que pour tout les $k\geq N$ $\abs{\norm{A^k}_p^{1/k} - \rho(A)} < \epsilon$.
\end{ex}
\begin{ex}
Premièrement, notons que le définition du Laplacian:
$$\Delta u = u_{xx} + u_{yy}.$$
Écrire les polynômes de Taylor pour $u$ sur les quatre points \`a cote de $u(x_i,y_j)=u_{i,j}$.
\begin{align*}
u_{i+1,j} & = u_{i,j} & + u_x(x_i,y_j) h & + u_{xx}(x_i,y_j) h^2/2 & + u_{xxx}(x_i,y_j) h^3/6 & + \order{h^4}, \\
u_{i-1,j} & = u_{i,j} & - u_x(x_i,y_j) h & + u_{xx}(x_i,y_j) h^2/2 & - u_{xxx}(x_i,y_j) h^3/6 & + \order{h^4}, \\
u_{i,j+1} & = u_{i,j} & + u_y(x_i,y_j) h & + u_{yy}(x_i,y_j) h^2/2 & + u_{yyy}(x_i,y_j) h^3/6 & + \order{h^4}, \\
u_{i,j-1} & = u_{i,j} & - u_y(x_i,y_j) h & + u_{yy}(x_i,y_j) h^2/2 & - u_{yyy}(x_i,y_j) h^3/6 & + \order{h^4}, \\ \hline
u_{i+1,j} + u_{i-1,j} + u_{i,j+1} + u_{i,j-1} & = 4 u_{i,j} & + u_{xx}(x_i,y_j)h^2 & + u_{yy}(x_i,y_j)h^2 & & + \order{h^4}.
\end{align*}
Réarranger \`a isoler $\Delta u$:
\begin{equation*}
\Delta u_{i,j} = \frac{u_{i+1,j} + u_{i-1,j} + u_{i,j+1} + u_{i,j-1} - 4 u_{i,j}}{h^2} + \order{h^2}.
\end{equation*}
\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
......@@ -97,7 +97,7 @@ et que ce schéma est d'ordre deux. % quelle que soit la dimension du domaine $d
la dimension du domaine et qui retourne la matrice du Laplacien discret avec
l'en-tête suivante,
\VerbatimInput[firstline=0,lastline=5]{FDLaplacian.m}
\VerbatimInput[firstline=0,lastline=5]{FDLaplacianTrue.m}
\end{enumerate}
......
StudentSolns/
\ No newline at end of file
......@@ -2,7 +2,7 @@
n = 10;
x = linspace(0,1,n+2);
u = sin(pi*x')*sin(pi*x) + 1;
u = sin(pi*x')*sin(pi*x)+1;
% II= eye(n+2); I0 = [zeros(1,n+2); II(2:end-1,:) ; zeros(1,n+2)]; II = II-I0;
% II= kron(II,I0) + kron(I0,II); II = sum(II,2);
......@@ -12,7 +12,15 @@ u = sin(pi*x')*sin(pi*x) + 1;
x = x(2:end-1);
g = ones(4*n,1);
f = -2*pi^2 * sin(pi*x')*sin(pi*x); f = f(:);
% g = @(x,y) 1;
% f = @(x,y) -2*pi^2 * sin(pi*x)*sin(pi*y);
uu = Heat2D(f,g,1);
uu = Heat2D(f,g,0,100);
uu = reshape(uu,n,n);
surf(x,x',uu - u(2:end-1,2:end-1))
pause
uu = Heat2D(f,g,1,100);
uu = reshape(uu,n,n);
surf(x,x',uu - u(2:end-1,2:end-1))
\ No newline at end of file
function [x,X] = BisectSolve(f,x1,x2,maxit,tol)
% BISECTSOLVE solves f(x)=0 using bisection method on interval [x1,x2].
% Function terminates when interval is smaller than tol or number of
% iterations is larger than maxit. Output X is a row vector of the
% iterates.
iter = 1;
a = x1; b = x2;
sa=sign(f(a)); sb=sign(f(b));
X = [];
while iter < maxit && abs(b-a) > tol
x = (a+b)/2;
X = [X x];
sf= sign(f(x));
if sf==0
a=x; b=x;
elseif sf==sa
a=x;
elseif sf==sb
b=x;
end
iter = iter+1;
end
\ No newline at end of file
function [x,X] = FixedPointSolve(A,f,x0,alpha,maxit,tol)
% FIXEDPOINTSOLVE solves a linear system Ax=f using a fixed point iteraiton
% The function terminates when number of iterations exceeds maxit or
% when the relative step size becomes smaller than tol.
% Output X is a matrix with the kth column the kth iterate.
X = x0;
x = x0 - alpha*(A*x0 - f);
iter = 2;
while norm(x-x0)/norm(x) > tol && iter < maxit
x0= x;
x = x0 - alpha*(A*x0 - f);
X = [X x];
iter = iter+1;
end
\ No newline at end of file
function [x,X,C] = GSSolve(A,f,x0,maxit,tol)
% GSSOLVE solves a linear system Ax=f using the Gauss-Seidel method
% The function terminates when number of iterations exceeds maxit or
% when the relative step size becomes smaller than tol.
% Output X is a matrix with the kth column the kth iterate.
M = tril(A);
N = M - A;
C = max(abs(eig(inv(M)*N)));
X = x0;
x = M \ (N*x0 + f);
iter = 2;
while norm(x-x0)/norm(x) > tol && iter < maxit
x0= x;
x = M \ (N*x0 + f);
X = [X x];
iter = iter+1;
end
\ No newline at end of file
function [x,X] = GaussNewton(f,df,x0,tol,maxiter)
%GAUSSNEWTON minimizes norm(f) using the Gauss-Newton method.
% Terminates if step length less than tol or number of iterations exceeds
% maxiter. Output X contains the kth iterate in the kth column.
x = x0 - ( df(x0) \ f(x0) );
X = x;
iter = 2;
while norm(x - x0)>tol && iter < maxiter
x0= x;
x = x0 - ( df(x0) \ f(x0) );
X = [X x];
iter = iter+1;
end
\ No newline at end of file
function [x,X,C] = JacobiSolve(A,f,x0,maxit,tol)
% JACOBISOLVE solves a linear system Ax=f using Jacobi's method
% The function terminates when number of iterations exceeds maxit or
% when the relative step size becomes smaller than tol.
% Output X is a matrix with the kth column the kth iterate.
M = diag(diag(A));
N = M - A;
C = max(abs(eig(inv(M)*N)));
X = x0;
x = M \ (N*x0 + f);
iter = 2;
while norm(x-x0)/norm(x) > tol && iter < maxit
x0= x;
x = M \ (N*x0 + f);
X = [X x];
iter = iter+1;
end
\ No newline at end of file
function [x,X] = Newton(f,df,x0,tol,maxiter)
%NEWTON sovles the equation f(x)=0 using Newton's method with initial guess
%x0.
% Terminates if step length less than tol or number of iterations exceeds
% maxiter. Output X contains the kth iterate in the kth column.
x = x0 - ( df(x0) \ f(x0) );
X = x;
iter = 2;
while norm(x - x0)>tol && iter < maxiter
x0= x;
x = x0 - ( df(x0) \ f(x0) );
X = [X x];
iter = iter+1;
end
\ No newline at end of file
%% Run file for TP8
%% Exercise 1
f1 = @(x) x.^2 - 3*x - 4; % sign error
f2 = @(x) (x+1).*(x+5).^3 .* (x-5).^5;
[x1,X1] = BisectSolve(f1,-3.8,0.5,30,1e-10);
[x2,X2] = BisectSolve(f2,-3.8,0.5,30,1e-10);
figure(1)
plot(1:length(X1),X1,'bo',1:length(X2),X2,'r*')
xlabel('Iteration')
ylabel('x')
figure(2)
semilogy(1:length(X1),abs(f1(X1)),'bo',1:length(X2),abs(f2(X2)),'r*')
xlabel('Iteration')
ylabel('|f(x)|')
figure(3)
semilogy(1:length(X1),abs(X1-x1),'bo',1:length(X2),abs(X2-x2),'r*')
xlabel('Iteration')
ylabel('Error')
P = fit((1:length(X1))',log(abs(X1'-x1)),'poly1');
order1 = exp(P.p1);
P = fit((1:length(X2))',log(abs(X2'-x2)),'poly1');
order2 = exp(P.p1);
%% Exercise 2
n = 100;
lambda_min = 3;
lambda_max = 17;
[U,~] = qr(rand(n));
d = lambda_min + (lambda_max-lambda_min) * rand(n,1);
A = U * diag(d) * U';
x = 2*rand(n,1)-1;
f = A*x;
I = eye(n);
figure(1)
clf
for alpha = -0.1:0.05:0.1
[~,X] = FixedPointSolve(A,f,zeros(n,1),alpha,50,1e-15);
C = norm(I - alpha*A);
[rowX,colX] = size(X);
errX = zeros(colX,1); ind = 1:colX;
for k=ind
errX(k) = norm(X(:,k) - x)/norm(x);
end
figure(1)
semilogy(ind,errX,'.--',ind,C.^ind,'-')
hold on
end
xlabel('Iteration')
ylabel('Relative Error')
legend('-0.1, ex','-0.1, th','\alpha=-0.05','-0.05, th','\alpha=0',...
'0, th','\alpha=0.05','0.05, th','\alpha=0.1','0.1,th')
hold off
figure(2)
clf
[~,XJ,CJ] = JacobiSolve(A,f,zeros(n,1),50,1e-15);
[~,XG,CG] = GSSolve(A,f,zeros(n,1),50,1e-15);
[~,colJ] = size(XJ); indJ = 1:colJ; errJ = zeros(colJ,1);
[~,colG] = size(XG); indG = 1:colG; errG = zeros(colG,1);
for k = indJ
errJ(k) = norm(XJ(:,k) - x)/norm(x);
end
for k = indG
errG(k) = norm(XG(:,k) - x)/norm(x);
end
semilogy(indJ,errJ,'r*--',indJ,CJ.^indJ,'r',indG,errG,'bo--',indG,CG.^indG,'b')
xlabel('Iteration')
ylabel('Relative error')
%% Exercise 3
clf
clear
clc
f = @(x) x.^3 - x - 3;
df= @(x) 3*x.^2 - 1;
[x0,X0] = Newton(f,df,0,1e-15,50);
[x1,X1] = Newton(f,df,1,1e-15,50);
figure(1)
semilogy(1:length(X0),abs(X0 - x1),'r*',1:length(X1),abs(X1-x1),'bo')
xlabel('Iteration')
ylabel('Absolute error')
f = @(x) [ x(1)^3-3*x(1)*x(2)^2-1;
3*x(1)^2*x(2)-x(2)^3 ];
df= @(x) [ 3*(x(1)^2 - x(2)^2), -6*x(1)*x(2);
6*x(1)*x(2), 3*(x(1)^2 - x(2)^2) ];
[x,X] = Newton(f,df,[-1;1],1e-10,100);
[~,col] = size(X); err = zeros(col,1);
for k = 1:col
err(k) = norm(X(:,k) - x);
end
figure(2)
semilogy(1:col,err)
xlabel('Iteration')
ylabel('Absolute error')
dg= @(x) [2*(x(1) - 1);4*(x(2) + 1)^3];
Hg= @(x) [2, 0; 0, 12*(x(2)+1)^2];
[x,X] = Newton(dg,Hg,[0;0],1e-10,100);
[~,col] = size(X); err = zeros(col,1);
for k=1:col
err(k) = norm(X(:,k) - x);
end
figure(3)
semilogy(1:col,err)
xlabel('Iteration')
ylabel('Absolute error')
\ No newline at end of file