Trouver rapidement les zéros d'une fonction

Bonjour,

dans le cadre de mon boulot, je cherche à accélérer le calcul du zéro d'une fonction (assez compliquée, je crois que c'est la formule de Manning en hydraulique), car elle est exécutée 5 millions de fois.

Pour l'heure, j'ai réussi à transformer cette fonction en polynome et ainsi rechercher des racines (roots) plutot que des zéros (fzero). Cela devrait accélérer un peu le programme.

Je me demande s'il est possible de vectoriser cette recherche de racines.
Pour info, la forme du polynome (d'ordre 5) dont je cherche les racines est :

Code:

---

P = [M(i)^5*Q(i)^1.5 0 0 2*Q(i) 3*M(i)*Q(i)^2 M(i)]

---

Merci !

Salut!

Citation:

---

j'ai réussi à transformer cette fonction en polynome

---

Sans vouloir te vexer, c'est la pire chose à faire. En effet, les polynômes sont parmi les premières fonctions qu'on apprend à l'école, et on en tire la conclusion fausse que ce sont de "gentilles fonctions"; or, pour tout dire, ce sont plutôt des "sales bêtes" avec lesquelles on a souvent les pires ennuis, du fait de leur instabilité: deux polynômes ayant des coefficients presque identiques peuvent avoir des zéros complètement différents.
Jean-Marc Blanc

Citation:

---

Envoyé par FR119492

Salut!

Sans vouloir te vexer, c'est la pire chose à faire. En effet, les polynômes sont parmi les premières fonctions qu'on apprend à l'école, et on en tire la conclusion fausse que ce sont de "gentilles fonctions"; or, pour tout dire, ce sont plutôt des "sales bêtes" avec lesquelles on a souvent les pires ennuis, du fait de leur instabilité: deux polynômes ayant des coefficients presque identiques peuvent avoir des zéros complètement différents.
Jean-Marc Blanc

---

J'explique ce que j'ai fait :
j'avais une fonction du type :

Code:

---

f(x) = A(i)*x*(x/(2*x+B(i))^(2/3)-Q(i)

---

que je devais annuler.

J'ai juste fait quelques opérations mathématiques de base, pour me retrouver avec la fonction polynomiale :

Code:

---

P(x) = C(i)*x^5 + D(i)*x^2 + E(i)*x + F(i)

---

qui s'annule aux même points que f(x).

Cela, afin de gagner en temps de calcul (roots est apparemment plus rapide que fzero).

Je ne pense pas que cela nuirait à la qualité des solutions (après quelques tests, les résultats sont équivalents).

Je voudrais surtout savoir s'il y a moyen de "vectoriser" roots ;)

Citation:

---

Envoyé par nahouto

Je voudrais surtout savoir s'il y a moyen de "vectoriser" roots ;)

---

Tiré de la documentation :

Citation:

---

The algorithm simply involves computing the eigenvalues of the companion matrix:

Code:

---

1 2 3

A = diag(ones(n-1,1),-1); A(1,:) = -c(2:n+1)./c(1); eig(A)

---

---

Tu peux donc diminuer le temps de calcul de ROOTS en créant ta propre fonction MYROOTS et en négligeant les tests et autres mises en forme de données présents au début de ROOTS.

Code:

---

1 2	edit roots

---

Maintenant, même si tout semble fonctionner, garde bien à l'esprit la remarque de Jean-Marc... on ne sait jamais... au cas où ;)

En suivant l'idée de Dut, j'ai épuré la fonction roots de Matlab, et je parviens à diviser par 2 le temps de calcul.

Mais j'ai particulièrement peur des conséquences des deux tests que j'ai retiré :

Code:

---

1 2 3 4 5 6 7

% c = le polynome à traiter % Strip leading zeros and throw away. % Strip trailing zeros, but remember them as roots at zero. inz = find(c); nnz = length(inz); c = c(inz(1):inz(nnz)); r = zeros(n-inz(nnz),1,class(c));

---

je ne sais pas si ce n'est que de la mise en forme...

Par contre, je ne pense pas tomber sur un polynome mal conditionné et devoir subir ça :

Code:

---

1 2 3 4 5 6 7 8 9

% Prevent relatively small leading coefficients from introducing Inf % by removing them. d = c(2:end)./c(1); ind = find(isinf(abs(d)), 1); while (~isempty(ind)) c = c(ind+1:end); d = c(2:end)./c(1); ind = find(isinf(abs(d)), 1); end

---

Sinon, vous me confirmez que ce genre d'opérations n'est pas vectorisable ?

Je ne suis pas un spécialiste donc les informations suivantes sont à prendre avec précaution... :aie:

Le premier test sert à diminuer la taille du problème à traiter en négligeant les coefficient 0 au début et à la fin du polynôme. Comme dans ton cas, l'ordre du polynôme est fixe, ce n'est pas important.

L'autre test porte sur le coefficient de l'ordre le plus élevé du polynôme. Il teste si ce coefficient n'est pas proche de 0. Dans ce cas, tout dépend de savoir si la valeur de ce coefficient est maîtrisée dans ton cas ou pas ?

Citation:

---

Envoyé par Dut

Je ne suis pas un spécialiste donc les informations suivantes sont à prendre avec précaution... :aie:

Le premier test sert à diminuer la taille du problème à traiter en négligeant les coefficient 0 au début et à la fin du polynôme. Comme dans ton cas, l'ordre du polynôme est fixe, ce n'est pas important.

L'autre test porte sur le coefficient de l'ordre le plus élevé du polynôme. Il teste si ce coefficient n'est pas proche de 0. Dans ce cas, tout dépend de savoir si la valeur de ce coefficient est maîtrisée dans ton cas ou pas ?

---

merci !

Pour le second test, étant donné que je traite des millions de polynomes, je ne sais pas si certains sont mal conditionnés (premier coeff trop petit). Je vais donc garder ce test.

Merci à tous !

Salut!

Quelques réflexions en vrac:

• Ta fonction f(x) présente une asymptote verticale en x=-B(i)/2, alors que P(x) est continue quelque soit x; étrange ...
• Que sais-tu a priori sur tes zéros:
  
  • combien y en a-t-il?
  • lequel ou lesquels cherches-tu?
  • sont-ils positifs, plus ou moins négatifs que -B(i)/2?
• Leibniz a écrit "Natura non facit saltus"; d'une valeur de i à la suivante, la variation des zéros doit être assez modérée. Si tu utilises fzero, tu peux gagner beaucoup en vitesse par un choix approprié de la variable x0, par exemple xo(i+1)=xz(i) ou mieux x0(i+1)=2*xz(i)-xz(i-1), xz étant le zéro trouvé précédemment.

Mais, une fois de plus, c'est de l'algorithmique numérique et plus du Matlab (salut Dut :D)

Jean-Marc Blanc

Salut !

• Je ne l'ai pas précisé, mais touts les variables du problème sont positives, donc pas de problème pour x=-B(i)/2.
• Je cherche un seul zéro, qui correspond à une hauteur d'eau en mètres. Il doit donc être réel positif, et compris entre 0 et 100 environ. Néanmoins, la résolution par roots donne 5 zéros, dont 4 complexes.
• Merci pour cette très bonne idée. Je parcours un cours d'eau, donc il serait en effet logique de supposer que les hauteurs d'eau varient peu le long du cours d'eau. Je vais explorer cette idée.

Je vais voir tout ça, mais en ce moment, je pars sur une autre idée : réaliser une abaque donnant des correspondances entre les valeurs de mes paramètres (M(i) et Q(i)) et le zéro recherché.

Salut!

Citation:

---

Néanmoins, la résolution par roots donne 5 zéros, dont 4 complexes

---

C'est normal, parce que c'est un polynôme du 5ème degré. Est-ce que je ne t'avais pas dit que les polynômes étaient des "sales bêtes"?
Jean-Marc Blanc

Citation:

---

Envoyé par FR119492

Est-ce que je ne t'avais pas dit que les polynômes étaient des "sales bêtes"?
Jean-Marc Blanc

---

si, mais je ne vois pas le rapport :?

Citation:

---

la résolution par roots donne 5 zéros

---

roots ne s'applique qu'à des polynômes.
Jean-Marc Blanc

Citation:

---

Envoyé par FR119492

roots ne s'applique qu'à des polynômes.
Jean-Marc Blanc

---

désolé, mais je ne comprends pas où tu veux en venir.

J'applique la commande roots sur le polynome de degré 5 puis je sélectionne le zéro réel.

Citation:

---

J'applique la commande roots sur le polynome de degré 5 puis je sélectionne le zéro réel.

---

C'est tout à fait correct. Je voulais seulement faire remarquer que tu trouves aussi des zéros qui ne veulent rien dire :P
Jean-Marc Blanc

Bon, je pense que je vais rester sur cette idée de roots, même non vectorisé.

Je vous donnerais des nouvelles après mes simulations.

A tout hasard, 'lsqnonlin' de l'optimization toolbox ? ou encore 'fminunc'

Je me suis pas plongé à fond dans tes équations, mais ces fonctions m'ont bien servis pour des problèmes d'optimisation des fois.

Citation:

---

Envoyé par Pierre845

A tout hasard, 'lsqnonlin' de l'optimization toolbox ? ou encore 'fminunc'

Je me suis pas plongé à fond dans tes équations, mais ces fonctions m'ont bien servis pour des problèmes d'optimisation des fois.

---

Merci, mais je ne dispose pas cette toolbox.

Je ne peux pas encore vous donner le gain en temps de calcul (entre roots et fzero) car mon matlab est planté. J'ai relancé le processus pour cette nuit.

Pour l'idée d'initialiser fzero avec le pas précédent, c'est assez rapide (environ 2 fois plus rapide), mais parfois il me renvoie une erreur :

Code:

---

Exiting fzero: aborting search for an interval containing a sign change because complex function value encountered during search.

---

En effet, lorsque je fais f(solution_précédente), j'obtiens un nombre complexe 8O.

Je pense que je vais rester sur le roots...

Je vous donne quelques nouvelles.

En remplaçant le fzero par un roots, je divise le temps de calcul de cette fonction par 10.
En reprenant la suggestion de Dut de supprimer les tests de roots, je divise le temps de calcul de cette fonction par 20.

Dans les deux cas, les résultats obtenus sont similaires à ceux obtenus précédemment :yaisse2:

Alors, merci à tous ! :king: