Intégrales et substitution avec giac

[chowy]

bonjour,

avant de commencer, comme c'est mon premier post, je tenais à vous féliciter pour votre projet... et j'espère pouvoir m'en servir dans mes progs...

mais pour cela, j'aurais besoin de calculer des intégrales simples (pour le moment, doubles par la suite) sur des polynomes, sur un intervalle borné... j'ai cherché dans les sources mais je n'ai pas réussi à trouver de méthode comme dans xcas, c'est à dire sous la forme :

Code : Tout sélectionner

integrate(polynome, x, a, b)

avec [a,b] mon intervalle et x ma variable...

je n'ai trouvé que :

Code : Tout sélectionner

integrate(polynome, x, c)

avec c le contexte. Cette méthode me renvoie alors la primitive de mon polynome.

Mes questions sont donc les suivantes :

1) existe-t-il une méthode "integrate(polynome, x, a, b)" que je n'aurais pas vu ?

2) actuellement, je calcule mon produit scalaire entre deux polynomes avec le code suivant :

Code : Tout sélectionner

giac::identificateur x("x");
giac::gen p1 = 15 * x * x + 25 * x;  // exemple de polynome 
giac::gen p2 = 23 * x + 5; // exemple de polynome

giac::gen primitive_p1_p2 = giac::integrate(p1 * p2, x, 0);
giac::gen integrale_p1_p2 = giac::subst(primitive_p1_p2, x, 1) - giac::subst(primitive_p1_p2, x, -1); // pour calculer sur l'intervalle [-1;1] : F(1) - F(-1) avec F(x) la primitive de mom polynome

std::cout << "Integrale : " << integrale_p1_p2 << std::endl;

ai-je complétement faux ? Parce que certains résultats semblent incorrects (mais seulement pour certains calculs)...

3) A quoi sert le contexte ? j'ai également essayé avec :

Code : Tout sélectionner

giac::context *c;
integrate(polynome, x, c);

sans plus de résultats... faut-il initialiser le "context" d'une certaine façon ?

Voila pour le moment...

D'avance, merci.

[parisse]

bonjour et bienvenue sur ce forum!

en general vous pouvez appeler une fonction xcas en C++ en utilisant le meme nom precede de _ et en groupant les arguments Xcas dans un vecteur, par exemple avec
using namespace giac;
_integrate(makevecteur(p1*p2,x,-1,1),0);
Bon il y a malheureusement plein de fonctions que je n'ai pas exportees dans les .h, donc il peut etre necessaire de rajouter la declaration avant utilisation.
Il y a aussi des fonctions C++ appelees par les fonctions avec _, dans votre cas integrate, qui calcule une primitive. Votre code m'a l'air correct, avez-vous des exemples ou ca ne marche pas?

Quant à la variable de contexte, elle servira (lorsque tout le code sera compatible) a pouvoir executer des threads en parallele sans avoir d'interactions.

[chowy]

bonjour,

effectivement, cela fonctionne avec un "makevecteur"...

par contre, j'ai des problèmes pour calculer des intégrales doubles. J'essaye de faire :

Code : Tout sélectionner

giac::identificateur x1("x1");
giac::identificateur x2("x2");

poly = 3*x1*x2 + x1 + 1; // exemple de polynome bivariable
poly_integrale = giac::_integrate(makevecteur(giac::_integrate(makevecteur(poly, x2, -1, 1), 0), x1, -1, 1), 0); // integrale double sur [-1,1]

ce code fournit le bon résultat dans certains cas, mais pas toujours...

en effet, lorsque "poly" est issu de différents calculs, certains résultats semblent incohérents. Par exemple, dans mon programme, j'ai deux polynomes :

Code : Tout sélectionner

p1 = (3/2*x1*x1*2-1)/2
p2 = (3/2*x2*x2*2-1)/2

les résultats des calculs des intégrales doubles des polynomes au carré sont :

Code : Tout sélectionner

4/5 // OK
4 // ???

de plus, si j'effectue les intégrales dans l'autre sens (d'abors "x1" puis "x2"), les résultats sont inversés.

Y a-t-il certaines précautions à prendre pour effectuer des calculs "enchainés" ? Existe-t-il un moyen direct de calculer ces intégrales doubles ?

merci

[parisse]

bon, j'ai fait un test avec p1^2 et p2^2 et j'obtiens 4/5 pour les 2, mais j'ai du transformer le 3/2*x1*x1*2 car sinon il y a erreur de compilation. De toutes facons, 3/2 poserait probleme a cause de la division qui serait interpretee comme division d'entiers et renverrait 1. J'obtiens 4, 4/5, 4/5, il ne semble pas yy avoir d'erreurs.

Code : Tout sélectionner

#include <iostream>
#include <giac/giac.h>

using namespace std;
using namespace giac;


int main(){
  identificateur x1("x1");
  identificateur x2("x2");
  
  gen poly = 3*x1*x2 + x1 + 1; // exemple de polynome bivariable
  gen poly_integrale = _integrate(makevecteur(_integrate(makevecteur(poly, x2, -1, 1), 0), x1, -1, 1), 0); // integrale double sur [-1,1]
  cout << poly_integrale << endl;
  gen a(rdiv(3,2));
  poly = (2*a*x1*x1-1)/2 ;
  poly_integrale = _integrate(makevecteur(_integrate(makevecteur(poly*poly, x2, -1, 1), 0), x1, -1, 1), 0); // integrale double sur [-1,1]
  cout << poly_integrale << endl;
  poly = (2*a*x2*x2-1)/2;
  poly_integrale = _integrate(makevecteur(_integrate(makevecteur(poly*poly, x2, -1, 1), 0), x1, -1, 1), 0); // integrale double sur [-1,1]
  cout << poly_integrale << endl;
}

[chowy]

bonjour,

voici mon programme qui me pose probleme (désolé, je n'ai pas de site où je pourrais le mettre en téléchargement, mais un simple copier-coller dans le fichier "baseLegendre2D.cpp" suffit...) :

pour info : la classe "IndexKey2D" n'est pas importante, elle sert juste à "simuler" un tableau 2D à partir de std::map...

Code : Tout sélectionner

#include <iostream>
#include <giac/giac.h>

#include <map>
#include <exception>
#include <string>

/**
 * Classe permettant simplement de gérer un "std::map" comme un tableau
 * à deux dimensions. Le paramètre "degree" permet de connaitre la taille
 * du tableau et donc de verifier la validité des index.
 *
 * utilisation :
 *
 * std::map<IndexKey2D, int> array;
 *
 * IndexKey2D index(3);
 *
 * array[index(0, 0)] = 5;
 * int value = array[index(1, 1)];
 */
class IndexKey2D
{
public:
  IndexKey2D(size_t degree)
    :
    m_degree(degree)
  {
    m_indexI = 0;
    m_indexJ = 0;
  }

  IndexKey2D(const IndexKey2D& indexKey2D)
    :
    m_degree(indexKey2D.getDegree()),
    m_indexI(indexKey2D.getIndexI()),
    m_indexJ(indexKey2D.getIndexJ())
  { }

  virtual ~IndexKey2D()
  {
    m_degree = 0; m_indexI = 0; m_indexJ = 0;
  }

  const IndexKey2D& IndexKey2D::operator = (const IndexKey2D& indexKey2D)
  {
    if(this != &indexKey2D)
      {
	m_degree = indexKey2D.getDegree();
	m_indexI = indexKey2D.getIndexI();
	m_indexJ = indexKey2D.getIndexJ();
      }
    
    return *this;
  }
  
  IndexKey2D& operator() (size_t indexI, size_t indexJ)
  {
    setIndexes(indexI, indexJ);
    return (*this);
  }

  IndexKey2D operator() (size_t indexI, size_t indexJ) const
  {
    IndexKey2D res(*this);

    res.setIndexes(indexI, indexJ);

    return res;
  }

  size_t getDegree(void) const { return m_degree; }
  size_t getIndexI(void) const { return m_indexI; }
  size_t getIndexJ(void) const { return m_indexJ; }

  void setIndexes(size_t indexI, size_t indexJ)
  {
    if((indexI + indexJ) <= m_degree)
      {
	m_indexI = indexI;
	m_indexJ = indexJ;
      }
    else
      {
	throw std::exception();
      }
  }

  bool operator < (const IndexKey2D& key) const
  {
    return (m_indexI * (m_degree + 1) + m_indexJ)
      < (key.getIndexI() * (m_degree + 1) + key.getIndexJ());
  }

  friend std::ostream& operator << (std::ostream& os,
				    const IndexKey2D& indexKey2D)
  {
    os << "(" << indexKey2D.getIndexI()
       << "," << indexKey2D.getIndexJ() << ")"; 

    return os;
  }

protected:
  size_t m_degree;
  size_t m_indexI, m_indexJ;
};




/**
 * Génère l'ensemble des polynomes bivariables de Legendre
 * jusqu'au degré "degre".
 *
 * Ces polynomes sont rangés dans une "base" (un std::map")
 * de la façon suivante (pour une base de degré 3 par exemple) :
 *
 *  P0,0   P0,1   P0,2   P0,3
 *  P1,0   P1,1   P1,2    --
 *  P2,0   P2,1    --     --
 *  P3,0    --     --     --
 *  
 *
 * Génération des polynomes de Legendre :
 * P0,0 = 1
 * P1,0 = x1
 * P0,1 = x2
 * Pi,j = (2*i - 1)/i * x1 * Pi-1,0 - (i-1)/i * Pi-2,0
 *
 * Fonction de poids associée : w = 1
 * Domaine de définition : [-1 ; 1]
 *
 * Remarque : la fonction de poids a volontairement été ignorée
 */
std::map<IndexKey2D, giac::gen>
genererBaseLegendre2D(const giac::identificateur& x1,
		      const giac::identificateur& x2,
		      size_t degre)
{
  giac::gen p;

  // Les polynomes sont rangés de la façon suivante :
  // map[0] = P0,0
  // map[1] = P0,1
  // map[2] = P0,2
  // map[i * (degre + 1) + j] = Pi,j
  IndexKey2D indexKey2D(degre);
  std::map<IndexKey2D, giac::gen> lstPoly;
  lstPoly.clear();

  // P0,0 = 1
  p = 1;
  lstPoly[indexKey2D(0, 0)] = p;

  if(degre == 0)
    return lstPoly;

  // P1,0 = x1
  p = x1;
  lstPoly[indexKey2D(1, 0)] = p;

  // P0,1 = x2
  p = x2;
  lstPoly[indexKey2D(0, 1)] = p;

  if(degre == 1)
    return lstPoly;

  // On génère tous les Px,0 à partir de P0,0 et P1,0
  for(int iI = 2; iI <= degre; iI++)
    {
      // Pi,0 = (2*i - 1)/i * x1 * Pi-1,0 - (i-1)/i * Pi-2,0
      lstPoly[indexKey2D(iI, 0)]
	= giac::rdiv(2 * iI - 1, iI) * x1 * lstPoly[indexKey2D(iI - 1, 0)]
  	- giac::rdiv(iI - 1, iI) * lstPoly[indexKey2D(iI - 2, 0)];
    }


  // On gènère tous les Px,1 à partir des Px,0 générés
  //
  // ATTENTION : comme on est sur la 2eme colonne -> on s'arrete une case avant
  // d'où iX1 <= (degre -1)
  for(int iI = 0; iI <= (degre - 1); iI++)
    {
      lstPoly[indexKey2D(iI, 1)] = x2 * lstPoly[indexKey2D(iI, 0)];
    }


  // On génère le reste à partir des Pi,0 et Pi,1
  for(int iI = 0; iI <= degre; iI++)
    {
      for(int iJ = 2; iJ <= degre; iJ++)
 	{
	  if((iI + iJ) <= degre)
	    {
	      // Pi,j = (2*j - 1)/j * x2 * P0,j-1 - (j-1)/j * P0,j-2
   	      lstPoly[indexKey2D(iI, iJ)]
		= giac::rdiv(2 * iJ - 1, iJ)
		* x2 * lstPoly[indexKey2D(iI, iJ - 1)]
		- giac::rdiv(iJ - 1, iJ) * lstPoly[indexKey2D(iI, iJ - 2)];
	    }
	}
    }

  // PROBLEME SI ON ACTIVE LA SIMPLIFICATION !!!
//    for(std::map<IndexKey2D, giac::gen>::iterator iter = lstPoly.begin();
//        iter != lstPoly.end();
//        iter++)
//      { 
//        (*iter).second = giac::simplify((*iter).second);
//      }
  

  return lstPoly;
}




/**
 * Normalisation des polynomes générés précédemment via la formule
 * de Graam-Schmidt :
 *
 * p_normalisé = p / sqrt(integrale(integrale(p * p, x2, -1, 1), x1, -1, 1))
 *
 */
void normaliserBase2D(std::map<IndexKey2D, giac::gen>& base,
 		      const giac::identificateur& x1,
 		      const giac::identificateur& x2)
{
  // On parcourt tous les polynomes de la base...
  for(std::map<IndexKey2D, giac::gen>::iterator iter = base.begin();
      iter != base.end();
      iter++)
    {   
      // On normalise
      // "(*iter).second" correspond ici au polynome courant
      (*iter).second = giac::rdiv((*iter).second, giac::sqrt(giac::_integrate(makevecteur(giac::_integrate(makevecteur((*iter).second *(*iter).second, x2, -1, 1), 0), x1, -1, 1), 0)));
      
      // PROBLEME SI ON ACTIVE LA SIMPLIFICATION !!!
//      (*iter).second = giac::simplify((*iter).second);
    }
}





/**
 * Affichage de tous les polynomes de la base
 */
void afficherBase2D(std::map<IndexKey2D, giac::gen>& base)
{
  // On parcourt séquentiellement le "std::map"
  for(std::map<IndexKey2D, giac::gen>::iterator iter = base.begin();
      iter != base.end();
      iter++)
    {
      // "(*iter).first" correspond à l'index de type (i,j)
      // "(*iter).second" correspond au polynome courant
      std::cout << "P" << (*iter).first << " : " << (*iter).second
		<< std::endl;
    }
}



/**
 * Teste l'orthonormalité des polynomes générés :
 *
 * < Pi1,j1 | Pi2,j2 > = 1 si i1 = i2 et j1 = j2
 * < Pi1,j1 | Pi2,j2 > = 0 sinon
 *
 * avec < Pi1,j1 | Pi2,j2 >
 *           = integrale(integrale(Pi1,j1 * Pi2,j2, x2, -1, 1), x1, -1, 1)
 *
 * Retourne "true" si la base est orthonormale, "false" sinon
 */
bool testOrthonormaliteBase2D(std::map<IndexKey2D, giac::gen>& base,
			      const giac::identificateur& x1,
			      const giac::identificateur& x2,
			      size_t degre)
{
  bool orthonormale = true;
  size_t cpt1 = 0, cpt2 = 0;

  IndexKey2D indexKey2D(degre);

  giac::gen poly_integrale_X1X2;
  
  for(int i1I = 0; i1I <= degre; i1I++)
    {
      for(int i1J = 0; i1J <= degre; i1J++)
	{
	  for(int i2I = i1I; i2I <= degre; i2I++)
	    {
	      for(int i2J = i1J; i2J <= degre; i2J++)
		{
 		  if(((i1I + i1J) <= degre) && ((i2I + i2J) <= degre))
 		    {
		      poly_integrale_X1X2 = giac::_integrate(makevecteur(giac::_integrate(makevecteur(base[indexKey2D(i1I, i1J)] * base[indexKey2D(i2I, i2J)], x2, -1, 1), 0), x1, -1, 1), 0);
		      
		      std::cout << "< P" << i1I << "," << i1J
				<< " | P" << i2I << "," << i2J << " > = "
				<< poly_integrale_X1X2
				<< std::endl;
	  
		      if((i1I == i2I) && (i1J == i2J))
			{
			  if(poly_integrale_X1X2 != 1)
 			    {
			      std::cerr << "Erreur : < P" 
					<< i1I << "," << i2I << " | "
					<< i1J << "," << i2J << " > = "
					<< poly_integrale_X1X2
					<< std::endl;

			      orthonormale = false;
			    }
			}
		      else
			{
			  if(poly_integrale_X1X2 != 0)
			    {
			      std::cerr << "Erreur : < P" 
					<< i1I << "," << i2I << " | "
					<< i1J << "," << i2J << " > = "
					<< poly_integrale_X1X2
					<< std::endl;

			      orthonormale = false;
			    }
			}
		    }
		}
	    }
	}
    }
  return orthonormale;
}




/**
 * Programme principal.
 *
 * usage :
 *
 * ./baseLegendre2D degre
 */
int main(int argc, char **argv)
{
  if(argc != 2)
    {
      std::cerr << "usage : " << argv[0] << " degre" << std::endl
		<< std::endl
		<< "degre : degre maximal des polynomes"
		<< std::endl;

      return EXIT_FAILURE;
    }

  size_t degre = atoi(argv[1]);

  giac::identificateur x1("x1");
  giac::identificateur x2("x2");

  std::map<IndexKey2D, giac::gen> base2D(genererBaseLegendre2D(x1, x2, degre));
  
  std::cout << "Avant normalisation : " << std::endl;
  afficherBase2D(base2D);
  std::cout << std::endl << std::endl;

  normaliserBase2D(base2D, x1, x2);

  std::cout << "Après normalisation : " << std::endl;
  afficherBase2D(base2D);
  std::cout << std::endl << std::endl;

  if(testOrthonormaliteBase2D(base2D, x1, x2, degre))
    {
      std::cout << "La base 2D de degré " << degre << " est orthonormale"
		<< std::endl;
    }
  else
    {
      std::cout << "La base 2D de degré " << degre << " N'est PAS orthonormale"
		<< std::endl;
    }

  return EXIT_SUCCESS;
}

utilisation :

Code : Tout sélectionner

./baseLegendre2D 5 2> ./res_erreur.txt

le programme calcule alors la base et affiche le produit de tous les polynomes entre eux.

Le probleme est que, dans certains cas, on obtient des "undef" au lieu de "0" ! (on peut les voir en regardant dans le fichier "res_erreur.txt")

Si on active les 2 blocs de simplidfication "giac::simplify" alors les "undef" disparaissent mais des "rootof..." apparaissent dès que l'on monte un peu en degré (par exemple : ./baseLegendre2D 15 2> ./res_erreur.txt)

voila, avez-vous une idée ? est-ce un bug ou une mauvaise utilisation de giac de ma part (au début, j'utilisais des giac::fraction à la place des giac::rdiv mais le résultat semble le même)

merci

[parisse]

je n'ai pas de undef avec la derniere version de giac, res_erreur.txt est vide, et voici le resultat affiche

Code : Tout sélectionner

Avant normalisation : 
P(0,0) : 1
P(0,1) : x2
P(0,2) : 3/2*x2*x2-1/2
P(0,3) : 5/3*x2*(3/2*x2*x2-1/2)-2/3*x2
P(0,4) : 7/4*x2*(5/3*x2*(3/2*x2*x2-1/2)-2/3*x2)-3/4*(3/2*x2*x2-1/2)
P(0,5) : 9/5*x2*(7/4*x2*(5/3*x2*(3/2*x2*x2-1/2)-2/3*x2)-3/4*(3/2*x2*x2-1/2))-4/5*(5/3*x2*(3/2*x2*x2-1/2)-2/3*x2)
P(1,0) : x1
P(1,1) : x2*x1
P(1,2) : 3/2*x2*x2*x1-1/2*x1
P(1,3) : 5/3*x2*(3/2*x2*x2*x1-1/2*x1)-2/3*x2*x1
P(1,4) : 7/4*x2*(5/3*x2*(3/2*x2*x2*x1-1/2*x1)-2/3*x2*x1)-3/4*(3/2*x2*x2*x1-1/2*x1)
P(2,0) : 3/2*x1*x1-1/2
P(2,1) : x2*(3/2*x1*x1-1/2)
P(2,2) : 3/2*x2*x2*(3/2*x1*x1-1/2)-1/2*(3/2*x1*x1-1/2)
P(2,3) : 5/3*x2*(3/2*x2*x2*(3/2*x1*x1-1/2)-1/2*(3/2*x1*x1-1/2))-2/3*x2*(3/2*x1*x1-1/2)
P(3,0) : 5/3*x1*(3/2*x1*x1-1/2)-2/3*x1
P(3,1) : x2*(5/3*x1*(3/2*x1*x1-1/2)-2/3*x1)
P(3,2) : 3/2*x2*x2*(5/3*x1*(3/2*x1*x1-1/2)-2/3*x1)-1/2*(5/3*x1*(3/2*x1*x1-1/2)-2/3*x1)
P(4,0) : 7/4*x1*(5/3*x1*(3/2*x1*x1-1/2)-2/3*x1)-3/4*(3/2*x1*x1-1/2)
P(4,1) : x2*(7/4*x1*(5/3*x1*(3/2*x1*x1-1/2)-2/3*x1)-3/4*(3/2*x1*x1-1/2))
P(5,0) : 9/5*x1*(7/4*x1*(5/3*x1*(3/2*x1*x1-1/2)-2/3*x1)-3/4*(3/2*x1*x1-1/2))-4/5*(5/3*x1*(3/2*x1*x1-1/2)-2/3*x1)


Aprs normalisation : 
P(0,0) : 1/2
P(0,1) : x2/sqrt(4/3)
P(0,2) : (3/2*x2*x2-1/2)/sqrt(4/5)
P(0,3) : (5/3*x2*(3/2*x2*x2-1/2)-2/3*x2)/sqrt(4/7)
P(0,4) : (7/4*x2*(5/3*x2*(3/2*x2*x2-1/2)-2/3*x2)-3/4*(3/2*x2*x2-1/2))/sqrt(4/9)
P(0,5) : (9/5*x2*(7/4*x2*(5/3*x2*(3/2*x2*x2-1/2)-2/3*x2)-3/4*(3/2*x2*x2-1/2))-4/5*(5/3*x2*(3/2*x2*x2-1/2)-2/3*x2))/sqrt(4/11)
P(1,0) : x1/sqrt(4/3)
P(1,1) : x2*x1/sqrt(4/9)
P(1,2) : (3/2*x2*x2*x1-1/2*x1)/sqrt(4/15)
P(1,3) : (5/3*x2*(3/2*x2*x2*x1-1/2*x1)-2/3*x2*x1)/sqrt(4/21)
P(1,4) : (7/4*x2*(5/3*x2*(3/2*x2*x2*x1-1/2*x1)-2/3*x2*x1)-3/4*(3/2*x2*x2*x1-1/2*x1))/sqrt(4/27)
P(2,0) : (3/2*x1*x1-1/2)/sqrt(4/5)
P(2,1) : x2*(3/2*x1*x1-1/2)/sqrt(4/15)
P(2,2) : (3/2*x2*x2*(3/2*x1*x1-1/2)-1/2*(3/2*x1*x1-1/2))/sqrt(4/25)
P(2,3) : (5/3*x2*(3/2*x2*x2*(3/2*x1*x1-1/2)-1/2*(3/2*x1*x1-1/2))-2/3*x2*(3/2*x1*x1-1/2))/sqrt(4/35)
P(3,0) : (5/3*x1*(3/2*x1*x1-1/2)-2/3*x1)/sqrt(4/7)
P(3,1) : x2*(5/3*x1*(3/2*x1*x1-1/2)-2/3*x1)/sqrt(4/21)
P(3,2) : (3/2*x2*x2*(5/3*x1*(3/2*x1*x1-1/2)-2/3*x1)-1/2*(5/3*x1*(3/2*x1*x1-1/2)-2/3*x1))/sqrt(4/35)
P(4,0) : (7/4*x1*(5/3*x1*(3/2*x1*x1-1/2)-2/3*x1)-3/4*(3/2*x1*x1-1/2))/sqrt(4/9)
P(4,1) : x2*(7/4*x1*(5/3*x1*(3/2*x1*x1-1/2)-2/3*x1)-3/4*(3/2*x1*x1-1/2))/sqrt(4/27)
P(5,0) : (9/5*x1*(7/4*x1*(5/3*x1*(3/2*x1*x1-1/2)-2/3*x1)-3/4*(3/2*x1*x1-1/2))-4/5*(5/3*x1*(3/2*x1*x1-1/2)-2/3*x1))/sqrt(4/11)


< P0,0 | P0,0 > = 1
< P0,0 | P0,1 > = 0
< P0,0 | P0,2 > = 0
< P0,0 | P0,3 > = 0
< P0,0 | P0,4 > = 0
< P0,0 | P0,5 > = 0
< P0,0 | P1,0 > = 0
< P0,0 | P1,1 > = 0
< P0,0 | P1,2 > = 0
< P0,0 | P1,3 > = 0
< P0,0 | P1,4 > = 0
< P0,0 | P2,0 > = 0
< P0,0 | P2,1 > = 0
< P0,0 | P2,2 > = 0
< P0,0 | P2,3 > = 0
< P0,0 | P3,0 > = 0
< P0,0 | P3,1 > = 0
< P0,0 | P3,2 > = 0
< P0,0 | P4,0 > = 0
< P0,0 | P4,1 > = 0
< P0,0 | P5,0 > = 0
< P0,1 | P0,1 > = 1
< P0,1 | P0,2 > = 0
< P0,1 | P0,3 > = 0
< P0,1 | P0,4 > = 0
< P0,1 | P0,5 > = 0
< P0,1 | P1,1 > = 0
< P0,1 | P1,2 > = 0
< P0,1 | P1,3 > = 0
< P0,1 | P1,4 > = 0
< P0,1 | P2,1 > = 0
< P0,1 | P2,2 > = 0
< P0,1 | P2,3 > = 0
< P0,1 | P3,1 > = 0
< P0,1 | P3,2 > = 0
< P0,1 | P4,1 > = 0
< P0,2 | P0,2 > = 1
< P0,2 | P0,3 > = 0
< P0,2 | P0,4 > = 0
< P0,2 | P0,5 > = 0
< P0,2 | P1,2 > = 0
< P0,2 | P1,3 > = 0
< P0,2 | P1,4 > = 0
< P0,2 | P2,2 > = 0
< P0,2 | P2,3 > = 0
< P0,2 | P3,2 > = 0
< P0,3 | P0,3 > = 1
< P0,3 | P0,4 > = 0
< P0,3 | P0,5 > = 0
< P0,3 | P1,3 > = 0
< P0,3 | P1,4 > = 0
< P0,3 | P2,3 > = 0
< P0,4 | P0,4 > = 1
< P0,4 | P0,5 > = 0
< P0,4 | P1,4 > = 0
< P0,5 | P0,5 > = 1
< P1,0 | P1,0 > = 1
< P1,0 | P1,1 > = 0
< P1,0 | P1,2 > = 0
< P1,0 | P1,3 > = 0
< P1,0 | P1,4 > = 0
< P1,0 | P2,0 > = 0
< P1,0 | P2,1 > = 0
< P1,0 | P2,2 > = 0
< P1,0 | P2,3 > = 0
< P1,0 | P3,0 > = 0
< P1,0 | P3,1 > = 0
< P1,0 | P3,2 > = 0
< P1,0 | P4,0 > = 0
< P1,0 | P4,1 > = 0
< P1,0 | P5,0 > = 0
< P1,1 | P1,1 > = 1
< P1,1 | P1,2 > = 0
< P1,1 | P1,3 > = 0
< P1,1 | P1,4 > = 0
< P1,1 | P2,1 > = 0
< P1,1 | P2,2 > = 0
< P1,1 | P2,3 > = 0
< P1,1 | P3,1 > = 0
< P1,1 | P3,2 > = 0
< P1,1 | P4,1 > = 0
< P1,2 | P1,2 > = 1
< P1,2 | P1,3 > = 0
< P1,2 | P1,4 > = 0
< P1,2 | P2,2 > = 0
< P1,2 | P2,3 > = 0
< P1,2 | P3,2 > = 0
< P1,3 | P1,3 > = 1
< P1,3 | P1,4 > = 0
< P1,3 | P2,3 > = 0
< P1,4 | P1,4 > = 1
< P2,0 | P2,0 > = 1
< P2,0 | P2,1 > = 0
< P2,0 | P2,2 > = 0
< P2,0 | P2,3 > = 0
< P2,0 | P3,0 > = 0
< P2,0 | P3,1 > = 0
< P2,0 | P3,2 > = 0
< P2,0 | P4,0 > = 0
< P2,0 | P4,1 > = 0
< P2,0 | P5,0 > = 0
< P2,1 | P2,1 > = 1
< P2,1 | P2,2 > = 0
< P2,1 | P2,3 > = 0
< P2,1 | P3,1 > = 0
< P2,1 | P3,2 > = 0
< P2,1 | P4,1 > = 0
< P2,2 | P2,2 > = 1
< P2,2 | P2,3 > = 0
< P2,2 | P3,2 > = 0
< P2,3 | P2,3 > = 1
< P3,0 | P3,0 > = 1
< P3,0 | P3,1 > = 0
< P3,0 | P3,2 > = 0
< P3,0 | P4,0 > = 0
< P3,0 | P4,1 > = 0
< P3,0 | P5,0 > = 0
< P3,1 | P3,1 > = 1
< P3,1 | P3,2 > = 0
< P3,1 | P4,1 > = 0
< P3,2 | P3,2 > = 1
< P4,0 | P4,0 > = 1
< P4,0 | P4,1 > = 0
< P4,0 | P5,0 > = 0
< P4,1 | P4,1 > = 1
< P5,0 | P5,0 > = 1
La base 2D de degr 5 est orthonormale

Donc je vous conseille de mettre a jour avec la derniere version de giac-0.6.2 et de recompiler pour voir si ca marche.
Pour les rootof, c'est normal car giac essaie de tout ecrire avec une unique extension algebrique, ce que vous ne souhaitez probablement pas.
Sinon, a propos des structures de donnees, si vous pensez travailler avec des n grands, vous avez peut-etre interet pour des raisons d'efficacite a utiliser les polynomes de giac, soit les polynomes listes a 1 variable dont les coefficients peuvent etre eux-memes des polynomes listes, soit directement les polynomes a plusieurs variables (mais ils sont d'usage moins evident). Les polynomes listes sont des vecteur (giac::vector<gen>=giac::vecteur) dont les coeffs sont donnes par ordre decroissant. Il existe pas mal de fonctions pour les manipuler dans modpoly.h, par exemple +, *, etc. ou operator_plus, operator_times (utiliser 0 comme pointeur d'environnement pour ne pas faire d'arithmetique modulaire), ou horner pour evaluer en un gen, integrate, derivative, etc. La conversion avec la forme symbolique se fait avec _e2r et _r2e.

[chowy]

bonjour,

effectivement, je n'ai plus de "undef" avec la version 0.6.2 de giac...

par contre, j'ai essayé de convertir mes "giac::gen" en polynomes mais je ne sais pas quelle classe utiliser entre "giac::polynome", "giac::modpoly", "giac::vecteur"...

par exemple, pour un polynome bivariable tel que : 3*x1*x2 - x1 + 2*x2 - 1, que me conseillerez-vous d'utiliser ? en gros, je veux avoir acces aux opérations arithmétiques classiques (+, -, *), au calcul d'intégrales (comme précedemment) et à la substitution de variables (ici x1 ou x2) par une valeur...

Merci

[parisse]

vecteur et modpoly designent le meme type (j'utilise 2 noms pour indiquer des fonctions de type polynomial ou vectoriel), des polynomes denses a 1 variable (cf. la doc succinte de giac.info), par ex. Mais les coeffs de ces polynomes peuvent aussi etre des polynomes a 1 variable, par exemple x2 dans votre cas.
Ainsi pour coder disons x1^4+(3x2^3-5x2+1)*x1+7*x2+5, on creerait un
vecteur p(makevecteur(1,0,0,makevecteur(3,0,-5,1),makevecteur(7,5));
Si vous l'utilisez comme un gen, il faut faire
gen g(p,_POLY1__VECT);
ensuite g*g ou g+g fera les operations necessaires.
Pour extraire le vecteur de g,
if (g.type==_VECT)
p=*g._VECTptr;
Pour integrer integrate(p,1), pour deriver derivate(p).
Pour evaluer p en x1=1, x2=2 on peut faire
_r2e(makevecteur(g,makevecteur(1,2)))

[chowy]

bonjour,

suite à vos conseils, j'essaye de passer les "giac::gen" en "giac::modpoly". Malheureusement, je n'arrive pas à faire certaines choses. Voici mon nouveau programme, le but étant toujours de calculer une base composée des polynomes de Legendre (ici, retour en 1D) :

Code : Tout sélectionner

#include <iostream>
#include <giac/giac.h>

/**
 * Génère l'ensemble des polynomes univariables de Legendre
 * jusqu'au degré "degre".
 *
 * Ces polynomes sont rangés dans un "std::vector"
 * de la façon suivante (pour une base de degré 3 par exemple) :
 *
 *  P0   P1   P2   P3
 *
 *
 * Génération des polynomes de Legendre :
 * P0 = 1
 * P1 = x
 * Pi = (2*i - 1)/i * x * Pi-1 - (i-1)/i * Pi-2
 *
 * Fonction de poids associée : w = 1
 * Domaine de définition : [-1 ; 1]
 *
 * Remarque : la fonction de poids a volontairement été ignorée
 */
std::vector<giac::modpoly> genererBaseLegendre1D(size_t degre)
{
  giac::modpoly p;
  std::vector<giac::modpoly> lstPoly;
  lstPoly.clear();

  // P0 = 1
  p = giac::makevecteur(1);
  lstPoly.push_back(p);

  if(degre == 0)
    return lstPoly;

  // P1 = x
  p = giac::makevecteur(1, 0);
  lstPoly.push_back(p);

  if(degre == 1)
    return lstPoly;

  // On génère le reste à partir des Pi,0 et Pi,1
  for(int i = 2; i <= degre; i++)
    {
      // Pi = (2*i - 1)/i * x * Pi-1 - (i-1)/i * Pi-2
      p = giac::makevecteur(giac::rdiv(2 * i - 1, i))
	* giac::makevecteur(1, 0) * lstPoly[i - 1]
  	- giac::makevecteur(giac::rdiv(i - 1, i)) * lstPoly[i - 2];

      lstPoly.push_back(p);
    }

  return lstPoly;
}



/**
 * Normalisation des polynomes générés précédemment via la formule
 * de Graam-Schmidt :
 *
 * p_normalisé = p / sqrt(integrale(p * p, x, -1, 1))
 *
 */
void normaliserBase1D(std::vector<giac::modpoly>& base)
{
  giac::gen x("x");
  giac::gen genPoly, genPoly_2;
  giac::gen genPoly_symb, genPoly_2_symb;
  giac::gen res;

  // On parcourt tous les polynomes de la base...
  for(std::vector<giac::modpoly>::iterator iter = base.begin();
      iter != base.end();
      iter++)
    {   
      // On normalise
      //      (*iter) = giac::rdiv((*iter), giac::sqrt(giac::_integrate(giac::makevecteur((*iter) * (*iter), 1), 0)));

      genPoly = giac::gen((*iter), giac::_POLY1__VECT);
      genPoly_2 = giac::gen((*iter) * (*iter), giac::_POLY1__VECT);
      genPoly_symb = giac::r2e(genPoly, x);
      genPoly_2_symb = giac::r2e(genPoly_2, x);
      
      res = giac::rdiv(genPoly_symb, giac::sqrt(giac::_integrate(giac::makevecteur(genPoly_2_symb, x, -1, 1), 0)));
     
      std::cout << "Res : " << res << std::endl;
      std::cout << "Poly : " << giac::e2r(res, x) << std::endl;
    }
}



/**
 * Affichage de tous les polynomes de la base
 */
void afficherBase1D(std::vector<giac::modpoly>& base)
{
  // On parcourt séquentiellement la base
  size_t cpt = 0;
  for(std::vector<giac::modpoly>::iterator iter = base.begin();
      iter != base.end();
      iter++)
    {
      std::cout << "P" << cpt << " : " << (*iter) << std::endl;
      cpt++;
    }
}




/**
 * Programme principal.
 *
 * usage :
 *
 * ./baseLegendre1D degre
 */
int main(int argc, char **argv)
{
  if(argc != 2)
    {
      std::cerr << "usage : " << argv[0] << " degre" << std::endl
		<< std::endl
		<< "degre : degre maximal des polynomes"
		<< std::endl;

      return EXIT_FAILURE;
    }

  size_t degre = atoi(argv[1]);


  std::vector<giac::modpoly> base1D(genererBaseLegendre1D(degre));
  
  std::cout << "Avant normalisation : " << std::endl;
  afficherBase1D(base1D);
  std::cout << std::endl << std::endl;

  normaliserBase1D(base1D);

  std::cout << "Après normalisation : " << std::endl;
  afficherBase1D(base1D);
  std::cout << std::endl << std::endl;

  return EXIT_SUCCESS;
}

Utilisation :

Code : Tout sélectionner

./baseLegendre1Dpoly 5

Le probleme se trouve dans la fonction "normaliserBase1D" :

1) je n'arrive pas à calculer (à partir des "giac::modpoly") mon integrale sur [-1;1] ni ma primitve. Peut-on le faire directement ? J'ai essayé un "giac::integrate(p, 1)" mais il me retourne des résultats que je n'arrive pas à interpréter...

2) j'essaye donc de passer par des "giac::gen" via "r2e" et "e2r" mais j'ai des problemes (je ne dois pas faire les choses correctement). Ici, le prog me renvoie un :

Code : Tout sélectionner

terminate called after throwing an instance of 'std::runtime_error'
  what():  Invalid dimension
Aborted

Pourriez-vous m'aiguiller ? Notamment pour la 1ere question, sachant que je souhaiterais utiliser ensuite des polynomes 2D (comme vous avez commencé à m'expliquer précédemment)

Merci

[parisse]

Bon, c'est de ma faute, la fonction integrate des modpoly n'est pas documentee correctement. Elle prend en effet en argument un vecteur mais avec les coeffs dans l'autre sens et elle n'ajoute pas le 0. Donc ca donne

Code : Tout sélectionner

void normaliserBase1D(std::vector<giac::modpoly>& base)
{
  giac::gen x("x");
  giac::gen genPoly, genPoly_2;
  giac::gen genPoly_symb, genPoly_2_symb;
  giac::gen res;

  // On parcourt tous les polynomes de la base...
  for(std::vector<giac::modpoly>::iterator iter = base.begin();
      iter != base.end();
      iter++)
    {   
      // On normalise
      giac::modpoly tmp((*iter) * (*iter));
      reverse(tmp.begin(),tmp.end());
      tmp=giac::integrate(tmp, giac::plus_one);
      reverse(tmp.begin(),tmp.end());      
      tmp.push_back(giac::zero);
      giac::gen t=giac::horner(tmp,1)-giac::horner(tmp,-1);
      (*iter) = inv(giac::sqrt(t)) * (*iter);
    }
}

Pour la generation des polynomes avant normalisation, vous pouvez ecrire un peu plus simplement

Code : Tout sélectionner

  // On gnre le reste  partir des Pi,0 et Pi,1
  for(int i = 2; i <= degre; i++)
    {
      // Pi = (2*i - 1)/i * x * Pi-1 - (i-1)/i * Pi-2
      p = giac::rdiv(2 * i - 1, i) * giac::makevecteur(1, 0) * lstPoly[i - 1]
     - giac::rdiv(i - 1, i) * lstPoly[i - 2];

      lstPoly.push_back(p);
    }

Pour les polynomes de Legendre sans poids, il est encore plus efficace de faire comme dans permu.cc en calculant directement chaque coeff non nul en fonction du precedent non nul, je ne sais pas si ca se generalise.

[chowy]

rebonjour,

mon nouveau programme est donc le suivant :

Code : Tout sélectionner

#include <iostream>
#include <giac/giac.h>

/**
 * Génère l'ensemble des polynomes univariables de Legendre
 * jusqu'au degré "degre".
 *
 * Ces polynomes sont rangés dans un "std::vector"
 * de la façon suivante (pour une base de degré 3 par exemple) :
 *
 *  P0   P1   P2   P3
 *
 *
 * Génération des polynomes de Legendre :
 * P0 = 1
 * P1 = x
 * Pi = (2*i - 1)/i * x * Pi-1 - (i-1)/i * Pi-2
 *
 * Fonction de poids associée : w = 1
 * Domaine de définition : [-1 ; 1]
 *
 * Remarque : la fonction de poids a volontairement été ignorée
 */
std::vector<giac::modpoly> genererBaseLegendre1D(size_t degre)
{
  giac::modpoly p;
  std::vector<giac::modpoly> lstPoly;
  lstPoly.clear();

  // P0 = 1
  p = giac::makevecteur(1);
  lstPoly.push_back(p);

  if(degre == 0)
    return lstPoly;

  // P1 = x
  p = giac::makevecteur(1, 0);
  lstPoly.push_back(p);

  if(degre == 1)
    return lstPoly;

  // On génère le reste à partir des Pi,0 et Pi,1
  for(int i = 2; i <= degre; i++)
    {
      // Pi = (2*i - 1)/i * x * Pi-1 - (i-1)/i * Pi-2
      p = giac::rdiv(2 * i - 1, i) * giac::makevecteur(1, 0) * lstPoly[i - 1]
  	- giac::rdiv(i - 1, i) * lstPoly[i - 2];

      lstPoly.push_back(p);
    }

  return lstPoly;
}



/**
 * Normalisation des polynomes générés précédemment via la formule
 * de Graam-Schmidt :
 *
 * p_normalisé = p / sqrt(integrale(p * p, x, -1, 1))
 *
 */
void normaliserBase1D(std::vector<giac::modpoly>& base)
{
  giac::gen x("x");
  giac::gen genPoly, genPoly_2;
  giac::gen genPoly_symb, genPoly_2_symb;
  giac::gen res;

  // On parcourt tous les polynomes de la base...
  for(std::vector<giac::modpoly>::iterator iter = base.begin();
      iter != base.end();
      iter++)
    {   
      // On calcule le carré du polynome
      giac::modpoly tmp((*iter) * (*iter));

      // On calcule la primitive
      reverse(tmp.begin(), tmp.end());
      tmp = giac::integrate(tmp, giac::plus_one);
      reverse(tmp.begin(), tmp.end());     
      tmp.push_back(giac::zero);

      // On calcule l'intégrale sur [-1;1]
      giac::gen t = giac::horner(tmp, 1) - giac::horner(tmp, -1);

      // On normalise
      //      (*iter) = giac::simplify(giac::inv(giac::sqrt(t))) * (*iter); 
      (*iter) = giac::inv(giac::sqrt(t)) * (*iter); 
    }
}



/**
 * Affichage de tous les polynomes de la base
 */
void afficherBase1D(std::vector<giac::modpoly>& base)
{
  // On parcourt séquentiellement la base
  size_t cpt = 0;
  for(std::vector<giac::modpoly>::iterator iter = base.begin();
      iter != base.end();
      iter++)
    {
      std::cout << "P" << cpt << " : " << (*iter) << std::endl;
      cpt++;
    }
}



/**
 * Teste l'orthonormalité des polynomes générés :
 *
 * < Pi1 | Pi2 > = 1 si i1 = i2
 * < Pi1 | Pi2 > = 0 sinon
 *
 * avec < Pi1 | Pi2 >
 *           = integrale(Pi1 * Pi2, x, -1, 1)
 *
 * Retourne "true" si la base est orthonormale, "false" sinon
 */
bool testOrthonormaliteBase1D(std::vector<giac::modpoly>& base,
			      size_t degre)
{
  bool orthonormale = true;
  size_t cpt1 = 0, cpt2 = 0;

  giac::gen poly_integrale;
  
  for(int i1 = 0; i1 <= degre; i1++)
    {
      for(int i2 = i1; i2 <= degre; i2++)
	{
	  giac::modpoly tmp(base[i1] * base[i2]);

 	  reverse(tmp.begin(),tmp.end());
	  tmp=giac::integrate(tmp, giac::plus_one);
 	  reverse(tmp.begin(),tmp.end());     
 	  tmp.push_back(giac::zero);
 	  giac::gen poly_integrale=giac::horner(tmp,1)-giac::horner(tmp,-1);

 	  std::cout << "< P" << i1 << " | P" << i2 << " > = "
 		    << poly_integrale << "     -> primitive : " << giac::simplify(tmp) << std::endl;

	  
	  if(((i1 == i2) && (poly_integrale != 1))
	     || ((i1 != i2) && (poly_integrale != 0)))
	    {
	      std::cerr << "Erreur : < P" << i1 << " | " << i2 << " > = "
			<< poly_integrale << std::endl;
	      
	      orthonormale = false;
	    }
	}
    }
  
  return orthonormale;
}

/**
 * Programme principal.
 *
 * usage :
 *
 * ./baseLegendre1D degre
 */
int main(int argc, char **argv)
{
  if(argc != 2)
    {
      std::cerr << "usage : " << argv[0] << " degre" << std::endl
		<< std::endl
		<< "degre : degre maximal des polynomes"
		<< std::endl;

      return EXIT_FAILURE;
    }

  size_t degre = atoi(argv[1]);


  std::vector<giac::modpoly> base1D(genererBaseLegendre1D(degre));
  
  std::cout << "Avant normalisation : " << std::endl;
  afficherBase1D(base1D);
  std::cout << std::endl << std::endl;

  normaliserBase1D(base1D);

  std::cout << "Après normalisation : " << std::endl;
  afficherBase1D(base1D);
  std::cout << std::endl << std::endl;

  if(testOrthonormaliteBase1D(base1D, degre))
    {
      std::cout << "La base 1D de degré " << degre << " est orthonormale"
		<< std::endl;
    }
  else
    {
      std::cout << "La base 1D de degré " << degre << " N'est PAS orthonormale"
		<< std::endl;
    }

  return EXIT_SUCCESS;
}

effectivement, les polynomes sont maintenant bien normalisés (je connais les 1er polynomes). Par contre, j'ai des résultats faux lors que je vérifie l'orthonormalité (alors qu'il sont bien orthonormaux !!!). Exemple pour un degré 2 :

Code : Tout sélectionner

Avant normalisation : 
P0 : Vector [1]
P1 : Vector [1,0]
P2 : Vector [3/2,0,-1/2]


Apres normalisation : 
P0 : Vector [1/(sqrt(2))]
P1 : Vector [1/(sqrt(2/3)),0]
P2 : Vector [3/2/sqrt(2/5),0,-1/2/sqrt(2/5)]


< P0 | P0 > = 1     -> primitive : [1/2,0]
< P0 | P1 > = 0     -> primitive : [1/4*sqrt(3),0,0]
< P0 | P2 > = 4*sqrt(5)/16+4*sqrt(5)/16+4*sqrt(5)/16+4*sqrt(5)/16     -> primitive : [1/4*sqrt(5),0,1/4*sqrt(5),0]
Erreur : < P0 | 2 > = 4*sqrt(5)/16+4*sqrt(5)/16+4*sqrt(5)/16+4*sqrt(5)/16
< P1 | P1 > = 1     -> primitive : [1/2,0,0,0]
< P1 | P2 > = 0     -> primitive : [3/16*sqrt(15),0,1/8*sqrt(15),0,0]
< P2 | P2 > = 1     -> primitive : [9/8,0,-5/4,0,5/8,0]
La base 1D de degre 2 N'est PAS orthonormale

1) Ici les résultats sont bons sauf pour <P0|P2> : il semble que la primitive calculée [1/4*sqrt(5),0,1/4*sqrt(5),0] soit incorrecte et devrait être : [1/4*sqrt(5),0, -1/4*sqrt(5),0] (le "-" devant le "1" du 3eme élément) ce qui founirait bien 0.

Bug ou mauvaise utilisation de ma part ? J'ai pourtant repris le même code de calcul d'intégrale que pour la normalisation.

2) Pour info : le "giac::plus_one" dans le "giac::integrate" veut-il dire que, dans le cas de calcul de primitives sur des vecteurs, il faut décaler les coeff d'une case (par exemple pour passer de x à x^2) ?

3) Concernant la génération directe des poids : effectivement, je pourrais les générer plus facilement, mais je ne sais pas trop comment générer ces poids avec des polynomes N-d (en fait, j'ai pas vraiment cherché...), je vais y regarder de plus pres...

Merci

[parisse]

il y a en effet un bug dans giac, je mets a jour le source. Mais vous devrez rajouter un normal pour verifier que le produit scalaire est nul, et le temps de calcul sera long a cause des racines carrees. Pour etre plus efficace, il ne faudrait pas normaliser la base avant la verification mais uniquement apres avoir fait les produits scalaires. Peut-etre devriez-vous envisager de conserver 2 donnees par polynome, le polynome et son coefficient de normalisation.