Programmation en langage C

Le programme ci-dessous implémente l'algorithme 4.2 en langage C.

#include <stdio.h>
#include <math.h>

#define NMAX 1001

/* 

programme elimin.c : elimination de Gauss (chap.4 du livre) 

Etant donne un entier N, deux reels a, c et un N-vecteur b,
appliquer l'agorithme d'elimination de Gauss,
au systeme lineaire Ax=b, ou A est la matrice tridiagonale 
symetrique de diagonale a et sous-diagonales -c.


parametres :

  N        : nombre d'inconnues du systeme lineaire
  a,c      : diagonale, sous-diagonale
  b        : N-vecteur, second membre du systeme lineaire en entree,
             solution du systeme lineaire Ax=b en sortie
  u        : N-vecteur utilise lors de l'elimination de Gauss

*/

void main (void){

  int N,i; double a,c,p; double b[NMAX],u[NMAX];
  
  N = 10;
  a = 2.;
  c = 1.;

/* initialisation du second membre de Ax=b */

  for (i=1;i<=N;i++) {
      b[i] = 1.;
  }

/* algorithme d'elimination de gauss */

  u[1] = -c/a;
  b[1] /= a;
  for (i=1;i<=N-1;i++) {
     p = a + c * u[i];        /* pivot */
     u[i+1] = -c / p;
     b[i+1] = (b[i+1] + c * b[i]) / p;
  }
  
/* resolution du systeme lineaire */

  for (i=N-1;i>=1;i--) {
    b[i] = b[i] - u[i] * b[i+1];
  }

/* impression des resultats */

  for (i=1;i<=N;i++) {
      printf(" i = %d x[i]-i*(N+1-i)/2 = %e \n",i,b[i]-i*(N+1.-i)/2.);
  }

}

Voilà le résultat du programme :

 i = 1 x[i]-i*(N+1-i)/2 = -8.881784e-16 
 i = 2 x[i]-i*(N+1-i)/2 = -1.776357e-15 
 i = 3 x[i]-i*(N+1-i)/2 = -1.776357e-15 
 i = 4 x[i]-i*(N+1-i)/2 = 0.000000e+00 
 i = 5 x[i]-i*(N+1-i)/2 = 0.000000e+00 
 i = 6 x[i]-i*(N+1-i)/2 = 0.000000e+00 
 i = 7 x[i]-i*(N+1-i)/2 = 0.000000e+00 
 i = 8 x[i]-i*(N+1-i)/2 = 0.000000e+00 
 i = 9 x[i]-i*(N+1-i)/2 = 0.000000e+00 
 i = 10 x[i]-i*(N+1-i)/2 = 0.000000e+00