Soit
M un entier positif, soit
u ⃗ 0 ∈ ℝ M , soit
A la
M × M matrice définie par
A = ( M + 1 ) 2 ⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ 2 − 1 − 1 2 ⋱ − 1 ⋱ − 1 ⋱ 2 − 1 − 1 2 ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ .
On cherche
u ⃗ : ℝ + → ℝ M la solution de
u ⃗ ˙ ( t ) = − A u ⃗ ( t ) ,
t > 0 ,
u ⃗ ( 0 ) = u ⃗ 0 . On considère le schéma d'Euler rétrograde pour approcher
u ⃗ : soit
N un entier positif,
h le pas de temps,
t n = n h ,
n = 0 , 1 , 2 , . . . , N . On pose
u ⃗ 0 = u ⃗ 0 , pour
n = 0 , 1 , 2 , . . . , N − 1 , étant donné
u ⃗ n ,
u ⃗ n + 1 est tel que
u ⃗ n + 1 − u ⃗ n h = − A u ⃗ n + 1 .
On propose l'algorithme (octave, matlab) suivant, que l'on complètera.
Vous pouvez sauver le fichier
sysdiff.m , le compléter et le tester avec octave, puis répondre au questions ci-dessous.
function [u,err] = sysdiff(M,h,N)
%
% entrees : M : nombre de composantes de u
% h : pas de temps
% N : nombre de pas de temps
%
% sortie : u : approximation de u au temps final
% err : erreur entre u et son approximation au temps final
% condition initiale
for i=1:M
u(i)=sin(i*pi/(M+1));
end
% decomposition de cholesky de la matrice I+hA : I+hA = LL^T
diag = 1+2*h*(M+1)*(M+1);
sdiag = -h*(M+1)*(M+1);
a(1)=sqrt(diag);
for i=1:M-1
c(i) = sdiag/???;
a(i+1) = sqrt(diag-c(i)*???);
end
% schema d Euler retrograde : (I+hA) u^{n+1} = u^{n}
for n=1:N
% resolution du systeme lineaire Ly = u^{n}
u(1)=u(1)/a(1);
for i=1:M-1
u(i+1) = (u(i+1)-???*u(i))/a(i+1);
end
% resolution du systeme lineaire L^T u^{n+1} = y
u(M)=u(M)/a(M);
for i=M-1:-1:1
u(i) = (u(i)-???*u(i+1))/a(i);
end
end
% calcul de l erreur au temps final
err=0.;
coeff = exp(-(2-2*cos(pi/(M+1)))*(M+1)*(M+1)*n*h);
for i=1:M
err=err+(sin(i*pi/(M+1))*coeff-
u(i))^2;
end
err=sqrt(err);
end
Question 1
Entrez les ??? de la ligne (pas de blancs, pas de ^) :
c(i) = sdiag/???;
Answer for Question 1
Question 2
Entrez les ??? de la ligne (pas de blancs, pas de ^) :
a(i+1) = sqrt(diag-c(i)*???);
Answer for Question 2
Question 3
Entrez les ??? de la ligne (pas de blancs) :
u(i+1) = (u(i+1)-???*u(i))/a(i+1);
Answer for Question 3
Question 4
Entrez les ??? de la ligne (pas de blancs) :
u(i) = (u(i)-???*u(i+1))/a(i);
Answer for Question 4
Question 5
On veut tester l'ordre de convergence du schéma d'Euler rétrograde. On fixe M = 5 , le temps final t = 0.1 , on fait tendre h vers zéro et N vers l'infini. Que vaut l'erreur lorsque vous tapez
[u,err]=sysdiff(5,0.01,10)
[u,err]=sysdiff(5,0.001,100)
[u,err]=sysdiff(5,0.0001,1000)
Entrez trois réels séparés d'un blanc.
Answer for Question 5
Question 6
Soit λ une valeur propre de A et x ⃗ ≠ 0 le vecteur propre correspondant. On a :