systèmes linéaires avec gsl

Version imprimable

17/09/2015, 11h28
hinaneee

systèmes linéaires avec gsl

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.
17/09/2015, 11h33
prgasp77

Voila : http://bfy.tw/1q0A !
Documentation §2.1 ... ho ! un exemple !
17/09/2015, 11h38
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.
17/09/2015, 12h29
hinaneee

Bonjour,
j'ai le schéma suivant
$Formule mathématique$
$Formule mathématique$ et $Formule mathématique$ et $Formule mathématique$ 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.
17/09/2015, 13h10
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.
17/09/2015, 14h24
prgasp77

Citation:

Envoyé par hinaneee

Bonjour,
j'ai le schéma suivant
$Formule mathématique$
$Formule mathématique$ et $Formule mathématique$ et $Formule mathématique$ 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] -> https://chart.googleapis.com/chart?c...{array}\right]

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

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:

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.

17/09/2015, 16h23
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:: ?
17/09/2015, 16h37
dragonjoker59

certainement pas de std...

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

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.

17/09/2015, 17h38
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:

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 };
Citation:

Envoyé par dragonjoker59

certainement pas de std...

(question rhétorique ;))
17/09/2015, 17h49
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.
17/09/2015, 19h02
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.

Voici le code pour une matrice à 3 composantes:

Code:

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.

17/09/2015, 21h44
dragonjoker59

Ta matrice tri-diagonale, tu la récupères à partir de quoi?
17/09/2015, 21h49
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?

Citation:

Envoyé par hinaneee

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:

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:

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 !

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

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.

18/09/2015, 12h59
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.
18/09/2015, 13h05
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.
18/09/2015, 13h21
dragonjoker59

Si tu regardes la déclaration de la classe matrix, tu vois qu'elle a 2 paramètres template, T et N.
Si tu regardes la déclaration de matrix::_data, celle-ci te permet de comprendre (au cas où tu n'aurais pas déjà compris) que N sert de dimension pour une matrice carrée, et T sert de type de données contenues dans la matrice.
Lorsque tu veux créer une instance de la classe matrix, il faut donc que tu spécifies ces 2 paramètres, ce que tu ne fais pas dans ta ligne matrix a_data(4,4);.
De plus si tu lis bien la déclaration de cette matrix, tu remarqueras qu'elle n'a aucun constructeur prenant 2 paramètres, donc définitivement, ta ligne de déclaration de a_data a des problèmes.

EDIT:

Je ne te mets pas la solution, je ne fais que t'aiguiller pour que tu la trouves, la solution.

D'accord, donc j'ai modifié et j'ai bien mis deux paramètres: le type (double ici), et la taille qui est 4 ici. Voici le fichier.c

Code:

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}) {}//constructeur par défaut
    matrix(size_t line=0, size_t column=0)//constructeur avec arguments
    // 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 <double> a_data(class double, size_t 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;
}

par contre il y'a toujours des erreurs. Je vous remercie pour votre aide.

18/09/2015, 14h14
dragonjoker59

Peux-tu m'expliquer textuellement ce qu'est censée faire cette ligne?
matrix <double> a_data(class double, size_t 4);
18/09/2015, 14h32
hinaneee

En fait, ca deverai être comme ca:

Code:

matrix a_data(class double, size_t 4);

et ca veut dire que a_data est un objet de la classe matrice, t.q ses éléments sont de type double, et sa taille est 4.
18/09/2015, 14h45
dragonjoker59

Citation:

Envoyé par hinaneee

En fait, ca deverai être comme ca:

Code:

matrix a_data(class double, size_t 4);

et ca veut dire que a_data est un objet de la classe matrice, t.q ses éléments sont de type double, et sa taille est 4.

Qu'est-ce qui te fait dire que ça devrait être comme ça?
18/09/2015, 15h08
hinaneee

C'est parce qu'on a mis template<class T, size_t N> , donc les éléments de la classe ont deux arguments: le type et la taille, je l'ai compris comme ca.
18/09/2015, 15h18
ternel

comment écrirais-tu l'appel de la fonction int somme(int a, int b)?
18/09/2015, 15h23
hinaneee

a_data . somme(2,3) par exemple.
18/09/2015, 15h51
ternel

Ca, ca aurait été int matrix::somme(int a, int b), mais l'idée y est.

Tu as bien consciencieusement supprimé les indications de type.
Concrètement, ce que tu as écrit est une expression d'appel de fonction.

Avec les templates de classes, c'est le type qu'il faut paramétrer.
c'est ce que tu fais entre les chevrons qui suivent le nom de la template, comme dans vector<double>.

Dans les parenthèses, on écrit des arguments de l'expression.
Tu as écris ceci:

Code:

matrix <double> a_data(class double, size_t 4);

L'idée derrière est précisément la même que dans int i=4; ou précisément std::string s("bonjour");.

Il s'agit de déclarer une variable nommée a_data, du type obtenu par la template matrix pour les paramètres double et 4, et avec le constructeur qui va bien.

Avec ces éclaircissement, comment l'écrirais-tu?
Et surtout pourquoi?

Donc ca doit être

Code:

1
2
 
 matrix  <double, 4> a_data

C'est bien ca?
et voici le code

Code:

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
 
#include <stdio.h>
#include <gsl/gsl_linalg.h>
#include<math.h>
#include <gsl/gsl_vector.h>
 
 
namespace dvp
{
template<class T, size_t N>
class matrix
{
    std::array<T,N*N> _data;
public:
    matrix() : _data({0}) {}
    matrix(size_t line=0, size_t column=0)//constructeur avec arguments
    // 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(); }
};
}
 
int 
main(void)
{
 
matrix<double, 4> a_data;
 
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;
}
 
}

Qu'en pensez-vous?

Désolée pour le multi-poste, mais il y'a du nouveau. En fait, j'ai fini par écrire le code sans les templates,
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
Code:

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 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 #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.

Citation:
Envoyé par hinaneee
C'est bien ca?
Code:

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 #include <stdio.h> #include <gsl/gsl_linalg.h> #include<math.h> #include <gsl/gsl_vector.h> namespace dvp { template<class T, size_t N> class matrix { std::array<T,N*N> _data; public: matrix() : _data({0}) {} matrix(size_t line=0, size_t column=0)//constructeur avec arguments // 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(); } }; } int main(void) { matrix<double, 4> a_data; 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; } }
Qu'en pensez-vous?
Ca m'a l'air bien ça (il va peut-être te manquer un dvp:: mais sinon c'est tout bon). Il ne reste plus qu'à utiliser la matrice crée de la même manière que dans ton message #14. Petit indice, il faut utiliser a_data.data() qui est homogène au type const double[] de la même manière que tu utilisais a_data dans gsl_matrix_view m = gsl_matrix_view_array (a_data, N, N);.
Citation:
Envoyé par hinaneee
Désolée pour le multi-poste, mais il y'a du nouveau. En fait, j'ai fini par écrire le code sans les templates,
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
Code:

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 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 #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.
Trop compliqué pour moi. Il faut avancer étape par étape.

25/09/2015, 15h08
hinaneee

Bonjour prgasp77. S'il te plaît, dans ta définition de la classe template matrix, où est le constructeur avec argument et le constucteur de copie? Et par quelle logique tu les as implémenté?
Merci par avance.
25/09/2015, 16h03
prgasp77

Citation:

Envoyé par hinaneee

Bonjour prgasp77. S'il te plaît, dans ta définition de la classe template matrix, où est le constructeur avec argument et le constucteur de copie? Et par quelle logique tu les as implémenté?
Merci par avance.

Bonjour, ils sont ligne 9 (laissés au lecteur). Ils sont assez triviaux à écrire, je te laisse essayer.
25/09/2015, 16h43
hinaneee
Oui, je les avais repéré; je parlais de la définition
Code:

1 2 3 4 matrix(size_t line=0, size_t column=0)//constructeur avec arguments T& operator()(size_t line, size_t column) { return _data[N * column + line]; }// constructeur de copie
Mon problème est qu'à chaque fois je cherche d'anciens programme pour faire copier-coller de ces deux définitions, mais je vais les retenir. Peux-tu m'expliquer qu'est ce qu'on doit mettre dans la définition d'un constructeur avec arguments? (moi j'essairai de traduire en C++).
Pour le constructeur de copie, vous l'avez déjà définie.
Dans une classe de vecteur template, on définie le constructeur de copie par:
Code:

1 2 3 4 5 6 7 8 template <typename Real> vector<Real>::vector(const vector& ref): N(ref.N) { v = new Real[N]; for(int i=0; i<N;i++) v[i]=ref.v[i]; }
c'est à dire qu'il prend un vecteur de reférence et il le copie, mais ca sert à quoi exactement le constructeur de copie? Et pour une matrice, est-ce qu'on peut faire
Code:

1 2 3 4 5 6 7 8 matrix (cont matrix& ref): N*N(ref.N*N) { m=new [N*N] for (int i=0; i<N; i++) for(int j=0; j<N;j++) a_data[i,j]= ref. a_data[i,j]; }
Je vous remercie par avance de me corriger et m'indiquer les bêtises et comment faire afin de les retenir.

Repartons sur ce code :

Code:

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(); }
};
}

Il n'y a nul besoin de constructeur avec argument.
Le constructeur par copie n'est pas défini. L'objectif du constructeur par copie est de pouvoir écrire

Code:

1
2
3
dvp::matrix<int, 4> source;
source(0,0) = 42;
dvp::matrix<int, 4> copy(source);

et d'avoir les objets source et copy dans le même état (c'est à dire copy(0,0) == 42).

1- Pourquoi on n'a pas besoin de constructeur avec arguments? Et quand est-ce qu'on en a besoin?
2- Dans les lignes du code
Code:

1 2 3 4 5 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(); }
peux-tu s'il te plaît, m'expliquer avec des mots, ce que tu fais dans chaque ligne? (nottament la première, que veut dire column + line?)
Je te remercie par avance.

27/09/2015, 23h19
hinaneee

Bonsoir,
peux-tu m'indiquer s'il te plaît, comment on définit le constructeur de copie pour une classe matrice? (sns utiliser dsvp).
Merci par avance.

Code:

1
2
3
4
5
6
7
8
9
namespace dvp
{
template<class T, size_t N>
class matrix
{
...
   matrix(const matrix& m){...};
...
};