with(Student:-NumericalAnalysis):
with(LinearAlgebra):
with(ArrayTools):
with(plots):# Question 1.1# A est suppos\303\251e inversible et triangulaire inf\303\251rieure
ResolInf := proc (A, b)
local n, x, i;
n := Size(A,1);
x := Vector (n);
for i from 1 to n do
x[i] := (b[i] - add(A[i,j]*x[j], j=1..(i-1)))/A[i,i];
od;
x;
end proc;AL := Matrix([[2,0,0],[4,5,0],[7,8,9]]);
b := <1,2,3>;
sL := ResolInf (AL, b);
evalm(AL &* sL);# Question 1.2# A est suppos\303\251e inversible et triangulaire sup\303\251rieure
ResolSup := proc (A, b)
local n, x, i;
n := Size(A,1);
x := Vector(n);
for i from n to 1 by -1 do
x[i] := (b[i] - add(A[i,j]*x[j], j=(i+1)..n))/A[i,i];
od;
x;
end proc;AU := Matrix([[2,5,4],[0,5,1],[0,0,9]]);
b := <1,2,3>;
sU := ResolSup (AU, b);
evalm(AU &* sU);# Question 1.3SolveLUP := proc (L,U,P,b)
local x, y, Pb;
Pb := evalm(P &* b);
y := ResolInf (L, Pb);
x := ResolSup (U, y);
return x;
end proc;A := Matrix([[7,9,25],[2,4,6],[0,7,3]]);
AP, AL, AU := MatrixDecomposition (A, method=PLU, output=[P,L,U]);
sA := SolveLUP(AL, AU, AP, <1,2,3>);
evalm(A &* sA);# Question 1.4CombineLignes := proc (A, dest, lambda, src)
local n, j;
n := Size(A,1);
for j from 1 to n do
A[dest,j] := A[dest,j] + lambda * A[src,j];
od;
return A;
end proc;A;
CombineLignes(A, 3, -5, 2);
A;SwapLignes := proc (A, i1, i2)
local j, tmp, n;
n := Size(A,1);
for j from 1 to n do
tmp := A[i1,j];
A[i1,j] := A[i2,j];
A[i2,j] := tmp;
od;
return A;
end proc;A;
SwapLignes(A, 3, 2);
A;TrouverLUP := proc (A)
local i, j, n, idest, L, U, P, lambda;
n := Size(A,1);
L := Matrix(n);
U := Matrix(n);
P := Matrix(n);
# Initialisation des matrices
# (pourrait \303\252tre fait plus succintement avec des fonctions)
for i from 1 to n do
for j from 1 to n do;
if (i = j) then
L[i,j] := 0;
U[i,j] := A[i,j];
P[i,j] := 1;
else
L[i,j] := 0;
U[i,j] := A[i,j];
P[i,j] := 0;
fi;
od;
od;
# Echange les lignes si n\303\251cessaire (0 en haut \303\240 gauche)
# Pour la stabilit\303\251 num\303\251rique, il faudrait \303\251changer avec la plus grande valeur
for i from 1 to n do
for j from i+1 to n do
if (U[i,i] = 0 and U[j,i] <> 0) then
SwapLignes (U, i, j);
SwapLignes (P, i, j);
# Il faut aussi \303\251changer les lignes dans L parce qu'au final on travaille sur PA
SwapLignes (L, i, j);
fi;
od;
for idest from i+1 to n do
lambda := -U[idest,i] / U[i,i];
L[idest, i] := -lambda; # Stocke le facteur de la transvection au bon endroit dans L
CombineLignes (U, idest, lambda, i);
od;
od;
for i from 1 to n do
L[i,i] := 1;
od;
return (P, L, U);
end proc:#A := Matrix([[7,9,25,84],[6,3,4,17],[7,9,4,16],[2,4,6,2]]);
A := Matrix([[7,9,25,84],[7,9,4,16],[6,3,4,17],[2,4,6,2]]);
P, L, U := TrouverLUP(A);
MatrixDecomposition (A, method=PLU);evalm(P &* A);
evalm(L &* U);# Question 1.5SolveLin := proc (A, b)
local L, U, P;
P, L, U := TrouverLUP(A);
return SolveLUP(L, U, P, b);
end proc:A := Matrix([[7,9,25,84],[7,9,4,16],[6,3,4,17],[2,4,6,2]]);
b := <1,2,3,4>;
x := SolveLin (A, b);
evalm (A &* x);# Question 2.1DetParLUP := proc (A)
local L, U, P, n;
n := Size(A,1);
P, L, U := TrouverLUP(A);
print (U);
# Ici je triche pour avoir la signature de la permutation P
# Mais il est facile d'obtenir la signature en calculant
# une d\303\251composition en transpositions de la signature
# (ajouter un tableau sigma dans TrouverLUP o\303\271 sigma[i] = sigma(i)
# qui permet de faire le calcul facilement)
return Determinant(P) * mul(U[i,i], i=1..n);
end proc:A := Matrix([[7,9,25,84],[7,9,4,16],[6,3,4,17],[2,4,6,2]]);
DetParLUP(A);
Determinant(A);# Question 2.2InverseMat := proc (A)
local L, U, P, n, b, Ainv, i;
P, L, U := TrouverLUP(A);
n := Size(A,1);
Ainv := Array(1..n);
for i from 1 to n do
b := Vector(n, 0);
b[i] := 1;
Ainv[i] := SolveLUP(L,U,P, b);
od;
return Matrix(n, (i,j) -> Ainv[j][i]);
end proc:A := Matrix([[7,9,25,84],[7,9,4,16],[6,3,4,17],[2,4,6,2]]);
InverseMat (A);
MatrixInverse(A);# Question 2.3Interpol := proc (fs, x, y)
# fs est un vecteur de fonctions f_i de base de l'espace
# x et y sont des vecteurs avec pour \303\251l\303\251ments les x_i et y_i
local L, U, P, A, At, c, m, n;
n := Size(fs,1);
m := Size(x,1);
A := Matrix (m, n, (i,j) -> fs[j](x[i]));
At := Transpose(A);
c := SolveLin(evalm(At &* A), evalm(At &* y));
return c;
end proc:n := 5;
m := 50;
unassign('x');
# Base et donn\303\251es
base := map (e -> unapply(e,x), <seq(x^i, i = 1..n)>):
roll := rand(1..50):
# Abscisses
x := < seq(i/m + roll()/50 * (-1)^(rand(0..1)()), i=1..m) >:
xmax := max(x):
xmin := min(x):
# Fonction th\303\251orique
cth := <seq(0.1*roll()/50 * (-1)^(rand(0..1)()), i=1..n)>:
f := x -> add(cth[i] * base[i](x), i=1..n):
# Donn\303\251es mesur\303\251es (avec bruit)
y := < seq(f(x[i]) + roll()/400 * (-1)^(rand(0..1)()), i=1..m) >:
# Donn\303\251es
data := plot(x, y, style = point, symbol = solidbox, symbolsize=15, color = blue):
# Trouver le Least Square approximant
c := evalf(Interpol(base, x, y)):
unassign('x');
# Affiche le r\303\251sultat
interpolantplot := plot (add(c[i] * base[i](x), i=1..n), x=xmin..xmax):
# Th\303\251orique
thplot := plot (f(x), x=xmin..xmax, color=green):
#display({thplot, data});
#display({data});
display({data, interpolantplot, thplot});TTdSMApJNlJUQUJMRV9TQVZFLzE5NTc5NzU0OFgsJSlhbnl0aGluZ0c2IjYiW2dsISIlISEhIyoiJCIkIiIjIiIlIiIoIiIhIiImIiIpRipGKiIiKkYmTTdSMApJNlJUQUJMRV9TQVZFLzE2MzI5NDE5MlgqJSlhbnl0aGluZ0c2IjYiW2dsISMlISEhIiQiJCIiIiIiIyIiJEYmTTdSMApJNlJUQUJMRV9TQVZFLzE2MzI5NDAzMlgqJSlhbnl0aGluZ0c2IjYiW2dsISMlISEhIiQiJCMiIiIiIiMiIiEjISIiIiM9RiY=TTdSMApJNlJUQUJMRV9TQVZFLzE2MzI5MjMzNlgsJSlhbnl0aGluZ0c2IjYiW2dsISIlISEhIyoiJCIkIiIjIiIhRigiIiZGKUYoIiIlIiIiIiIqRiY=TTdSMApJNlJUQUJMRV9TQVZFLzE2MzI5MTIxNlgqJSlhbnl0aGluZ0c2IjYiW2dsISMlISEhIiQiJCIiIiIiIyIiJEYmTTdSMApJNlJUQUJMRV9TQVZFLzE2MzI4Mzc5MlgqJSlhbnl0aGluZ0c2IjYiW2dsISMlISEhIiQiJCEiIiMiIiIiIiRGKEYmTTdSMApJNlJUQUJMRV9TQVZFLzE2MzI4MzY5NlgsJSlhbnl0aGluZ0c2IjYiW2dsISIlISEhIyoiJCIkIiIoIiNfISM1IiIqIiNwISM4IiNEIiRUIiEjRkYmTTdSMApJNlJUQUJMRV9TQVZFLzE4NDM2NzI1MlgsJSlhbnl0aGluZ0c2IjYiW2dsISIlISEhIyoiJCIkIiIiIiIhRihGKEYoRidGKEYnRihGJg==TTdSMApJNlJUQUJMRV9TQVZFLzE2MzYwNjAwMFgsJSlhbnl0aGluZ0c2IjYiW2dsISIlISEhIyoiJCIkIiIiIiIhIyIiIyIiKEYoRicjIiM1IiNcRihGKEYnRiY=TTdSMApJNlJUQUJMRV9TQVZFLzE2MzYwNjcwNFgsJSlhbnl0aGluZ0c2IjYiW2dsISIlISEhIyoiJCIkIiIoIiIhRigiIipGJ0YoIiNEIiIkIyEjJykiI1xGJg==TTdSMApJNlJUQUJMRV9TQVZFLzE4NDM1OTczMlgqJSlhbnl0aGluZ0c2IjYiW2dsISMlISEhIiQiJCMiI2siI1YjIiNJRikjISNGRilGJg==TTdSMApJNlJUQUJMRV9TQVZFLzE2MzI4MzY5NlgsJSlhbnl0aGluZ0c2IjYiW2dsISIlISEhIyoiJCIkIiIoIiNfISM1IiIqIiNwISM4IiNEIiRUIiEjRkYmTTdSMApJNlJUQUJMRV9TQVZFLzE2MzI4MzY5NlgsJSlhbnl0aGluZ0c2IjYiW2dsISIlISEhIyoiJCIkIiIoIiNfISM1IiIqIiNwISM4IiNEIiRUIiEjRkYmTTdSMApJNlJUQUJMRV9TQVZFLzE2MzI4MzY5NlgsJSlhbnl0aGluZ0c2IjYiW2dsISIlISEhIyoiJCIkIiIoIiNfISM1IiIqIiNwISM4IiNEIiRUIiEjRkYmTTdSMApJNlJUQUJMRV9TQVZFLzE2MzI4MzY5NlgsJSlhbnl0aGluZ0c2IjYiW2dsISIlISEhIyoiJCIkIiIoIiNfISM1IiIqIiNwISM4IiNEIiRUIiEjRkYmTTdSMApJNlJUQUJMRV9TQVZFLzE2MzI4MzY5NlgsJSlhbnl0aGluZ0c2IjYiW2dsISIlISEhIyoiJCIkIiIoIiNfISM1IiIqIiNwISM4IiNEIiRUIiEjRkYmTTdSMApJNlJUQUJMRV9TQVZFLzE2MzI4MzY5NlgsJSlhbnl0aGluZ0c2IjYiW2dsISIlISEhIyoiJCIkIiIoIiNfISM1IiIqIiNwISM4IiNEIiRUIiEjRkYmTTdSMApJNlJUQUJMRV9TQVZFLzIwNzUyNDg4MFgsJSlhbnl0aGluZ0c2IjYiW2dsISIlISEhIzEiJSIlIiIoRiciIiciIiMiIipGKiIiJCIiJSIjREYsRixGKCIjJSkiIzsiIzxGKUYmTTdSMApJNlJUQUJMRV9TQVZFLzIwOTcxMjcwNFgsJSlhbnl0aGluZ0c2IjYiW2dsISIlISEhIzEiJSIlIiIiIiIhRihGKEYoRihGJ0YoRihGJ0YoRihGKEYoRihGJ0YmTTdSMApJNlJUQUJMRV9TQVZFLzIwNzUyNTEzMlgsJSlhbnl0aGluZ0c2IjYiW2dsISIlISEhIzEiJSIlIiIiIyIiJyIiKEYnIyIiI0YqIiIhRidGLSMhIzUiI0xGLUYtRicjIiQ3IyIkJHBGLUYtRi1GJ0YmTTdSMApJNlJUQUJMRV9TQVZFLzIwNTI1NDA0MFgsJSlhbnl0aGluZ0c2IjYiW2dsISIlISEhIzEiJSIlIiIoIiIhRihGKCIiKiMhI0xGJ0YoRigiI0QjISRBIkYnISNARigiIyUpISNiISNvIyEmIVE3IiQkcEYmTTdSMApJNlJUQUJMRV9TQVZFLzE4MzkzMDcxMlgsJSlhbnl0aGluZ0c2IjYiW2dsISIlISEhIzEiJSIlIiIiIiIhRihGKEYoRihGJ0YoRihGJ0YoRihGKEYoRihGJ0YmTTdSMApJNlJUQUJMRV9TQVZFLzIwNDc0MjEyOFgsJSlhbnl0aGluZ0c2IjYiW2dsISIlISEhIzEiJSIlIiIiIyIiJyIiKEYnIyIiI0YqIiIhRidGLSMhIzUiI0xGLUYtRicjIiQ3IyIkJHBGLUYtRi1GJ0YmTTdSMApJNlJUQUJMRV9TQVZFLzIwNTI1NDI5MlgsJSlhbnl0aGluZ0c2IjYiW2dsISIlISEhIzEiJSIlIiIoIiIhRihGKCIiKiMhI0xGJ0YoRigiI0QjISRBIkYnISNARigiIyUpISNiISNvIyEmIVE3IiQkcEYmTTdSMApJNlJUQUJMRV9TQVZFLzE4Mzk5MTMzMlgsJSlhbnl0aGluZ0c2IjYiW2dsISIlISEhIzEiJSIlIiIoRiciIiciIiMiIipGKiIiJCIiJSIjREYsRixGKCIjJSkiIzsiIzxGKUYmTTdSMApJNlJUQUJMRV9TQVZFLzE4Mzk5MTI2MFgqJSlhbnl0aGluZ0c2IjYiW2dsISMlISEhIiUiJSIiIiIiIyIiJCIiJUYmTTdSMApJNlJUQUJMRV9TQVZFLzE4NDE4NDU4OFgqJSlhbnl0aGluZ0c2IjYiW2dsISMlISEhIiUiJSMiJWVFIiUmNCQjISV3NUYpIyIlSEBGKSMhJC4oRilGJg==TTdSMApJNlJUQUJMRV9TQVZFLzE5NjcxNTY0NFgsJSlhbnl0aGluZ0c2IjYiW2dsISIlISEhIzEiJSIlIiIoRiciIiciIiMiIipGKiIiJCIiJSIjREYsRixGKCIjJSkiIzsiIzxGKUYmTTdSMApJNlJUQUJMRV9TQVZFLzIwOTU4NDEyOFgsJSlhbnl0aGluZ0c2IjYiW2dsISIlISEhIzEiJSIlIiIoIiIhRihGKCIiKiMhI0xGJ0YoRigiI0QjISRBIkYnISNARigiIyUpISNiISNvIyEmIVE3IiQkcEYmTTdSMApJNlJUQUJMRV9TQVZFLzE4NTAxNjQ4NFgsJSlhbnl0aGluZ0c2IjYiW2dsISIlISEhIzEiJSIlIiIoRiciIiciIiMiIipGKiIiJCIiJSIjREYsRixGKCIjJSkiIzsiIzxGKUYmTTdSMApJNlJUQUJMRV9TQVZFLzIwOTk4MzcyMFgsJSlhbnl0aGluZ0c2IjYiW2dsISIlISEhIzEiJSIlIyEkJEciJSE+JyMiI2hGKSMiIzgiJSY0JCMiIyQpRikjISQ9I0YuIyIkaCZGLiMhJD4kRi4jIiNgRi4jIiRiJCIlUTcjISRwI0Y7IyIjTSIkPicjISNARjsjIiR0JiImIVE3IyEkNiNGRUYzIyEkJHBGRUYmTTdSMApJNlJUQUJMRV9TQVZFLzIwNTQ2MzY0NFgsJSlhbnl0aGluZ0c2IjYiW2dsISIlISEhIzEiJSIlIyEkJEciJSE+JyMiI2hGKSMiIzgiJSY0JCMiIyQpRikjISQ9I0YuIyIkaCZGLiMhJD4kRi4jIiNgRi4jIiRiJCIlUTcjISRwI0Y7IyIjTSIkPicjISNARjsjIiR0JiImIVE3IyEkNiNGRUYzIyEkJHBGRUYm