Algorithme d’Analyse numérique sur octave

  Informatique

Factorisation LU

Est utilisé pour résoudre des systèmes de n équation et n inconue. Attention dans certains cas cet algorithme n’est pas stable, il faut utiliser a la place PA=LU

function [L U]=LU(A)
  %auteur: Nathan De Smedt
  [m n]=size(A);
  if m!=n
    "erreur"
  else
    for k=1:n
      for j=k:n
        U(k,j)=A(k,j);
      endfor
      L(k,k)=1;
      for i=k+1:n
        L(i,k)=A(i,k)/A(k,k);
      endfor
      for i=k+1:n
        for j=k+1:n
          A(i,j)=A(i,j)-L(i,k)*A(k,j);
        endfor
      endfor
    endfor
  endif

Factorisation QR

Algorithme de factorisation d’une matrice A en une matrice carrée (Q) et une triangulaire supérieure (L). Est utilisé pour résoudre des systèmes de n équation et m inconnue

function [Q R]=QR(A)
    %auteur: Nathan De Smedt
  [m n]= size(A);
  R=A;
  Q=eye(max(m,n));
  for k=1:n
    x=R(k:m,k);
    v=x;
    v(1)=v(1)+sign(x(1))*norm(x);
    v=v/norm(v);
    R(k:m,k:n)=R(k:m,k:n)-v*(2*(v'*R(k:m,k:n)));
    q=eye(max(m,n));
    q(k:max(m,n),k:max(m,n))=q(k:max(m,n),k:max(m,n))-2*v*v'/norm(v)^2;
    Q=Q*q;
  endfor

Méthode de Jacobi

La méthode de Jacobi permet de trouver une solution d’un système d’équations par iteration

function [x,it]= jacobi(A,b,x,maxit,dp)
  %auteur: Nathan De Smedt
  [m,n]=size(A)
  y=x;
  for it=1:maxit
     for i=1:n
       s=0;
       for j=1:n
         if j!=i
           s=s+A(i,j)*x(j);
         endif
       endfor
       y(i)=(b(i)-s)/A(i,i);
     endfor
     x=y;
     p=norm(A*x-b);
     if p<dp
       break
     endif
  endfor

Méthode de Gauss-Seidel

Cette méthode est similaire à celle de Jacobi, mais elle converge plus rapidement vers une solution.

function [x,it]= gausssidel(A,b,x,maxit,dp)
  %auteur: Nathan De Smedt
  [m,n]=size(A)
  for it=1:maxit
     for i=1:n
       s=0;
       for j=1:n
         if j!=i
           s=s+A(i,j)*x(j);
         endif
       endfor
       x(i)=(b(i)-s)/A(i,i);
     endfor
     p=norm(A*x-b);
     if p<dp
       break
     endif
  endfor

Dichotomie

Résolution d’équations non linéaire par dichotomie

function x=dichotomie(f,a,b,prec)
  %auteur: Nathan De Smedt
  x=a;
  if f(a)*f(b)>0
    "erreur"
  end   
  if f(a)==0
    x=a;
  end
  if f(b)==0
    x=b;
  end
  if f(a)*f(b)<0
    while abs(f(x))>=prec
      x=a/2+b/2;
      if f(x)*f(a)<0
        b=x;
      elseif f(x)*f(b)<0
        a=x;
      endif
    endwhile
  end
end 

Fausse position

Méthode de résolution d’équation non linéaire par fausse position

function x=faussepos(f,a,b, prec)
  %auteur: Nathan De Smedt
  if f(a)*f(b)<0
    x=a
    while f(x)*f(x)>prec*prec
      x=a-f(a)*(b-a)/(f(b)-f(a));
      if f(a)*f(x)<0
        b=x;
      elseif f(b)*f(x)<0  
        a=x;
      endif
      if f(x)*f(x)<prec*prec
        break
      endif
    endwhile
  endif

Newton-Raphson

Méthode de résolution d’équation non linéaire par Newton-Raphson.

function x=newtonraphson(f,fp,x, prec)
    %auteur: Nathan De Smedt
  while abs(f(x))>prec
    if fp(x)==0
      break
    endif
    x=x-f(x)/fp(x);
  endwhile

Newton-Raphson avec recherche linéaire

Méthode similaire a celle de Newton-Raphson mais en s’assurant de ne pas diverger.

function x=newtonraphsonr(f,fp,x, prec)
  %auteur: Nathan De Smedt
  while abs(f(x))>prec
    if fp(x)==0
      break
    endif
    a=1;
    xp=x;
    while abs(f(x))>=abs(f(xp))
      xp=x-a*f(x)/fp(x);
    endwhile
    x=xp;
  endwhile

Intégration par la formule des trapèze

Méthode d’intégration d’une fonction en calculant l’air de trapèze.

function A=trapeze(f,a,b,n)
  %auteur: Nathan De Smedt 
  h=(-a+b)/n;
  A=(f(a)+f(b))*h/2;
  for i=1:n-1
    A=A+h*f(a+i*h);
  endfor

Intégration par la formule de Simpson

Intégration sur un intervalle grâce a la formule de Simpson

function A=simpson(f,a,b,n)
  %auteur: Nathan De Smedt
  h=(-a+b)/n;
  A=(f(a)+f(b))*h/6;
  for i=1:n-1
    A=A+h/6*(2*f(a+i*h)+4*f(a+(i+1/2)*h));
  endfor

Intégration par la méthode de Romberg

Intégration sur un intervalle grâce a la méthode de Romberg

function x=romberg(f,a,b,k)
  %auteur: Nathan De Smedt
  I=zeros(k,k);
  for i=1:k
    I(1,i)=trapeze(f,a,b,2^(i-1));
  endfor
  for j=1:k-1
    for i=j+1:k
      I(j+1,i)=(4^j*I(j,i)-I(j,i-1))/(4^j-1);
      x=I(k,k);
      I
    endfor
  endfor

Euler progressif

Résolution d’équations différentiels aux conditions initials grâce a la méthode d’Euler progressive

function y=eulerp(f,y0,a,b,n)
  h=(b-a)/n;
  y=zeros(1,n);
  y(1)=y0;
  for i=1:n-1
    y(i+1)=y(i)+h*f(a+h*i,y(i));
  endfor

Laisser un commentaire