Discussion :

[hinaneee]

Bonjour,
je cherche un exemple qui montre comment utiliser la bibliothèque GSL pour résoudre un système linéaire.
Je vous remercie par avance.

[prgasp77]

Voila : http://bfy.tw/1q0A !
Documentation §2.1 ... ho ! un exemple !

[hinaneee]

vous pensez vraiment que je n'ai pas essayé ca avant de venir vous poser la question?
je ne trouve pas d'exemple qui montre comment résoudre un système linéaire (résolution au complet). Merci.

[hinaneee]

Bonjour,
j'ai le schéma suivant

et et sont données.
Je l'ai écris à la main sous la forme d'un système matriciel. Comment l'écrire en code c++ pour ensuite le résoudre en utilisant GSL? (en fait on commence par son écriture) s'il vous plaît. Je vous remercie par avance.

[ternel]

Conformément à la charte que nous avons tous acceptés, je ne peux que te demander de nous montrer ce que tu as déjà fait.
Même si ca ne fonctionne pas.

C'est la seule façon de nous permettre de t'aider.

[prgasp77]

Bonjour,
j'ai le schéma suivant

et et sont données.
Je l'ai écris à la main sous la forme d'un système matriciel. Comment l'écrire en code c++ pour ensuite le résoudre en utilisant GSL? (en fait on commence par son écriture) s'il vous plaît. Je vous remercie par avance.

Admettons. Est-ce qu'on peut voir ton système sous forme matriciel ? Je ne suis pas sûr de comprendre à quoi correspond ton supscript dans ta notation ; c'est de l'exponentielle ? developpez.com n'a pas de moteur TeX, utilise celui de google chart https://chart.googleapis.com/chart?cht=tx&chl=\left[\begin{array}a%26b\\c%26d\end{array}\right] ->

Pour un exemple de résolution de système par décomposition LU : http://www.gnu.org/software/gsl/manu...gebra-Examples

[hinaneee]

Bon, pour le moment, je commence par écrire mon système matriciel en C++. J'ai réussi à écrire la matrice A, par contre pour l'inconnue X et le second membre F qui sont des vecteur, je ne sais pas comment faire rentrer l'iconue X et le second membre F. Voici mon fichier .cpp

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
 
#include<iostream>
#include<boost/numeric/ublas/matrix.hpp>
#include<boost/numeric/ublas/io.hpp>
using namespace std;
#include<vector>
 
int main()
{
using namespace boost::numeric::ublas;
int N=6;
matrix<double> m(N,N);
for (unsigned i=0;i<6;++i)
{
m(i,i)=1;
m(i+1,i)=2;
m(i,i+1)=3;
std::cout<<m<<std::endl;
}
};

je vous remercie par avance.

[prgasp77]

Je suis désolé, mais tu as tout faut. Il faut utiliser les types fournis par la bibliothèque. Regarde donc l'exemple.

De plus, vire donc ces using namespace qui sentent vraiment mauvais ... ton type matrix<double, double> il provient de std:: ou de boosl::numeric::ublas:: ?

[dragonjoker59]

certainement pas de std...

[hinaneee]

Voici ce que j'ai écrit pour résoudre le système $A u^{n+1}=b^n$

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
 
#include <stdio.h>
#include <gsl/gsl_linalg.h>
 
int
main (void)
{
 int N=4;
int M=4;
 
  double a_data[N,N]; 
for (int i=1;i<N;i++)
{
for (int j=1;j<N;j++)
{
if(i==j) a_data[i,j]=1;
else
if(j==i-1) a_data[i,j]=60;
else
if (j==i+1) a_data[i,j]=7;
else
a_data[i,j]=0;
}
}
 
double b_data[N];
b_data[1]=1;
for(i=2: i<N;i++)
b_data[i]=0;
 
 
for(int n=0; n<M; n++)
{
  gsl_matrix_view m 
    = gsl_matrix_view_array (a_data, N, N);
 
  gsl_vector_view b
    = gsl_vector_view_array (b_data, N);
 
  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, &b.vector, x);
 
  printf ("x = \n");
  gsl_vector_fprintf (stdout, x, "%g");
 
  gsl_permutation_free (p);
  gsl_vector_free (x);
  b_data=x;
  return 0;
}
}

il me renvoie ces erreurs:
matrice.c:10:10: error: array size missing in ‘a_data’
double a_data[];
^
matrice.c:11:1: error: ‘for’ loop initial declarations are only allowed in C99 or C11 mode
for (int i=1;i<N;i++)
^
matrice.c:11:1: note: use option -std=c99, -std=gnu99, -std=c11 or -std=gnu11 to compile your code
matrice.c:13:1: error: ‘for’ loop initial declarations are only allowed in C99 or C11 mode
for (int j=1;j<N;j++)
^
matrice.c:15:18: warning: left-hand operand of comma expression has no effect [-Wunused-value]
if(i==j) a_data[i,j]=1;
^
matrice.c:17:20: warning: left-hand operand of comma expression has no effect [-Wunused-value]
if(j==i-1) a_data[i,j]=60;
^
matrice.c:19:21: warning: left-hand operand of comma expression has no effect [-Wunused-value]
if (j==i+1) a_data[i,j]=7;
^
matrice.c:21:9: warning: left-hand operand of comma expression has no effect [-Wunused-value]
a_data[i,j]=0;
^
matrice.c:25:8: error: array size missing in ‘b_data’
double b_data[];
^
matrice.c:27:8: error: expected ‘;’ before ‘:’ token
for(i=2: i<N;i++)
^
matrice.c:27:17: error: expected ‘;’ before ‘)’ token
for(i=2: i<N;i++)
^
matrice.c:31:1: error: ‘for’ loop initial declarations are only allowed in C99 or C11 mode
for(int n=0; n<M; n++)
^
matrice.c:54:9: error: assignment to expression with array type
b_data=x;
^
matrice.c:57:1: warning: control reaches end of non-void function [-Wreturn-type]
}
^

Ce que je remarque, c'est qu'il n'accepte pas la boucle for, pourtant il faut calculer à chaque fois u^{n+1} en utilisant u^n, et en partant du vecteur b_data qui a N composante qui sont toute nulles à part la première.
Comment l'arranger? S'il vous plaît. Je vous remercie par avance.

[prgasp77]

Aïe ! Il va falloir revoir les bases. double[M,N] n'a aucune chance de fonctionner. Regarde donc l'exemple : ils se basent sur un tableau unidimensionnel. Oublie les boucles for imbriquées pour l'instant et écris en dur la valeur du tableau :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
double a_data[] = { 0.18, 0.60, 0.57, 0.96,
                      0.41, 0.24, 0.99, 0.58,
                      0.14, 0.30, 0.97, 0.66,
                      0.51, 0.13, 0.19, 0.85 };

certainement pas de std...

(question rhétorique )

[hinaneee]

En fait, j'avais écris double m[], mais il l'a refusé, j'ai donc écris cette bêtise.
Si j'écris les valeurs de A de manière explicite, il n'y a aucun problème, ca marche, je l'ai essayé, mais ce n'est pas ce qui m'interesse. La matrice A est de taille N où N=99, elle est tri-diaginale, on ne peut pas la rentrer à la main, donc comment on fait pour la rentrer si on ne peut pas utiliser les boucles for?
Je vous remercie par avance.

[prgasp77]

Je comprends, mais dans la découverte d'une bibliothèque, il faut avancer à petits pas. Mieux vaut un petit système entré en dur pour commencer afin de se simplifier la vie et vérifier que les résultats obtenus sont bien ceux que l'on attendait. Donc dans un premier temps, je te propose de résoudre un système de dimension 3 avec la lib. Montre nous le code et le résultat de son exécution. Ensuite on verra pour résoudre n'importe quel système.

[hinaneee]

Voici le code pour une matrice à 3 composantes:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
 
#include <stdio.h>
#include <gsl/gsl_linalg.h>
 
int 
 
main (void)
{
 int N=3;
double a_data[] = { 1,3, 0,
                      5, 1, 3,
                      0, 5, 1};
 
double b_data[] = { 1, 0, 0};
 
  gsl_matrix_view m 
    = gsl_matrix_view_array (a_data, N, N);
 
  gsl_vector_view b
    = gsl_vector_view_array (b_data, N);
 
  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, &b.vector, x);
 
  printf ("x = \n");
  gsl_vector_fprintf (stdout, x, "%g");
 
  gsl_permutation_free (p);
  gsl_vector_free (x);
  return 0;
 
}

et voici la solution obtenue:
x =
0.482759
0.172414
-0.862069

Maintenant, comment généraliser à un système n*n? J'ai lu qu'il faut une fonction spécifique pour introduire une matrice tri-diagonale, mais je ne comprend pas comment l'utiliser. Je vous remercie pour votre aide.

[dragonjoker59]

Ta matrice tri-diagonale, tu la récupères à partir de quoi?

[hinaneee]

Pardon, je ne comprend pas ta question. on construit la matrice tri-diagonale à parir de 3 vecteurs qui formes les 3 diagonales, et le reste c'est des zéros, j'avais essayé de la définir comme ca dans un de mes messages, mais il n'a pas accepté le for comme je l'ai dis. Vous pouvez m'aider?

[prgasp77]

Maintenant, comment généraliser à un système n*n? J'ai lu qu'il faut une fonction spécifique pour introduire une matrice tri-diagonale, mais je ne comprend pas comment l'utiliser. Je vous remercie pour votre aide.

Bon début. Tu as vérifié que le résultat obtenu était bon ?
Maintenant, pour traiter le cas n*n, je te propose de transformer ceci

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
double a_data[] = { 1,3, 0,
                      5, 1, 3,
                      0, 5, 1};

en code dymanique. Pour cela, il te faut une classe sachant accéder à l'élément [ligne, colonne] d'une matrice et sachant retourner les éléments de cette matrice sous la forme d'un tableau unidimensionnel.

Je te propose un début d'implémentation (utilisant C++11, si ce n'est pas une option pour toi, remplace std::array<T,N> par std::vector<T> et _data({0}) par _data(N*N, 0)), non relu, non essayé :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
namespace dvp
{
template<class T, size_t N>
class matrix
{
    std::array<T,N*N> _data;
public:
    matrix() : _data({0}) {}
    // copy & move constructors & operator= left to lecter
    T& operator()(size_t line, size_t column) { return _data[N * column + line]; }
    const T& operator()(size_t line, size_t column) const { return _data[N * column + line]; }
    T* data() { return _data.data(); }
    const T* data() const { return _data.data(); }
};
}

Et y'a qu'à qu'on !

[hinaneee]

Voilà, j'ai pris une matrice 4*4 pour bien montrer la boucle for, et voici le code:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
 
 
#include <stdio.h>
#include <gsl/gsl_linalg.h>
 
 
int
 
 
namespace dvp
{
template<class T, size_t N>
class matrix
{
    std::array<T,N*N> _data;
public:
    matrix() : _data({0}) {}
    // copy & move constructors & operator= left to lecter
    T& operator()(size_t line, size_t column) { return _data[N * column + line]; }
    const T& operator()(size_t line, size_t column) const { return _data[N * column + line]; }
    T* data() { return _data.data(); }
    const T* data() const { return _data.data(); }
};
}
 
main (void)
{
 matrix a_data(4,4);
  double a_data(1,1)=1;
  double a_data(1,2)=3;
  double a_data(4,3)=5;
  double a_data(4,4)=1;
for (int i=2;i<4;i++)
{
  a(i,i)=1;
  a(i,i+1)=3;
  a(i+1,i)=5;
}
 
  double b_data[] = { 1.0, 2.0, 3.0, 4.0 };
 
  gsl_matrix_view m 
    = gsl_matrix_view_array (a_data, 4, 4);
 
  gsl_vector_view b
    = gsl_vector_view_array (b_data, 4);
 
  gsl_vector *x = gsl_vector_alloc (4);
 
  int s;
 
  gsl_permutation * p = gsl_permutation_alloc (4);
 
  gsl_linalg_LU_decomp (&m.matrix, p, &s);
 
  gsl_linalg_LU_solve (&m.matrix, p, &b.vector, x);
 
  printf ("x = \n");
  gsl_vector_fprintf (stdout, x, "%g");
 
  gsl_permutation_free (p);
  gsl_vector_free (x);
  return 0;
}

mais il me renvoie les erreurs suivantes:

exemple2.c:26:16: error: expected declaration specifiers or ‘...’ before numeric constant
matrix a_data(4,4);
^
exemple2.c:26:18: error: expected declaration specifiers or ‘...’ before numeric constant
matrix a_data(4,4);
^
exemple2.c:27:17: error: expected declaration specifiers or ‘...’ before numeric constant
double a_data(1,1)=1;
^
exemple2.c:27:19: error: expected declaration specifiers or ‘...’ before numeric constant
double a_data(1,1)=1;
^
exemple2.c:28:17: error: expected declaration specifiers or ‘...’ before numeric constant
double a_data(1,2)=3;
^
exemple2.c:28:19: error: expected declaration specifiers or ‘...’ before numeric constant
double a_data(1,2)=3;
^
exemple2.c:29:17: error: expected declaration specifiers or ‘...’ before numeric constant
double a_data(4,3)=5;
^
exemple2.c:29:19: error: expected declaration specifiers or ‘...’ before numeric constant
double a_data(4,3)=5;
^
exemple2.c:30:17: error: expected declaration specifiers or ‘...’ before numeric constant
double a_data(4,4)=1;
^
exemple2.c:30:19: error: expected declaration specifiers or ‘...’ before numeric constant
double a_data(4,4)=1;
^
exemple2.c:33:3: warning: implicit declaration of function ‘a’ [-Wimplicit-function-declaration]
a(i,i)=1;
^
exemple2.c:33:9: error: lvalue required as left operand of assignment
a(i,i)=1;
^
exemple2.c:34:11: error: lvalue required as left operand of assignment
a(i,i+1)=3;
^
exemple2.c:35:11: error: lvalue required as left operand of assignment
a(i+1,i)=5;
^
exemple2.c:41:30: error: ‘a_data’ undeclared (first use in this function)
= gsl_matrix_view_array (a_data, 4, 4);
^







. Je vous remercie par avance.

[dragonjoker59]

As-tu déjà utilisé les templates? Parce que j'ai bien l'impression que tu n'as pas tout à fait compris comment ça marche.

[hinaneee]

En fait non, pas dans des cas comme celu ci. Je sais qu'on utilise des classe templates pour que la classe soit défini pour n'importe quel type d'attribut, mais là comment on l’appelle pour construire l'objet a_data? S'il vous plaît. Je vous remercie par avance.