Cours de l'IMT Atlantique sur les Méthodes Numériques

A100_Octave.md

Méthodes Numériques avec Octave

A101_Octave.md

Octave

clear all % Efface les variables en mémoire
close all % Fermeture de toutes les figures
tic, toc % chronometrage

Variables

val = 15; # ( ; disable the console output )
disp(val) # print the value

log(val) % logarithme naturel
exp(val) % e puissance val
abs(val) % valeur absolue
floor(val) % arrondi inferieur
ceil(val) % arrondu superieur

Vecteurs et Matrices

x = [1; 2; 3]; % vecteur colonne
y = [4 5 6]; % vecteur ligne

x' % transposée de x. Dans le cas d'une matrice carrée, on inverse le triangle supérieur et inférieur!

x + y % broadcasting pour réaliser l'op => les valeurs sont dupliquées

x * y % multiplication de matrices. il faut des tailles compatibles

x .* y % même résultat que précédent, mais en fait c'est elem par elem

max(x) % maximum
min(x) % minimum
sum(x) % somme
prod(x) % produit

A = [4  2  1; 
    4 -5  1; 
    2 -1  6] % matrice carrée

A(5); % 5ème valeur indiquée dans A
A(1,:); % première ligne de A
A(:,2); % 2eme colonne de A
A(2:end, 1:2); % sous matrice de ligne 2 à end et de colonne 1 à 2

A(:,2)=A(:,3); % affecter la 3ème colonne de A à la 2ème colonne de A

B = [A(:,1) A(:,3)]; % concaténer 2 vecteurs colonne

ones(4); % matrice carrée de taille 4 avec que des 1
eye(4); % idem que one mais la diagonale contient que des 1
diag(ones(4)); % prend la diagonale d'une matrice et en fait un vecteur colonne
[ones(4)-eye(4) zeros(4,2) ones(4,2)*2]; % fonction préféfinies

Fonctions et Graphiques

clc; % reset l'affichage console (et pas la mémoire)
clear all; % reset la mémoire (et pas l'affichage console)

function vecteur = produit(x)
    A = [4 2 1; 1 -5 1;2 -1 6]; % fonction
    vecteur = A*x;
endfunction
b = produit(2); % appel fonction

function v = fibonacci(n)
    v = zeros(n,1); % initialiser le vecteur colonne
    v(1) = 1; v(2) = 1;
    for i = 3:n
      v(i) = v(i-1)+v(i-2);
    endfor
endfunction

x = fibonacci(5)';

x = 0:2:10; % discretisation de l'intervale [0,10] avec un pas de 2
n = length(x); % récupérer la longueur d'un vecteur
y = zeros(n); % y doit être initialiser à la même longueur que x
for i = 1:n
   y(i) = exp(x(i))+1; % affectation de y
endfor

plot(x,y); 

B100_Resolution_systeme_lineaire.md

1. Resolution des systèmes linéaires

A = [1  -5   1; 
     2  -1   6; 
     4   2  -3 ]

b = [ 24; 
      37; 
      10 ]

linsolve(A, b) # Automatisation

B111_EliminationGauss.md

1.1 Elimination de Gauss (Méthode directe)

• Etape 1) Triangulisation de la matrice carrée A avec b
• Etape 2) Substitution rétrograde pour l'obtention des valeurs de x

B112_EliminationGauss.m

1;

function x = EliminationGauss(A, b)
[A, b] = triangulisation(A, b)
x = substitution(A, b);
endfunction

# Transformation en matrice triangulaire superieure
function [A, b] = triangulisation(A, b)
  n = rows(b); # ordre de la matrice
  for k = 1:n-1
    if abs(A(k, k)) < 0.00001
      return; # Pas de triangulisation possible
    else
      for i = k+1:n
        m = A(i,k)/A(k,k);
        b(i) = b(i) - m*b(k);
        A(i, k) = 0;
        for j = k+1:n
          A(i, j) -= m*A(k, j);
        endfor
      endfor
    endif
  endfor
endfunction

# Substitution retrograde
function x = substitution(A, b)
n = rows(b); # ordre de la matrice
x = zeros(n, 1);
for i = 0:n-1
  j = n-i;
  x(j) = b(j);
  for k = j+1:n
   x(j)-=A(j,k)*x(k);
  endfor
  x(j)/=A(j,j);
endfor
endfunction

A = [1 1 1; 1 0.4 0; 10 5 1];
b = [24; 14; 152];
x = EliminationGauss(A, b)

#{
A =

   1.0000   1.0000   1.0000
        0  -0.6000  -1.0000
        0        0  -0.6667

b =

   24.0000
  -10.0000
   -4.6667

x =

   12.0000
    5.0000
    7.0000
#}

B121_MethodeJacobi.md

1.2 Méthode de Jacobi (Méthode itérative)

On sait que la méthode fonctionne si la diagonale de la matrice A est strictement dominante.

=> ne pas hésiter à modifier la matrice pour obtenir cette configuration

B122_MethodeJacobi.m

1;

function x = Jacobi(A, b, x0, epsilon)
  diff = epsilon + 1;
  x = x0; # valeur initiale 

  while (diff > epsilon)
    x_old = x;
    x = JacobiIter(A, b, x);
    diff = norm(x - x_old, inf);
  endwhile
endfunction

function x1 = JacobiIter(A, b, x)
  n = rows(b); # ordre de la matrice
  x1 = x;
  
  for i = 1:n
    x1(i) = b(i);
    for j = 1:n
      if j != i
        x1(i) -= A(i,j)*x(j);
      endif
    endfor
    x1(i) /= A(i,i);
  endfor
endfunction

# Ne converges pas
A = [1 1 1; 1 0.4 0; 10 5 1]
b = [24; 14; 152]
x = Jacobi(A, b, [1, 1, 1], 1e-1)

#{
A =

    1.0000    1.0000    1.0000
    1.0000    0.4000         0
   10.0000    5.0000    1.0000

b =

    24
    14
   152

x =

  -Inf   NaN  -Inf
#}

# Avec diagonale strictement dominante
A = [4 -1 0; -1 4 -1; 0 -1 4]
b = [6; 4; 6]
x = Jacobi(A, b, [1, 1, 1], 1e-3)

#{
A =

   4  -1   0
  -1   4  -1
   0  -1   4

b =

   6
   4
   6

x =

   1.9998   1.9998   1.9998
#}
  

B131_GaussSiedel.md

1.3 Méthode de Gauss-Seidel (Méthode itérative)

On sait que la méthode de Gauss-Seidel converge deux fois plus vite que celle de Jacobi si la matrice A est à diagonale strictement dominante

=> ne pas hésiter à modifier la matrice pour obtenir cette configuration

B132_GaussSiedel.m

1;

function x = GaussSeidel(A, b, x0, epsilon)
  diff = epsilon + 1;
  x = x0; # valeur initiale

  while (diff > epsilon)
    x_old = x;
    x = GaussSeidelIter(A, b, x);
    diff = norm(x - x_old, inf);
  endwhile
endfunction

function x1 = GaussSeidelIter(A, b, x)
  n = rows(b); # ordre de la matrice
  x1 = x;

  for i = 1:n
    x1(i) = b(i);
    for j = 1:n
      if (j > i)
        x1(i) -= A(i,j) * x(j);
      elseif (j < i)
        x1(i) -= A(i,j) * x1(j);
      endif
    endfor
    x1(i) /= A(i,i);
  endfor
endfunction

A = [1 1 1; 1 0.4 0; 10 5 1]
b = [24; 14; 152]
x = GaussSeidel(A, b, [1, 1, 1], 1e-1)

#{
A =

    1.0000    1.0000    1.0000
    1.0000    0.4000         0
   10.0000    5.0000    1.0000

b =

    24
    14
   152

x =

   12    5    7
#}

C111_Lagrange.md

2. Interpolation

L’interpolation consiste à construire une fonction fˆ(t) qui passe par les points (ti , yi ).

L’interpolation est une approximation (la seule méthode d’approximation que l’on va étudier dans ce cours, mais il y en a d’autres).

2.1 Lagrange

• Etape 1) On a un ensemble des noeuds par laquel la fonction doit passer:

 p1 = (x1, y1)
 p2 = (x2, y2)
     ...
 pn = (xn, yn)

• Etape 2) On obtient les équations de base du polynome:

 P(x1) = y1
 P(x2) = y2
     ...
 P(xn) = yn

• Etape 3) On obtient les équations du polynome de degré n:

ex (degré 3): Ax^3 + bX^2 + cX + d

• Etape 4) On remplace X par X1, X2, ..., Xn et on obtiens un système d'équation

• Etape 5) On résout grace à une des méthodes de la partie 1

Remarque:

On sait que pour n noeud, il existe:

• une infinité de polynomes de degré n (ou plus)
• un unique polynome de (de degré n-1 ou moins)

C112_Lagrange.m

1;

function L = Lagrange(x, y, X)
  n = columns(x);
  L = 0; 
  for i = 1:n
    L += PolyLagrange(x, y, i, X);
  endfor
endfunction

function Li = PolyLagrange(x, y, i, X)
  n = columns(x);
  Li = 1;
  for j = 1:n
    if (j != i)
      Li .*= (X - x(j));
    endif
  endfor
  for j = 1:n
    if (j != i)
      Li /= (x(i) - x(j));
    endif
  endfor
  Li *= y(i);
endfunction

x = [0 1 2 3];
y = [1 1.1 1.3 4];

plot(x, y, "o", "linewidth", 2); # points a interpoler

t = linspace(x(1), x(end), 1000); # echantilloner un interval
L = Lagrange(x, y, t); 

hold on;
plot(t, L); # polynome sur l'interval

C113_Lagrange.png

C114_Lagrange.png

C121_Spline.md

2.2 spline

En mathématiques appliquées et en analyse numérique, une spline est une fonction définie par morceaux par des polynômes.

C122_Spline.m

clear all % Efface les variables en mémoire
close all % Fermeture de toutes les figures
tic, toc % chronometrage

X = [ 0 1 2 3 ];
Y = [ 1 1.1 1.3 4 ];
hold on
plot(X, Y, "o");

function res = P(t)
  res = 0.4 * t.^3 - 1.15 * t.^2 + .85 * t + 1;
endfunction

t =linspace(X(1), X(end));
plot(t, P(t), "r");

function [A, b] = spline_mat(x,y)
  b = [diff(y) 0]';
  A = zeros(length(y));
  line = diff(x)./2;
  A = diag([line 1]) + diag(line,1);
endfunction

% je pense qu'on peut mieux faire
function y = spline_eval(x, y, m, t)
   xi_idx = interp1(x,1:length(x),t,'previous'); % voisin précédent
   xi_idxnext = interp1(x,1:length(x),t,'next'); % voisin suivant
   delta_xi = x(xi_idxnext)-x(xi_idx);
   mi = m(xi_idx);
   mi_next = m(xi_idxnext);
   y = (mi_next/(2*delta_xi))*(t-x(xi_idx))^2-(mi/(2*delta_xi))*(t-x(xi_idxnext))^2+y(xi_idx)+(delta_xi/2)*mi;
endfunction

[A, b] = spline_mat(X,Y);
A
b
m = linsolve(A,b) % peut être à remplacer par le pivot de Gauss....

t=linspace(X(1), X(end));

% @(t) est une fonction anonyme. arrayfun permet d'évaluer une fonction avec un tableau
y = arrayfun(@(t) spline_eval(X, Y, m, t), t);

plot(t, y, "b");

C123_Spline.png

gistfile2.txt