Commit 66a9cbf0 authored by Conor McCoid's avatar Conor McCoid

Updates to MI course

parent 3c3e075b
function [u,F]=Heat2D(f,g,p)
function u=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
......@@ -12,22 +12,25 @@ 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);
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 );
N = triu(D,1) + tril(D,-1);
M = 1 / (4 * (n+1)^2);
elseif p==1
N = triu(D,1);
M = D - N; M = -M \ kron(I,I);
end
u0 = kron(zeros(n,1),zeros(n,1));
u0 = kron(ones(n,1),ones(n,1));
err= 1;
tol= 1e-8;
tol= eps;
itermax = 1000;
iter = 1;
while err > tol && iter < itermax
u = D*LU*u0 - D*F;
err = norm(u-u0);
u0= u;
u = M*N*u0 - M*F;
u0 = u;
err= norm(D*u - F);
iter = iter+1;
end
\ 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
x = linspace(0,1,n+2);
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);
% g = u(:); g = g(II==1);
% f = -2*pi^2 * ( u(2:end-1,2:end-1) - 1 ); f = f(:);
x = x(2:end-1);
g = ones(4*n,1);
f = -2*pi^2 * sin(pi*x')*sin(pi*x); f = f(:);
uu = Heat2D(f,g,1);
uu = reshape(uu,n,n);
surf(x,x',uu - u(2:end-1,2:end-1))
\ No newline at end of file
......@@ -48,15 +48,24 @@
\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}
\item (4 pts) Consid\'erer le Laplacien discret $A \in \R^{n^2 \times n^2}$ en dimension deux pour le carr\'e unitaire avec $n^2$ points intérieures (alors, $h = \frac{1}{n+1}$).
Soit $u = \sin(\pi x) \sin(\pi y) + 1$.
R\'esoudre le syst\`eme $Au = f$ avec la fonction \texttt{Heat2D} de la s\'erie pr\'ec\'edent, en utilisant les méthodes de Jacobi et Gauss-Seidel.
Utiliser les erreurs \`a estimer $\rho(M^{-1}N)$ comme une fonction de $h$ pour les deux m\'ethodes.
\textit{Indice:} Modifier \texttt{Heat2D} telle qu'elle compute l'erreur et le taux de convergence asymptotique.
Pour le dernier, c'est le plus facile \`a computer la pente d'une ligne.
Les $\rho(M^{-1}N)$ ont la forme $1 - \alpha h^2 + \mathcal{O}(h^4)$; $\alpha=?$
% \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}
......
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
h=1/(n+1);
D=1/h^2*spdiags(ones(n,1)*[1 -2 1],[-1 0 1],n,n);
if d==1
A=D;
elseif d==2
I = eye(n);
A=kron(I,D) + kron(D,I);
elseif d==3
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
%if d==1
% A=D;
%elseif d==2
% A=[];
% for i=1:n
% A=blkdiag(A,D);
% end
% I2=1/h^2*spdiags(ones(n^2,1)*[1 1],[-n n],n^2,n^2);
% A=A+I2;
%elseif d==3
% C=[];
% for i=1:n
% C=blkdiag(C,D);
% end
% A=[];
% for i=1:n
% A=blkdiag(A,C);
% end
% 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=A+I2+I3;
%end
function [u,pest]=Heat2D(uu,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,m] = size(uu);
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);
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;
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);
elseif p==1
N = triu(D,1);
M = D - N; M = -M \ kron(I,I);
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;
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
% Test suite for Serie 3
%% Heat2D, estimate spectral radius
nn = 10:5:50;
pest = zeros(length(nn),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,0);
end
% plot(log(1./(nn+1)),log(pest-1),'r*--')
P = fit(1./(nn'+1),pest,'poly2');
%%
n = 10;
A = rand(n); A = -(A+A')/2;
A = diag(-sum(A,2)+0.1) + A;
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment