[Zangdar MortPartout]

Bonjour,

J'ai à résoudre le type d'équation suivante:

-x^2 - x + (z-1) = 0

C'est très facile avec la toolbox symbolique mais je voudrais le faire sans. Concrètement, j'ai 4 polynômes p1,...p4 en z qui se représentent sous la forme de vecteurs [a(n), a(n-1), ... , a1, a0] ... et l'équation que je cherche à résoudre est:

-x^2*p1 + x*(p2-p3) + p4 = 0

Encore une fois avec la toolbox symbolique le problème est très simple, mais je voudrais faire sans. Avez vous des idées ?

[Eden Paradise]

Citation :

Publié par Zangdar MortPartout

Bonjour,

J'ai à résoudre le type d'équation suivante:

-x^2 - x + (z-1) = 0

C'est très facile avec la toolbox symbolique mais je voudrais le faire sans. Concrètement, j'ai 4 polynômes p1,...p4 en z qui se représentent sous la forme de vecteurs [a(n), a(n-1), ... , a1, a0] ... et l'équation que je cherche à résoudre est:

-x^2*p1 + x*(p2-p3) + p4 = 0

Encore une fois avec la toolbox symbolique le problème est très simple, mais je voudrais faire sans. Avez vous des idées ?

Sans pour quelle raison ?
Parce que le type en sortie c'est du symb ? Si c'est juste ça, tu peux utiliser la conversion double tout simplement.

Sinon, tu devrais expliciter pour quelle raison tu n'as pas accès à toolbox... Car devoir recoder un truc déjà existant c'est tout de même dommage.

Si n<=3 tu as accès aux formules littérales, de Cardan, de Ferrari et de jesaisplusqui. Après quoi Galois t'annonce que cette méthode est terminée.

[Zangdar MortPartout]

Je recode pas, je suis en train de coder. Et sinon, le code avec la symbolic toolbox c'est 3 lignes.

Je veux m'en passer parce que :

1 - je veux que mon programme soit utilisable sur n'importe quel matlab
2 - pour tracer des courbes, l'évaluation des valeurs à partir de subs est beaucoup plus lente qu'avec de simples vecteurs numériques.

Exemple:

Code:

wt=0:.0001:2*pi;
x=exp(i*wt);

tic
FE_1 = (sqrt(2*x-1) - 1); 
FE_2 = (-sqrt(2*x-1) - 1); 
toc

z=sym('z');
FEz_1 = (sqrt(2*z-1) - 1); 
FEz_2 = (-sqrt(2*z-1) - 1); 

tic
FE_1 = double(subs(FEz_1,x));
FE_2 = double(subs(FEz_2,x));
toc

Citation :

Elapsed time is 0.000517 seconds.
Elapsed time is 13.765495 seconds.

Vu que je risque de travailler avec des fonctions d'ordre 4 - 5 ça va vraiment être galère.

Qu'est-ce que la conversion double ?

[Zangdar MortPartout]

Sinon j'ai vu que pour tracer une fonction symbolique il y a ezplot qui marche bien.

Mais ça ne donne pas le résultat attendu, car je dois tracer la partie réelle en abscice, et la partie imaginaire en ordonnées. ezplot me donne la norme en fonction de z.

[Ecthelion]

Ton polynôme en z est-il fixé pour chaque cas dont tu cherches la solution ou bien cherches-tu les solutions telles-que f(x,z)=0 ?

[kermo]

De ce que je comprends il faut résoudre l'équation en x et z.

Il y a une infinité de solutions a priori, et tu veux afficher le lieu des racines en z pendant que x varie entre exp(-i pi) et exp(+i pi) ?

Edit : si le problème est de pouvoir distribuer du code, tu peux le faire en Python

Code:

import sympy as sp
import pylab as pl
import time
from sympy.utilities.lambdify import lambdify

x = sp.var('x')
z = sp.var('z')

p1 = z**2+3*z-1
p2 = z+1
p3 = -2*z+3
p4 = 2*z**2-8*z+2

P = -x**2*p1 + x*(p2-p3) + p4

t = time.time()
sol = sp.solve_poly_system([P], z)

Z = [lambdify(x, s[0], "numpy") for s in sol]
X = pl.exp(1j*pl.arange(0,2*pl.pi, 0.001))

pl.close('all')
for expr in Z:
    v = expr(X)
    pl.plot(v.real, v.imag)
pl.show()
print'Elapsed time', time.time() - t

Là ça affiche le tout en 0.27 sec.

si le problème est le temps de calcul tu peux diminuer le pas, 0.0001 ça me semble très petit pour un lieu des racines.

[Zangdar MortPartout]

C'est exactement ça Kermo, je trace le lieu de racines en fonction de Z ! (c'est z qui varie entre 1 et exp(i*2pi), pas X, donc on résous le polynome en X, puis on trace les racines en faisant varier z)

Merci, idéalement, j'aimerais le faire en matlab pur s'il y a une possibilité de le faire.

Merci

[kermo]

Ok donc c'est l'inverse de ce que j'ai fait, mais ça marche pareil.

Mais pour le coup si tu ne résous qu'en x alors tu n'as qu'un polynôme de degré 2, non ? Tu peux facilement faire une fonction numérique qui te donne la solution en fonction de la valeur des polynômes en un z donné, et du coup pas besoin de symbolic toolbox.

Ça et un pas un peu plus grand et ça tournera sans trop de souci je pense.
Sinon c'est le moment idéal pour comparer avec Sympy ou bien chercher des moyens d'exporter du code Matlab, mais je pense qu'il faut quand même des librairies sur la machine cible.

[Zangdar MortPartout]

Effectivement je ne résous qu'un polynôme de degré 2, donc je pourrais faire varier z puis résoudre.

[Ecthelion]

En matlab, tu as la fonction roots qui, je pense, fais ce que tu veux. Il te faudrait ensuite faire une boucle sur z. Après selon la précision, il est possible qu'il ne soit pas judicieux du tout de faire une boucle.

[Zangdar MortPartout]

Il faut aussi que je puisse ordonner les racines correctement. Je ne sais pas comment fonctionne roots exactement mais il faudrait que les racines soient renvoyées toujours dans le même ordre (du genre -b + d puis -b - d et jamais d'inversion).
Aussi l'idéal serait que je puisse trier le vecteur résultant en fonction des parties réelles croissantes, de façon à éviter les lignes de jointures indésirables sur le graphe.

[kermo]

Citation :

Publié par Zangdar MortPartout

Aussi l'idéal serait que je puisse trier le vecteur résultant en fonction des parties réelles croissantes, de façon à éviter les lignes de jointures indésirables sur le graphe.

Pour ça il suffit de faire un plot en points
Il doit y avoir des fonctions sort() ou équivalent sous Matlab.

[Zangdar MortPartout]

Oui c'est pas très compliqué:

[sorted, idx] = sort(real(a));

idx c'est la liste ordonnées des index, ce qui règle le problème.

[Melchiorus]

Pour résoudre l'équation une fois que tu as donné une valeur à z, tu peux comparer avec des formules explicites pour voir si c'est plus rapide.

Si ton équation est -x^2*p1(z) + x*(p2(z)-p3(z)) + p4(z) = 0 alors :

Pose a=(p2(z)-p3(z))/(2*p1(z)).

Puis b=sqrt(a^2+p4(z)/(p1(z)). (Matlab sait trouver une racine carré pour un nombre complexe)

Les solutions sont x1=a-b et x2=a+b.

[Zangdar MortPartout]

bah ça c'est pas un problème il suffit d'utiliser la fonction roots.
Leur méthode subs par contre laisse clairement à désirer.

[Melchiorus]

Citation :

Publié par Zangdar MortPartout

bah ça c'est pas un problème il suffit d'utiliser la fonction roots.

Je parlais de comparer avec la fonction roots.
J'ose espérer qu'ils ont optimisé leur fonction roots quand c'est un polynôme de degré 2 mais on ne sait jamais.

[Ecthelion]

Du coup, avec roots et sort, c'est bon non, ton problème doit être résolu ?

Sinon, côté optimisation de temps, il pourrait être intéressant de coder toi même des résolutions analytiques plutôt que d'utiliser roots. Comme intégrer la méthode de Ferrari jusqu'à l'ordre 4. A voir si c'est utile de passer du temps à coder cela pour gagner peut-être très peu au final.

[Zangdar MortPartout]

Ça marche, le timing est pas idéal mais l'avantage c'est qu'il est indépendant de la complexité du problème puisque ce qui prend du temps c'est de calculer les racines, en gros 7 secondes par requête.

Par contre ce qui me gène plus c'est le graphique que je sors:

- Le premier graphe c'est la façon "artisanale" de faire, on détermine les racines en z, et ensuite on fait un petit script dans laquelle on passe l'expression des racines, que l'on a écrite soit même. Ensuite, on fait le plot(real(r),imag(r))



- Le second graphe, c'est avec la méthode ou on boucle sur z pour sortir les valeurs de chaque racine.



Idéalement, je veux juste le contour. Je pense que c'est tout rempli parce que c'est la façon dont les racines sont triées. Ce sont surement juste des lignes qui s'entrecroisent et qui donnent cet effet. Donc il me faudrait revoir le tri. Bien sur je ne peux pas garder un tracé rempli puisque je vais placer des pôles à l'intérieur. En plus, à cause des jointures, il est faux.

PS: Si je ne trie pas les racines, j'obtiens le bon tracé, mais en plein. J'ai pas encore essayé avec juste des points mais ça peut donner l'illusion effectivement.

[kermo]

Ça remplit car en triant sur la partie réelle tu oscilles entre les points conjugués, donc t'arrêtes pas de faire des lignes verticales.

Les racines sont continues en fonction de z sauf en quelques points (d'où tes droites). Sauf à trouver ce point de discontinuité, faut mieux plotter sous forme de points pour ne pas avoir l'artefact.

[Zangdar MortPartout]

Merci du coup de main, effectivement avec les points pas besoin de se prendre la tête ni de trier ...

Plus besoin d'intervention manuelle pour sortir la courbe je vais donc pouvoir l'interfacer avec le process que j'ai déjà et qui place les pôles dans la région déterminée ... ça va me simplifier la vie

[kermo]

Cool, c'est pour quoi au fait, du développement ou un cours d'automatique ?

Je suis surpris en tout cas du temps que ça prend.
Je pensais Python plus lent que Matlab sur ce genre de truc, et même en passant par le module symbolique pour ensuite traduire les solutions en numérique c'est largement inférieur à la seconde.

[Zangdar MortPartout]

La fonction subs n'a pas l'air vraiment performante, alors que pourtant ça ne doit pas être très compliqué de générer du C à partir de l'expression littérale ... ou alors je m'en sers mal je sais pas.

Sinon c'est un cours purement théorique sur les techniques de simulation. Donc on n'est pas vraiment dans de l'automatique. C'est plus pour savoir comment choisir sa technique de numérisation d'un système, être capables de différencier ODE45 et ODE3 etc etc ...

[kermo]

Tu fais aussi le lien avec les approches discrètes en z ?
Sous Python il y a une fonction lambdify qui semble générer du code à partir d'une expression, j'avais comparé avec subs et il n'y a pas photo.

[Zangdar MortPartout]

Oui justement, on va du plus simple au plus complexe, on commence par simplement discrétiser à l'entrée puis à la sortie, jusqu'aux méthodes de Runge-Kutta d'ordre 6.

C'est un bon cours car celui que j'avais eu en France en 3ème année était centré uniquement sur l'aspect mathématique et j'avais été complètement largué. Je n'ai rien retenu de ce cours sur les "techniques d'intégration numériques".
Là c'est l'inverse, l'aspect mathématique est totalement "out" et on ne démontre rien, mais par contre je saurais m'en servir au taf. Il devrait bien exister un juste milieu ...