Ubuntu-fr

wifiii

gsl pour résoudre un système linéaire

Bonjour,
je souhaite résoudre un système linéaire $Au^{n+1}= f^n + w$
ce système linéaire vient du schéma implicite
$$-a u^{n+1}_{j-1} + b u^{n+1}_j - c u^{n+1}_{j+1}=u^n_j$$
avec $u^0_j=0$ pour tout $j$ de 1 à $N$, $u^n_0=a$ et u^n_{N+1}=0$ pour tout $n$ de 0 à $M$.
où $A$ est une matrice tridiagonale de taille N*N, les éléments de sa diagonales sont $b$, de sa sous-diagonales: $-a$, et de sa sur-diagonale $-c$.
$(f^n)_j=(u^n)_j, j=1,...,N$, $w$ est un vecteur nul partout sauf pour sa première composante qui vaut a.
$n$ va de 0 à M, et $(u^0_j)_j = 0$ pour $j$ de 1 à N.

Je souhaite résoudre ce système en utilisant gsl. Mon problème est à quel moment je met la boucle d'itération sur n? Voici le programme que j'ai écris

#include <stdio.h>
#include <gsl/gsl_linalg.h>
#include<math.h>
#include <gsl/gsl_vector.h>


 int main (void)
{
    int N=99;
    int M=99;


//les data
double L=1000.;
double T=3860.;
double phi=0.3;
double D=0.1;
double v=0.05;


double delta=T/(M+1);
double h=L/(N+1);
double lambda1=(D/phi)*(delta/pow(h,2));
double lambda2=(v/phi)*(delta/h);
double a=lambda1+lambda2;
double b=(2*lambda1)+1+lambda2;
double c=lambda1;



double *a_data; //allouer un espace pour a_data
a_data =(double*) malloc(N*N*sizeof(double));
a_data[0]=b; 
a_data[1]= -c;
a_data[N*N-2]= -a;
a_data[N*N-1]=b;
for(int i=2;i<N;i++)
{
int j=(i-1)*N + i-1;
printf("j=%i\n",j);
a_data[j]=b;
a_data[j-1]=-a;
a_data[j+1]=-c;
printf("a_data[j]=%i\n",j);
}

double *f_data; //allouer un espace pour b_data
    f_data =(double*) malloc(N*sizeof(double));

    f_data[0]=a;
    for(int i=1; i<N;i++) f_data[i]=0;


  double *w_data; //allouer un espace pour b_data
    w_data =(double*) malloc(N*sizeof(double));

    w_data[0]=a;

   gsl_matrix_view m = gsl_matrix_view_array (a_data, N, N);

   gsl_vector_view f = gsl_vector_view_array (f_data, N);

   gsl_vector_view w = gsl_vector_view_array (w_data, N);
       

for (int i=1;i< M;99)
{

 
      gsl_vector *x = gsl_vector_alloc (N);

     int s;

     gsl_permutation * p = gsl_permutation_alloc (N);

     gsl_linalg_LU_decomp (&m.matrix, p, &s);


     gsl_linalg_LU_solve (&m.matrix, p, &f.vector, x);

     printf ("x = \n");
     gsl_vector_fprintf (stdout, x, "%g");

     gsl_permutation_free (p);
     
     
    int gsl_vector_memcpy (gsl_vector *f, const gsl_vector *x);  
     
    int gsl_vector_add(gsl_vector *f, const gsl_vector *w);
  
    gsl_vector_free (x);

   
    }
 return 0; 

}

Il ne donne que des 0, sûrement à cause de la boucle qui est mal placée. Comment l'aranger? S'il vous plaît. Je vous remercie par avance.

pingouinux

Re : gsl pour résoudre un système linéaire

Bonjour,
Je vois déjà ceci de curieux à la ligne 66 :

for (int i=1;i< M;99)

wifiii

Re : gsl pour résoudre un système linéaire

oui, c'est pour faire l'itération sur n. Comme le système est $A u^{n+1}=f^n$, ce qui nous interess, c'est la solution pour n=M-1 (trouver $u^{M-1}$, donc à chaque fois que l'on trouve un $u^n$, on l'utilise pour calculer $u^{n+1}$ et cela pour n de 1 à M-1 (c'est un schéma implicite). Mais je ne sais pas à quel moment mettre l'itération. Il y'a un souci pour faire l'itération avec gsl?
Je vous remercie par avance.

Dernière modification par wifiii (Le 19/09/2015, à 14:36)

pingouinux

Re : gsl pour résoudre un système linéaire

Ce que je veux dire en #2, c'est que l'indice de boucle i n'est pas incrémenté.

wifiii

Re : gsl pour résoudre un système linéaire

oui, pardon; donc le code devient comme ci-dessous, mais le résultat est toujours éronné.
En fait ce que je chercher exactement, c'est calculer la solution de $Au^1= f^0 $ où $f^0 = u^0 + w$ ainsi, je trouve $u^1$, puis je remplace le second membre $f^0$ par $u^1$, je lui rajoute $w$^et je recalculer la solution du système $A u^2 = f^1$ où $f^1 = u^1 + w$ et ainsi de suite, jusqu'à $M-1$.
Comment faut-il arranger le programme pour obtenir ca? S'il vous plaît.
Je vous remercie par avance.

wifiii

Re : gsl pour résoudre un système linéaire

En fait, j'ai avancé depuis tout à l'heure. Maintenant j'ai les 3 questions suivantes s'il vous plaît:
1- on sait que x[0]=1 et x[N+1]=0. Comment introduire ces deux valeurs dans les programme?

2- S'il te  plaît, crois-tu que la définition que j'ai implémenté pour la matrice A est correcte ? C'est une matrice tridiagonale où les éléments de la diagonale sont b, de la sous diagonale sont -a, et ceux de la sur diagonale sont -c. Peut-être que tu peux me proposer une meilleure façon d'implémenter cette matrice.

3- Comment fait-on pour mettre le dernier x obtenu dans un fichier.txt avec les x_i qui sont égalent à i*h ? (h est défini dans le programme).  J'ai essayé ceci:

#include <stdio.h>
#include <gsl/gsl_linalg.h>
#include<math.h>
#include <gsl/gsl_vector.h>


 int main (void)
{
    int N=99;
    int M=99;


//les data
double L=1000.;
double T=3860.;
double phi=0.3;
double D=0.1;
double v=0.05;


double delta=T/(M+1);
double h=L/(N+1);
double lambda1=(D/phi)*(delta/pow(h,2));
double lambda2=(v/phi)*(delta/h);
double a=lambda1+lambda2;
double b=(2*lambda1)+1+lambda2;
double c=lambda1;



double *a_data; //allouer un espace pour a_data
a_data =(double*) malloc(N*N*sizeof(double));
a_data[0]=b; 
a_data[1]= -c;
a_data[N*N-2]= -a;
a_data[N*N-1]=b;
for(int i=2;i<N;i++)
{
int j=(i-1)*N + i-1;
printf("j=%i\n",j);
a_data[j]=b;
a_data[j-1]=-a;
a_data[j+1]=-c;
printf("a_data[ j]=%i\n",j);
}

double *f_data; //allouer un espace pour b_data
    f_data =(double*) malloc(N*sizeof(double));

    f_data[0]=a;
    for(int i=1; i<N;i++) f_data[ i]=0;


  double *w_data; //allouer un espace pour b_data
    w_data =(double*) malloc(N*sizeof(double));

    w_data[0]=a;

   gsl_matrix_view m = gsl_matrix_view_array (a_data, N, N);

   gsl_vector_view f = gsl_vector_view_array (f_data, N);

   gsl_vector_view w = gsl_vector_view_array (w_data, N);
       


   for (int i=1;i< M;i++)
 {
   gsl_vector *x = gsl_vector_alloc (N);
   int s;
   gsl_permutation * p = gsl_permutation_alloc (N);
   gsl_linalg_LU_decomp (&m.matrix, p, &s);
   gsl_linalg_LU_solve (&m.matrix, p, &f.vector, x);
   printf ("x = \n");
   gsl_vector_fprintf (stdout, x, "%g");
   gsl_permutation_free (p);
   gsl_vector_memcpy (&f.vector,x);
   gsl_vector_add(&f.vector, &w.vector);
   

   {
FILE * f = fopen ("resultat.txt", "w");
gsl_vector_fprintf (f, x, "%.5g");
fclose (f);
}

  
   gsl_vector_free (x);
  }



 return 0; 

}

il n'affiche que x, mais ce que je cherche, c'est qu'il affiche x_i=i*h et x[ i] pour tout i de 0 à N+1, pour que je puisse faire le graphe après. Qu'est ce qu'il faut utiliser pour ça? S'il vous plaît.

Je vous remercie par avance.

pingouinux

Re : gsl pour résoudre un système linéaire

Ne connaissant pas gsl, je ne sais pas répondre aux questions 1 et 3.
Pour la question 2, je pense que ceci devrait être correct :

for(i=0;i<N*N;i++)    a_data[i]=0;
for(i=0;i<N*N;i+=N+1) a_data[i]=b;
for(i=1;i<N*N;i+=N+1) a_data[i]=-c;
for(i=N;i<N*N;i+=N+1) a_data[i]=-a;

Ajouté : Initialisation à zéro

Dernière modification par pingouinux (Le 19/09/2015, à 20:30)

wifiii

Re : gsl pour résoudre un système linéaire

Non, cette écriture est incompatible avec gsl. Voici ce qu'il me renvoie:
In function ‘main’:
matrice.c:34:5: error: ‘i’ undeclared (first use in this function)
for(i=0;i<N*N;i+=N+1) a_data[ i]=b;
     ^
matrice.c:34:5: note: each undeclared identifier is reported only once for each function it appears in
matrice.c:76:4: error: incompatible type for argument 1 of ‘gsl_vector_get’
    gsl_vector_get(v, 0)=1;
    ^
In file included from /usr/include/gsl/gsl_vector_complex_double.h:28:0,
                 from /usr/include/gsl/gsl_vector.h:5,
                 from /usr/include/gsl/gsl_linalg.h:25,
                 from matrice.c:2:
/usr/include/gsl/gsl_vector_double.h:166:20: note: expected ‘const struct gsl_vector *’ but argument is of type ‘double’
INLINE_DECL double gsl_vector_get (const gsl_vector * v, const size_t i);

pingouinux

Re : gsl pour résoudre un système linéaire

matrice.c:34:5: error: ‘i’ undeclared (first use in this function)
for(i=0;i<N*N;i+=N+1) a_data[ i]=b;
     ^
matrice.c:34:5: note: each undeclared identifier is reported only once for each function it appears in

J'indiquais simplement une méthode relativement simple pour initialiser la matrice. Il est évident que si une variable n'est pas déclarée, il faut le faire avant (ou dans le for, comme tu l'as fait).

wifiii

Re : gsl pour résoudre un système linéaire

Bonjour,
en fait ca compile bien. Mais comment faire pour afficher la matrice a_data pour que je vois si c'est bon? J'ai essayé
printf("a_data[ i]=%i\n");

mais ca m'affiche un nombre bizare a_data[ i]=11747472
Je vous remercie par avance pour votre aide.

wifiii

Re : gsl pour résoudre un système linéaire

j'ai essayé d'implémenter une matrice tridiagonale 4*4, et voici le code

#include <stdio.h>
#include <gsl/gsl_linalg.h>
#include<math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>

int
main (void)
{
int i;
gsl_vector * diag = gsl_vector_alloc (4);
for (i = 0; i < 4; i++)
{
gsl_vector_set (diag, i, 2);
}
{
FILE * f = fopen ("test.dat", "w");
gsl_vector_fprintf (f, diag, "%.5g");
fclose (f);
}
gsl_vector_free (diag);

//////////////////////

gsl_vector * e = gsl_vector_alloc (4);
for (i = 0; i < 3; i++)
{
gsl_vector_set (e, i, -1);
gsl_vector_set(e,3,0);
}
{
FILE * f = fopen ("test.dat", "w");
gsl_vector_fprintf (f, e, "%.5g");
fclose (f);
}
gsl_vector_free (e);

//////////////////


gsl_vector * w = gsl_vector_alloc (4);
for (i = 1; i <=3; i++)
{
gsl_vector_set (w, i, -1);
gsl_vector_set(w,0,0);
}
{
FILE * f = fopen ("test.dat", "w");
gsl_vector_fprintf (f, w, "%.5g");
fclose (f);
}
gsl_vector_free (w);


//////////////////



gsl_vector * b = gsl_vector_alloc (4);
for (i = 1; i <=3; i++)
{
gsl_vector_set (b, i, 0);
gsl_vector_set(b,0,1);
}
{
FILE * f = fopen ("test.dat", "w");
gsl_vector_fprintf (f, b, "%.5g");
fclose (f);
}
gsl_vector_free (b);


/////////////////////////

//gsl_vector *x = gsl_vector_alloc (4);



/////////////////////////


int j;
gsl_matrix * m = gsl_matrix_alloc (4, 4);

gsl_matrix_set (m, 0, 0,2);
gsl_matrix_set (m, 0, 1, -1);
gsl_matrix_set (m, 3, 2, -1);
gsl_matrix_set (m, 3, 3, 2);

for (i = 1; i < 2; i++)
for (j = 0; j < 3; j++)
{
  if (j==i-1) 
  gsl_matrix_set (m, i, j,-1);
else
  if (j==i)
  gsl_matrix_set (m, i, j,2);
else
 if (j==i+1)
  gsl_matrix_set (m, i, j, -1);
else 
gsl_matrix_set (m, i, j, 0);
}




for (i = 0; i < 3; i++) /* OUT OF RANGE ERROR */
for (j = 0; j < 3; j++)
printf ("m(%d,%d) = %g\n", i, j,
gsl_matrix_get (m, i, j));
gsl_matrix_free (m);


return 0;
}

et en fait il ne m'affiche que ca
m(0,0) = 2
m(0,1) = -1
m(0,2) = 0
m(1,0) = -1
m(1,1) = 2
m(1,2) = -1
m(2,0) = 0
m(2,1) = 0
m(2,2) = 0

il y'a des éléments de la matrice qui manquents. Je ne vois vraiment pas comment arranger. Je vous remercie par avance.

Dernière modification par wifiii (Le 20/09/2015, à 16:38)

wifiii

Re : gsl pour résoudre un système linéaire

pingouinux peux-tu m'expliquer comment tu as raisonné pour construire la matrice que tu me propose?
Merci par avance.

pingouinux

Re : gsl pour résoudre un système linéaire

Voici un petit programme (fait à partir d'un extrait du tien), qui montre comment imprimer la matrice, et qui compare nos deux méthodes.
Tu remarqueras que j'y utilise calloc, qui initialise la matrice à 0, au lieu de malloc qui ne le fait pas (même si le hasard fait que ça marche dans certains cas).

$ cat prog.cc
#include <stdio.h>
#include <stdlib.h>

 int main (void)
{
    int N=5;

  //les data
    double a=10;
    double b=20;
    double c=30;

    double *a_data; //allouer un espace pour a_data
 // Ta méthode
    printf("a_data (ta méthode)\n");
    a_data =(double*) malloc(N*N*sizeof(double));

    a_data[0]=b; 
    a_data[1]= -c;
    a_data[N*N-2]= -a;
    a_data[N*N-1]=b;
    for(int i=2;i<N;i++)
    {
       int j=(i-1)*N + i-1;
       printf("j=%i\n",j);
       a_data[j]=b;
       a_data[j-1]=-a;
       a_data[j+1]=-c;
    }
 // Impression de la matrice
    int k;
    for(k=0;k<N*N;k++){
       printf(" %10g",a_data[k]);
       if(!((k+1)%N)) printf("\n");
    }
    printf("\n");
    free(a_data);

 // Ma méthode
    printf("a_data (ma méthode)\n");
    a_data =(double*) calloc(N*N,sizeof(double));
    for(k=0;k<N*N;k+=N+1) a_data[k]=b;     // Diagonale principale
    for(k=1;k<N*N;k+=N+1) a_data[k]=-c;    // Diagonale du dessus
    for(k=N;k<N*N;k+=N+1) a_data[k]=-a;    // Diagonale du dessous
 // Impression de la matrice
    for(k=0;k<N*N;k++){
       printf(" %10g",a_data[k]);
       if(!((k+1)%N)) printf("\n");
    }
    printf("\n");

    return 0; 
}

à utiliser ainsi

$ make prog && ./prog
g++     prog.cc   -o prog
a_data (ta méthode)
j=6
j=12
j=18
         20        -30          0          0          0
        -10         20        -30          0          0
          0        -10         20        -30          0
          0          0        -10         20        -30
          0          0          0        -10         20

a_data (ma méthode)
         20        -30          0          0          0
        -10         20        -30          0          0
          0        -10         20        -30          0
          0          0        -10         20        -30
          0          0          0        -10         20

La matrice a_data est stockée sous la forme d'un tableau mono-dimensionnel, les lignes se succédant (ligne 0 à ligne N-1).
L'indice k (dans le grand tableau) de l'élément situé à la ligne i et à la colonne j est : k = i*N + j.
Quand k augmente de N, on passe d'un élément à celui situé juste en-dessous. Pour rester sur une même diagonale, il faut incrémenter k de N+1.

Édité : Petites corrections

Dernière modification par pingouinux (Le 21/09/2015, à 10:31)

pingouinux

Re : gsl pour résoudre un système linéaire

J'ai repris le programme en #13, et y ai ajouté deux autres façons d'initialiser la matrice. L'impression se fait maintenant dans la fonction impression_matrice.

#include <stdio.h>
#include <stdlib.h>

void impression_matrice(double* tab, int N) {
    for(int k=0;k<N*N;k++){
       printf(" %10g",tab[k]);
       if(!((k+1)%N)) printf("\n");
    }
}
int main (void)
{
    int N=5;
    int i, j, k;

  //les data
    double a=10;
    double b=20;
    double c=30;

    double *a_data; //allouer un espace pour a_data
 // Ta méthode
    printf("a_data (ta méthode)\n");
    a_data =(double*) malloc(N*N*sizeof(double));
 // Initialisation de la matrice
    a_data[0]=b; 
    a_data[1]= -c;
    a_data[N*N-2]= -a;
    a_data[N*N-1]=b;
    for(int i=2;i<N;i++)
    {
       int j=(i-1)*N + i-1;
       printf("j=%i\n",j);
       a_data[j]=b;
       a_data[j-1]=-a;
       a_data[j+1]=-c;
    }
 // Impression de la matrice
    impression_matrice(a_data,N);

 // Ma méthode
    free(a_data);
    printf("a_data (ma méthode)\n");
    a_data =(double*) calloc(N*N,sizeof(double));
 // Initialisation de la matrice
    for(k=0;k<N*N;k+=N+1) a_data[k]=b;     // Diagonale principale
    for(k=1;k<N*N;k+=N+1) a_data[k]=-c;    // Diagonale du dessus
    for(k=N;k<N*N;k+=N+1) a_data[k]=-a;    // Diagonale du dessous
 // Impression de la matrice
    impression_matrice(a_data,N);

 // Autre méthode 1
    free(a_data);
    printf("a_data (autre méthode 1)\n");
    a_data =(double*) calloc(N*N,sizeof(double));
 // Initialisation de la matrice
    for(i=0;i<N;i++) for(j=0;j<N;j++) {
       k=i*N+j;
       if     (i==j)   a_data[k]=b;     // Diagonale principale
       else if(i+1==j) a_data[k]=-c;    // Diagonale du dessus
       else if(i-1==j) a_data[k]=-a;    // Diagonale du dessous
    }
 // Impression de la matrice
    impression_matrice(a_data,N);

 // Autre méthode 2
    free(a_data);
    printf("a_data (autre méthode 2)\n");
    a_data =(double*) calloc(N*N,sizeof(double));
 // Initialisation de la matrice
    for(k=0;k<N*N;k++) {
       if     ((k%(N+1))==0) a_data[k]=b;     // Diagonale principale
       else if((k%(N+1))==1) a_data[k]=-c;    // Diagonale du dessus
       else if((k%(N+1))==N) a_data[k]=-a;    // Diagonale du dessous
    }
 // Impression de la matrice
    impression_matrice(a_data,N);

    return 0; 
}