Commit eb28a7ee authored by conmccoid's avatar conmccoid

Update for June:

Coding updates to Edge Intersectins
Three months of assistant work added
Initial commit on Stochastic Networks research
parent 6907e333
function A = FDLaplacian(n,d)
% FDLAPLACIAN compute a finite difference Laplacian
%A = FDLaplacian(n,d) computes a finite difference Laplacian on the
%unit interval/square/cube (d=1,2,3) using n interior points
m=n;
if d==1
A=diag(2*ones(1,m)) + diag(ones(m-1,1),1) + diag(ones(m-1,1),-1) ;
end
if d==2
A=diag(4*ones(1,m^2)) + diag(ones(m^2-1,1),1) + diag(ones(m^2-1,1),-1) + diag(ones(m^2-m,1),m)+ diag(ones(m^2-m,1),-m);
end
if d==3
A=diag(6*ones(1,m^3)) + diag(ones(m^3-1,1),1) + diag(ones(m^3-1,1),-1) + diag(ones(m^3-m,1),m)+ diag(ones(m^3-m,1),-m)+diag(ones(m^3-m^2,1),m^2)+ diag(ones(m^3-m^2,1),-m^2);
end
end
A=FDLaplacian(5,1)
A=FDLaplacian(5,2)
A=FDLaplacian(5,3)
\ No newline at end of file
......@@ -6,28 +6,27 @@ function u=Heat2D(f,g,p)
% 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;
n = length(g); n = n/4; ind=1:n; % taille du systeme
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;
D = (n+1)^2*spdiags(ones(n,1)*[1 -2 1],[-1 0 1],n,n);
D = kron(I,D) + kron(D,I);
G = kron(g(ind), I(:,1)) + ... % bord a y=0
kron(g(ind+n),I(:,end)) + ... % bord a y=1
kron(I(:,1), g(ind+2*n)) + ... % bord a x=0
kron(I(:,end),g(ind+3*n)); % bord a x=1
F = f - (n+1)^2 * G; % cote droit
D = (n+1)^2*spdiags(ones(n,1)*[1 -2 1],[-1 0 1],n,n); % Lap. en dim. 1
D = kron(I,D) + kron(D,I); % Lap. en dim. 2
if p==0
if p==0 % Jacobi
N = triu(D,1) + tril(D,-1);
M = 1 / (4 * (n+1)^2);
elseif p==1
elseif p==1 % Gauss-Seidel
N = triu(D,1);
M = D - N; M = -M \ kron(I,I);
end
u0 = kron(ones(n,1),ones(n,1));
err= 1;
tol= eps;
itermax = 1000;
iter = 1;
u0 = kron(ones(n,1),ones(n,1)); % estime initiale
err= 1; iter = 1; % initialisation
tol= eps; itermax = 1e3; % parametres d'iteration
while err > tol && iter < itermax
u = M*N*u0 - M*F;
u0 = u;
......
......@@ -8,6 +8,7 @@ x0= sqrt(5)*[cos(x0);sin(x0)]; % anneau
x1= G*x0; x2= G*x1; % comportement de l'anneau
figure(1)
clf
hold on
plot(x0(1,:),x0(2,:),'b--',... % anneau
x1(1,:),x1(2,:),'r--',... % anneau apres 1 appl. de G
......@@ -51,7 +52,7 @@ ylabel('Vecteurs propres')
axis([0,1,-1,1])
%% Vecteurs propres en dimension deux (ex5)
n = 20; % discretization
n = 10; % discretization
x = linspace(0,1,n+2)'; % carr'e unitaire
k = 2; % nombre des vecteurs propres
......@@ -85,4 +86,34 @@ for i = 1:k^2
surf(x(2:end-1),x(2:end-1)',reshape(V(:,i),n,n))
end
%% Conditions au bord (ex6)
\ No newline at end of file
%% Heat2D (ex6)
n = 10; x = linspace(0,1,n+2); % discretisation
uex= -0.25*(x'-0.5).^2 * (x-0.5).^2; % solution exacte
g = [ uex(2:end-1,1) ;... % bord a y=0
uex(2:end-1,end);... % bord a y=1
uex(1,2:end-1)' ;... % bord a x=0
uex(end,2:end-1)']; % bord a x=1
x = x(2:end-1);
f = -0.5*( (x'-0.5).^2+(x-0.5).^2 ); % 'forcing' sur l'interieur
u = Heat2D(f(:),g,0); % solution approx.
u = reshape(u,n,n); % solution en dim. 2
uex= uex(2:end-1,2:end-1); % l'interieur
figure(1)
subplot(1,3,1) % solution exacte sur l'interieur
surf(x,x',uex)
xlabel('x')
ylabel('y')
zlabel('Exact solution')
subplot(1,3,2) % solution approx. sur l'interieur
surf(x,x',u)
xlabel('x')
ylabel('y')
zlabel('Approximate solution')
subplot(1,3,3) % erreur
surf(x,x',abs(u-uex))
xlabel('x')
ylabel('y')
zlabel('Error')
\ No newline at end of file
......@@ -11,36 +11,45 @@ if n~=m
disp('uu is not square')
end
n = n-2;
ind0 = [1 ; zeros(n,1) ; 1]; ind1 = [0 ; ones(n,1) ; 0];
ind = kron(ind0,ind1) + kron(ind1,ind0);
g = uu(:); g = g(ind==1);
I = speye(n);
ind=1:n;
% ind0 = [1 ; zeros(n,1) ; 1]; ind1 = [0 ; ones(n,1) ; 0];
% ind = kron(ind0,ind1) + kron(ind1,ind0);
% g = uu(:); g = g(ind==1);
%
% ind=1: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;
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;
D = (n+1)^2*spdiags(ones(n,1)*[1 -2 1],[-1 0 1],n,n);
D = kron(I,D) + kron(D,I);
F = uu(2:end-1,2:end-1); F = D*F(:);
if p==0
N = triu(D,1) + tril(D,-1);
M = 1 / (4 * (n+1)^2);
M = -(4 * (n+1)^2);
elseif p==1
N = triu(D,1);
M = D - N; M = -M \ kron(I,I);
M = D - N;
end
u0 = kron(ones(n,1),ones(n,1));
itermax = 100;
err= zeros(itermax,1);
for k = 1:itermax
u = M*N*u0 - M*F;
u = -M \ (N*u0 - F);
u0 = u;
err(k)= norm(reshape(u,n,n) - uu(2:end-1,2:end-1), 'fro');
end
P = fit((1:itermax)',log(err),'poly1');
pest = exp(P.p1);
\ No newline at end of file
pest = exp(P.p1);
% figure(2)
% plot(1:itermax,log(err),'b.',1:itermax,P(1:itermax),'r--')
% xlabel('Iteration')
% ylabel('log(error)')
% legend('log(e_k)',[num2str(P.p2) ' + (' num2str(P.p1) ') k'])
% pause
\ No newline at end of file
%% Laplacian convergence
n = 20;
A = FDLaplacianTrue(n,2); A = -A;
x = linspace(0,1,n);
u = cos(x') * sin(x); u = u(:);
f = A*u;
[v,vk] = SteepestDescent(A,f,ones(n^2,1),1e-10,100);
[vrow,vcol] = size(vk);
error = zeros(vcol,1);
for i=1:vcol
w = vk(:,i) - u;
error(i) = sqrt(w'*A*w);
end
ind = round(0.8*vcol):vcol;
P0= fit((1:vcol)',log(error),'poly1');
P1= fit(ind',log(error(ind)),'poly1');
C0 = exp(P0.p1);
C1 = exp(P1.p1);
Cex= cos(pi/(n+1));
ind = (1-vcol):0;
plot(1:vcol,log(error),'b.',...
1:vcol,log(error(end) * C0 .^ind),'m--',...
1:vcol,log(error(end) * C1 .^ind),'r--',...
1:vcol,log(error(end) * Cex.^ind),'k-')
xlabel('Iteration')
ylabel('log(Error)')
legend('|e_k|_A','Pente (entiere)','Pente (partielle)','Pente (exacte)')
%% Laplacian convergence - fonction de h
nn= 10:5:100;
C = zeros(1,length(nn));
for k = 1:length(nn)
n = nn(k);
A = FDLaplacianTrue(n,2); A = -A;
x = linspace(0,1,n);
u = cos(x') * sin(x); u = u(:);
f = A*u;
[v,vk] = SteepestDescent(A,f,ones(n^2,1),1e-10,100);
[vrow,vcol] = size(vk);
error = zeros(vcol,1);
for i=1:vcol
w = vk(:,i) - u;
error(i) = sqrt(w'*A*w);
end
ind = round(0.8*vcol):vcol;
P = fit(ind',log(error(ind)),'poly1');
C(k)= exp(P.p1);
end
h = 1./(nn+1);
plot(h,C,'b*',h,cos(pi*h),'k--')
xlabel('h')
ylabel('$\frac{\kappa(A)-1}{\kappa(A)+1}$','interpreter','latex')
legend('Approx.','Exacte')
\ No newline at end of file
......@@ -2,44 +2,87 @@
%% Heat2D, estimate spectral radius
nn = 10:5:50;
% Jacobi
nn = 10:5:100;
pest = zeros(length(nn),1);
p = 0;
for k=1:length(nn)
n = nn(k);
x = linspace(0,1,n+2);
u = sin(pi*x')*sin(pi*x) + 1;
[uu,pest(k)] = Heat2D(u,0);
[uu,pest(k)] = Heat2D(u,p);
end
% plot(log(1./(nn+1)),log(pest-1),'r*--')
P = fit(1./(nn'+1),pest,'poly2');
h = 1./(linspace(nn(1),nn(end),1000)+1);
figure(1)
subplot(1,2,1)
plot(1./(nn+1),pest,'b*',h,P(h),'r--',h,cos(pi*h).^(p+1),'k-',...
'MarkerSize',10,'Linewidth',2)
xlabel('h')
ylabel('\rho(M^{-1}N)')
legend('Taux de convergence asymptotique',...
[num2str(P.p3) ' + (' num2str(P.p2) ') h + (' num2str(P.p1) ') h^2'],...
'\rho(M^{-1}N) exacte')
title('Jacobi')
set(gca,'FontSize',20,'linewidth',2)
%%
figure(2)
subplot(1,2,1)
loglog(1./(nn+1),1-pest,'b*',h,(pi^2/2)*h.^2,'k-',...
'MarkerSize',10,'Linewidth',2)
xlabel('h')
ylabel('1-\rho')
legend('1-\rho (calcule)','1-\rho (exacte)')
title('Jacobi')
set(gca,'Fontsize',20,'Linewidth',2)
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);
% Gauss-Seidel
nn = 10:5:100;
pest = zeros(length(nn),1);
p = 1;
for k=1:length(nn)
n = nn(k);
x = linspace(0,1,n+2);
u = sin(pi*x')*sin(pi*x) + 1;
[uu,pest(k)] = Heat2D(u,p);
end
semilogy(1:vcol,error)
P = fit(1./(nn'+1),pest,'poly2');
h = 1./(linspace(nn(1),nn(end),1000)+1);
figure(1)
subplot(1,2,2)
plot(1./(nn+1),pest,'b*',h,P(h),'r--',h,cos(pi*h).^(p+1),'k-',...
'MarkerSize',10,'Linewidth',2)
xlabel('h')
ylabel('\rho(M^{-1}N)')
legend('Taux de convergence asymptotique',...
[num2str(P.p3) ' + (' num2str(P.p2) ') h + (' num2str(P.p1) ') h^2'],...
'\rho(M^{-1}N) exacte')
title('Gauss-Seidel')
set(gca,'FontSize',20,'linewidth',2)
figure(2)
subplot(1,2,2)
loglog(1./(nn+1),1-pest,'b*',h,(pi^2)*h.^2,'k-',...
'MarkerSize',10,'Linewidth',2)
xlabel('h')
ylabel('1-\rho')
legend('1-\rho (calcule)','1-\rho (exacte)')
title('Gauss-Seidel')
set(gca,'Fontsize',20,'Linewidth',2)
%% Laplacian convergence
n = 7;
n = 20;
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));
[v,vk] = SteepestDescent(A,f,ones(n^2,1),1e-10,100);
[vrow,vcol] = size(vk);
error = zeros(vcol,1);
for i=1:vcol
......@@ -47,9 +90,46 @@ for i=1:vcol
error(i) = sqrt(w'*A*w);
end
P = fit((1:vcol)',log(error),'poly1');
plot(1:vcol,log(error),'.-',1:vcol,P(1:vcol),'--')
ind = round(0.8*vcol):vcol;
P0= fit((1:vcol)',log(error),'poly1');
P1= fit(ind',log(error(ind)),'poly1');
C0 = exp(P0.p1);
C1 = exp(P1.p1);
Cex= cos(pi/(n+1));
ind = (1-vcol):0;
plot(1:vcol,log(error),'b.',...
1:vcol,log(error(end) * C0 .^ind),'m--',...
1:vcol,log(error(end) * C1 .^ind),'r--',...
1:vcol,log(error(end) * Cex.^ind),'k-')
xlabel('Iteration')
ylabel('log(Error)')
legend('|e_k|_A','Pente (entiere)','Pente (partielle)','Pente (exacte)')
%% Laplacian convergence - fonction de h
nn= 10:5:100;
C = zeros(1,length(nn));
for k = 1:length(nn)
n = nn(k);
A = FDLaplacianTrue(n,2); A = -A;
x = linspace(0,1,n);
u = cos(x') * sin(x); u = u(:);
f = A*u;
[v,vk] = SteepestDescent(A,f,ones(n^2,1),1e-10,1000);
[vrow,vcol] = size(vk);
error = zeros(vcol,1);
for i=1:vcol
w = vk(:,i) - u;
error(i) = sqrt(w'*A*w);
end
ind = round(0.8*vcol):vcol;
P = fit(ind',log(error(ind)),'poly1');
C(k)= exp(P.p1);
end
Cest = exp(P.p1); %Cest = (1+Cest)/(1-Cest);
\ No newline at end of file
h = 1./(nn+1);
plot(h,C,'b*',h,cos(pi*h),'k--')
xlabel('h')
ylabel('$\frac{\kappa(A)-1}{\kappa(A)+1}$','interpreter','latex')
legend('Approx.','Exacte')
\ No newline at end of file
\documentclass[a4paper, 12pt]{article}
\usepackage{preamble}
\usepackage[french]{babel}
\usepackage[utf8]{inputenc}
\usepackage{tikz}
\usepackage[margin=2cm]{geometry}
\renewcommand{\vec}[1]{\mathbf{#1}}
\title{S\'erie 3, Exercices 1c) et 2c)}
\date{}
\begin{document}
\maketitle
\section{Exercice 1c)}
\label{sec:ex1}
\newcommand{\diag}{\text{diag}}
Il est possible d'écrire les $\rho(M^{-1}N)$ pour Jacobi (que je note comme $\rho_J$) avec les m\'ethodes de l'exercice 1a).
Consid\'erer $A$, la matrice du Laplacien en dimension deux, sans les \'el\'ements sur la diagonale.
On sait que $A = I \otimes D + D \otimes I$ par les s\'eries pr\'ec\'edentes, o\`u $D$ est la matrice du Laplacien en dimension un.
Donc:
\begin{equation*}
\diag(A) = I \otimes \diag(D) + \diag(D) \otimes I \implies \ N = I \otimes (D - \diag(D)) + (D - \diag(D)) \otimes I,
\end{equation*}
et alors les valeurs propres de $N$ sont $-\frac{2}{h^2} \left ( \cos \left ( i \pi h \right ) + \cos \left ( j \pi h \right ) \right )$.
La matrice $M^{-1}$ n'est qu'une valeur scalaire $-\frac{h^2}{4}$ fois l'identit\'e et la plus grande valeur propre de $N$ en valeur absolue est quand $i=j=1$.
Alors, $$\rho_J = \cos( \pi h) = 1 - \frac{\pi^2}{2} h^2 + \order{h^4}.$$
De plus, pour les matrices avec une propriet\'e specifique, $\rho_{GS} = \rho_J^2$ (o\`u $\rho_{GS} = \rho(M^{-1}N)$ pour Gauss-Seidel).
Le Laplacien discret a cette propriet\'e.
Alors, $$\rho_{GS} = \cos(\pi h)^2 = 1 - \pi^2 h^2 + \order{h^4}.$$
Mais comment les approximer num\'eriquement?
D'abord, on doit trouver une relation entre les objets d'it\'eration qu'on peut calculer et les $\rho_J$ et $\rho_{GS}$.
Par exemple, $$R_\infty(G) = \ln(\rho(G)). \quad \text{(eqn. 2.16 dans le polycopi\'e)}$$
$R_\infty(G)$ est le taux de convergence asymptotique.
\c{C}a veut dire qu'il est la pente $m$ de la ligne $$\ln\norm{e_k} = m k + \ln\norm{e_0},$$ quand $k$ tend vers l'infini.
(Peut-\^etre cette d\'efinition est un peu bizarre.
J'explique plus.
On sait que $\norm{e_k} \leq \norm{G^k} \norm{e_0}$ (la norme de la forme d'erreur, eqn. 2.6).
Si on prend le logarithme de cette relation on a
\begin{align*}
\ln \norm{e_k} & \leq \ln \norm{G^k} + \ln \norm{e_0} \\
& \leq k \ln \norm{G^k}^{1/k} + \ln \norm{e_0}. && \text{(une astuce du logarithme)}
\end{align*}
Par la s\'erie 1, quand $k$ tend vers l'infini, $\norm{G^k}^{1/k}$ tend vers $\rho(G)$.
Alors, la pente de cette ligne, qu'on note par $R_\infty$, devient $\ln \rho(G)$.
)
Donc, pour estimer la valeur de $\rho(G)$ on peut construire cette ligne.
Il est n\'ecessaire de calculer $e_k$, les erreurs exactes, pour beaucoup de $k$.
On consid\`ere une solution $u$, calcule une fonction $f$ par $f=Au$, et r\'eso\^ut le syst\`eme $Au = f$ par les m\'ethodes Jacobi et Gauss-Seidel.
\begin{figure}[h!]
\centering
\includegraphics[width=0.5\textwidth]{TauxAsymp.png}
\caption{Jacobi, $n=10$.}
\label{fig:taux}
\end{figure}
Regarder la figure \ref{fig:taux}.
Les erreurs $\ln \norm{e_k}$ cr\'eent une ligne.
La pente de cette ligne est -0.04135.
On peut la calculer par $$\frac{\ln \norm{e_{100}} - \ln \norm{e_0}}{100} = \frac{1}{100} \ln \left ( \frac{\norm{e_{100}}}{\norm{e_0}} \right )$$ ou r\'esoudre un syst\`eme d'\'equations normales:
\begin{equation*}
\begin{bmatrix} 0 & 1 \\ 1 & 1 \\ \vdots & \vdots \\ 100 & 1 \end{bmatrix} \begin{bmatrix} m \\ c \end{bmatrix} = \begin{bmatrix} \ln \norm{e_0} \\ \ln \norm{e_1} \\ \vdots \\ \ln \norm{e_{100}} \end{bmatrix} .
\end{equation*}
La derni\`ere m\'ethode est appell\'e 'a line of best fit'.
La valeur $m$ approxime $R_\infty$ et $c$ approxime $\norm{e_0}$.
\begin{figure}[h!]
\centering
\includegraphics[width=\textwidth]{ex1c.png}
\caption{$\rho_J$ et $\rho_{GS}$ comme fonctions en $h$.}
\label{fig:ex1c}
\end{figure}
On r\'ep\`ete pour plusieurs valeurs de $h$.
Apr\`es un certain nombre de $h$ on peut construire une figure comme figure \ref{fig:ex1c}.
Il y a une comparaison entre les taux de convergence asymptotiques (qui approximent $\rho(G)$) et $\rho(G)$ exacte.
Il y a aussi un autre 'fit', de polyn\^ome de d\'egree deux.
Si on ne veut pas utiliser un 'fit' on peut calculer la pente et la constante de la ligne
\begin{equation*}
\ln (1 - \rho(G)) = \ln \alpha + 2 \ln h.
\end{equation*}
\c{C}a veut dire que si on a les $y = \ln (1 - \rho(G))$ et les $x = \ln h$ alors on devrait avoir une ligne $y=mx + c$ avec $m = 2$ et $c = \ln \alpha$.
On peut aussi estimer $\alpha$ en prenant $(1 - \rho(G))/h^2$ pour un $h$ quelconque.
Par exemple, pour la m\'ethode de Jacobi et $n=10$ on sait que $\rho \approx \exp(-0.04135)$.
Alors, $$\alpha \approx (1-0.9595)(10+1)^2 = 4.9013.$$
La valeur exacte est $\pi^2/2 = 4.9348$, et nous avons une bonne estimation.
\begin{figure}
\centering
\includegraphics[width=\textwidth]{ex1c_alt.png}
\caption{Pr\'esentation alternative des r\'esultats, telle qu'ils forment des lignes droites.}
\label{fig:ex1c droites}
\end{figure}
\newpage
\section{Exercice 2c)}
\label{sec:ex2}
Par exercice 1a) on sait que $$\kappa(A) = \frac{\lambda_{max}(A)}{\lambda_{min}(A)} = \frac{1 + \cos(\pi h)}{1 - \cos(\pi h)}.$$
Donc $$\frac{\kappa(A) - 1}{\kappa(A) + 1} = \frac{1 + \cos(\pi h) - (1 - \cos(\pi h))}{1 + \cos(\pi h) + (1 - \cos(\pi h))} = \cos(\pi h).$$
Alors, la valeur que nous devons estimer change avec $n$.
On choisit un $n$ sp\'ecifique (par exemple, $n=20$) et proc\'eder.
Par Th\'eor\`eme 19 dans chapitre 3 du polycopi\'e on sait que $$\norm{e_k}_A \leq \left ( \frac{\kappa(A) - 1}{\kappa(A) + 1} \right )^k \norm{e_0}_A = C^k \norm{e_0}_A$$ pour la m\'ethode de descente maximale.
Rappeler que $\norm{\vec{v}}_A = \sqrt{\vec{v}^\top A \vec{v}}$.
Alors, nous pouvons r\'ep\'eter l'id\'ee que nous avons utilis\'e pour l'exercice 1c):
Trouver la pente de la ligne $$\ln \norm{e_k}_A = mk + \ln \norm{e_0}_A.$$
On a besoin d'une solution exacte.
Je choisis $u(x,y) = \cos(y) \sin(x)$ et je construis le vecteur $\vec{u}_{ex} \in \bbr^{n^2}$ qui est la fonction $u(x,y)$ \'evalu\'ee sur les points int\'erieures et discr\`etes.