Trigonométrie

[CoucouHibou]

Bonjour;

Je suis nouvel utilisateur de XCAS (ou plutôt de QCAS). je cherche à utiliser ce logiciel pour m'aider à calculer des expressions analytiques pour les amplitudes des différents termes de décomposition d'un tenseur quelconque auquel est appliqué une rotation quelconque. Je travaille sous Debian jessie, avec qcas.

Je tombe très rapidement sur de gros soucis pour tout ce qui concerne les expressions trigonométriques. Déjà, lorsque je veux définir une fonction exponentielle complexe, par exemple f(x):=exp(i*x), le rappel de la fonction avec un argument non défini fait disparaître i : pour a non défini dans la feuille de calcul f(a) retourne e^a. Par contre, curieusement, f(pi/2*a) retourne bien e^{i a \frac{pi}{2}}.

Il me suffit donc de "tricher" en n'utilisant que des arguments multiples de pi/2... J'ai donc défini mes matrices de rotation de cette manière, et les ai appliquées à des tenseurs d'ordre 2. Hormis la forme (sous forme de somme d'exponentielle complexe, exp2trig ne donnant rien de bien utile), les résultats me paraissaient justes jusqu'à ce que je tombe sur un autre résultat inexact. Je joins ici les entrées de la feuille de calcul pour que vous puissiez en juger, mais ma question est la suivante : y a-t-il quelque chose que je devrais faire pour éviter ces soucis.

Je précise que je ne suis pas mathématicien, donc il y a sans doute un certain manque de rigueur dans tout ce qui va suivre. Le contexte est le suivant : je veux calculer la modulation de certains termes dans un Hamiltonien d'interaction entre deux spins lorsque j'applique des rotations sur ces spins.

Première entrée : définition des matrices des observables projection du moment cinétique suivant un axe de 2 spins I et S

Code : Tout sélectionner

Iz:=1/2*[[1,0,0,0],[0,1,0,0],[0,0,-1,0],[0,0,0,-1]];
Ix:=1/2*[[0,0,1,0],[0,0,0,1],[1,0,0,0],[0,1,0,0]];
Iy:=1/2*[[0,0,-i,0],[0,0,0,-i],[i,0,0,0],[0,i,0,0]];
Sz:=1/2*[[1,0,0,0],[0,-1,0,0],[0,0,1,0],[0,0,0,-1]];
Sx:=1/2*[[0,1,0,0],[1,0,0,0],[0,0,0,1],[0,0,1,0]];
Sy:=1/2*[[0,-i,0,0],[i,0,0,0],[0,0,0,-i],[0,0,i,0]];

Seconde entrée : définition des tenseurs d'ordre 2 associés :

Code : Tout sélectionner

T20:=sqrt(2/3)*(2*Iz*Sz-Ix*Sx-Iy*Sy);
T2p1:=-1*(Iz*Sx+Ix*Sz+i*(Iz*Sy+Iy*Sz));
T2m1:=(Iz*Sx+Ix*Sz-i*(Iz*Sy+Iy*Sz));
T2p2:=(Ix*Sx-Iy*Sy+i*(Iy*Sx+Ix*Sy));
T2m2:=(Ix*Sx-Iy*Sy-i*(Iy*Sx+Ix*Sy));

Troisième entrée : définition des matrices de passage pour la diagonalisation de Ix+Sx et Iy+Sy

Code : Tout sélectionner

Rp90x:=simplify(eval(exp(-i*pi/2*Ix)))*simplify(eval(exp(-i*pi/2*Sx)));
Rm90x:=simplify(eval(exp(i*pi/2*Ix)))*simplify(eval(exp(i*pi/2*Sx)));
Rp90y:=simplify(eval(exp(-i*pi/2*Iy)))*simplify(eval(exp(-i*pi/2*Sy)));
Rm90y:=simplify(eval(exp(i*pi/2*Iy)))*simplify(eval(exp(i*pi/2*Sy)));

Quatrième entrée : définition des opérateurs de rotation suivant un angle quelconque autour de x, y ou z :

Code : Tout sélectionner

Rpz(c):=[[exp(-i*c),0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,exp(i*c)]];
Rpy(c):=Rm90x*[[exp(-i*c),0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,exp(i*c)]]*Rp90x;
Rpx(c):=Rp90y*[[exp(-i*c),0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,exp(i*c)]]*Rm90y;

À noter ici que : utiliser Rpz(c) pour la définition de Rpx et Rpy ne donne rien de bon à cause du problème d'exponentielle complexe signalé plus haut. L'utilisation de la forme exp(i*c)=cos(c)+i*sin(c) ne fonctionne pas non plus (i est "effacé")

Cinquième entrée : évaluation sur un exemple facile à calculer à la main :

Code : Tout sélectionner

simplify(dotprod(Rpy(c*pi/2)*(Ix*Sx)*Rpy(-c*pi/2),(Ix*Sz+Iz*Sy))

Le résultat de la rotation autour de y de l'opérateur IxSx devrait être cos^2(c) IxSx+sin^2(c) IzSz - sin(c)cos(c) (IxSz+IzSx), or le résultat retourné par cette cinquième entrée est (retranscription complète) :

Code : Tout sélectionner

\frac{(e^{i c pi/2})^4-i}{16(e^{i c pi/2})^2}

ce qui, sauf erreur de ma part n'est pas du tout égal à -sin2c, même s'il y a une petite ressemblance (le facteur 16 n'est pas à prendre en compte, aucun des opérateurs de cet exemple n'étant normalisé)

J'espère que tout ceci n'est qu'une bévue de ma part, mais je vous saurais gré de m'éclairer là dessus, cordialement.

[frederic han]

bonjour,

pour le : f(x):=exp(i*x)
le i manquant dans f(a) est un probleme d'affichage specifique a qcas. (si l'on fait click droit copier sur la reponse et qu'on le colle dans une entree le i est bien present.

Je vais regarder, ca vient surement d'une modif recente dans la conversion mathml.

Frederic

[parisse]

Je trouve -1/8*sin(c*pi) en faisant

Code : Tout sélectionner

normal(exp2trig(lin(simplify(dotprod(Rpy(c*pi/2)*(Ix*Sx)*Rpy(-c*pi/2),(Ix*Sz+Iz*Sy)))

Sinon, je n'observe rien d'anormal si je fais f(x):=exp(i*x) puis f(a) dans Xcas.

[CoucouHibou]

En effet, je m'étais laissé séduire par l'aspect graphique plus moderne de Qcas, mais Xcas n'a aucun des bugs signalés. Par contre, il faut faire attention et signaler que dans qcas, en plus du bug sus-décrit, sqrt(-1) est noté i, alors que dans xcas, c'est I.

Ceci dit, je n'arrive pas simplement à lui faire rassembler les expressions exponentielles complexes en sin et cos : l'utilisation imbriquée de exp2trig et tcollect me donne des fractions avec un polynôme en sin et cos au dénominateur, et non au numérateur.

Merci pour votre aide rapide en tout cas, cordialement,

Hibou

Edit : oups je n'avais pas fais gaffe au code donné ci-dessus pour la simplification, ça marche !

[frederic han]

Donc je confirme, d'ou venait le probleme:
Pour l'affichage qcas utilise la fonction mathml de la librairie giac. Les modifications recentes que j'avais demandé avaient un petit Pb qui faisait disparaitre des i.

J'ai envoye une nouvelle version a bernard pour qu'il mette a jour giac.
a bientot

Frederic

[frederic han]

Voila, Bernard vient de me dire qu'il integre la correction dans la prochaine version de giac (1.1.0-48)
il faudra donc mettre a jour votre paquet giac lorsqu'il sera disponible pour que cela corrige le probleme d'affichage dans qcas.

a bientot

NB: Pour l'affichage de i pour sqrt(-1) dans le mode par defaut de calcul (mode xcas) le sqrt(-1) est i alors que c'est I en mode maple. L'affichage se fait donc avec cette logique. J'ai l'impression que si vous travaillez avec des I c'est que vous etes en mode maple non?

Frederic