fsolve

[Guillaume]

Bonjour,

j'ai un petit problème avec fsolve que je ne rencontrais pas avant :

g(x):=-exp(2-x)-x*(-(exp(2-x)))+1/x

fsolve(g(x),x,2.1,newton_solver)

renvoie une approximation de l'une ou l'autre des 2 solutions ou undef ou 5.90931518923+4.54385779848*i
ou

terminate called after throwing an instance of 'std::bad_alloc'
what(): std::bad_alloc
Abandon

[parisse]

je n'ai pas de plantage, j'obtiens des racines complexes, ou 0.58... ou undef. C'est normal que newton ne renvoie pas toujours la même chose, étant donné le point de départ, la 1ère recherche échoue et la suivante utilise le générateur aléatoire pour déterminer un nouveau point de départ. Pour pouvoir reproduire l'abandon, il faudrait connaitre la valeur de seed du generateur aleatoire (il suffit de faire rand() juste avant la commande fsolve pour le connaitre).

[Guillaume]

Il a fallu une quinzaine d'essais :

Code : Tout sélectionner

31>> rand()
2108070049

// Time 0
32>> fsolve(g(x),x,2.1,newton_solver)
terminate called after throwing an instance of 'std::bad_alloc'
  what():  std::bad_alloc
Abandon

[parisse]

je viens d'essayer
for j from 1 to 15 do rand() od
rand() me renvoie alors 2108070049
puis g(x):=-exp(2-x)-x*(-(exp(2-x)))+1/x
fsolve(g(x),x,2.1,newton_solver) me renvoie alors 5.90931518923+4.54385779848*i
Vous avez une config spéciale par rapport à la config par défaut?

[Guillaume]

C'est bizarre de tomber exactement sur la même "graine", non ?
Je n'ai pas de config spéciale. Un autre utilisateur m'avait signalé ce problème.

[parisse]

non, c'est normal car je ne fais pas d'initialisation du générateur aléatoire à l'ouverture de session (on peut bien sur le faire explicitement avec srand(NULL) par exemple).

[Guillaume]

Bonjour,

je remplace newton par steffenson mais de nouveaux problèmes arrivent avec for ... in ... do ...end_for

Code : Tout sélectionner

TVapp(L,F,f):={
local  k,j,S,m,fp;

f:=unapply(f,x);
fp:=function_diff(f);

S:=NULL;

j:=[seq(k,k=floor(L[0])..100)] minus F;
  for k in j do for(m:=-5;m<=5;m++){ S:=S,resoudre_numerique(fp(x),x,k+m*0.1,steffenson_solver)};end_for

return(S)
}:;

Alors TVapp([0,+infinity],[0],ln(x))

ne fonctionne pas.

Quand je fais ligne par ligne avec giac ça marche

Code : Tout sélectionner

26>> L:=[0,+infinity];F:=[0];S:=NULL;
[0,+(infinity)],[0]

// Time 0
27>> j:=[seq(k,k=floor(L[0])..100)] minus F:;
"Done"

// Time 0
28>> for k in j do for(m:=-5;m<=5;m++){ S:=S,resoudre_numerique(ln(x),x,1+m*0.1,steffenson_solver)};end_for;

[parisse]

je suppose que l'appel à la GSL est abusé par le changement de signe en x=0. Je pense qu'il vaudrait mieux utiliser une méthode à priori moins efficace de résolution, par exemple dichotomie (bisection_solver) si un changement de signe est repéré entre 2 points de la discrétisation. Les méthodes efficaces ont en effet un inconvénient c'est qu'il faut déjà être assez prêt d'une solution pour que ça converge (et ça converge alors très vite).

[Guillaume]

Ce qui me trouble c'est que la même action mise dans une procédure bloque alors qu'elle fonctionne en dehors (dans l'une, un bad argument value et pas dans l'autre) : cela semble indépendant de la méthode de résolution.

[parisse]

oui, bizarre... Je n'ai pas observé le bad arg value (mais je n'ai pas testé avec la 0.8.5, seulement avec la 0.9.0)

[Guillaume]

N'y a-t-il pas de version 0.9 disponible ?....

[parisse]

Seulement en source. Comme il contient une implementation de l'algorithme de PGCD modulaire pour les extensions algebriques, j'attend que le code soit bien stabilise avant d'en faire des versions binaires.

[parisse]

je viens de tester avec la 0.8.5, et j'obtiens en reponse une liste de 0.0. Il y a donc bien quelque chose de mysterieux!