kernelopts(assertlevel=1); ############################################################################ # fonction qui calcul une suite de Sturn de f par rapport à la variable var# # !!! f doit être sans facteur carré # ############################################################################ StHaSeq := proc(f,var) local p, SH; # p est la dérivée de f par rapport à var p := diff(f,var); # On initialise la suite SH := [f,p]; while(degree(p,var) <> 0) do p := -rem(SH[nops(SH)-1],SH[nops(SH)],var); SH := [op(SH),p]; od; RETURN(SH); end proc: ############################################## # fonction qui découpe un intervalle en deux # ############################################## Split := J -> [[J[1],(J[1]+J[2])/2],[(J[1]+J[2])/2,J[2]]]: ################################################# # Calcul de la variation d'une liste de nombres # ################################################# Variation := proc(L) local c,v,i,j; v := 0; j := 1; while ( (L[j] = 0 ) and (j < nops(L)) ) do j := j+1; od; c := L[j]; for i from j+1 to nops(L) do if (c*L[i] < 0) then v := v+1; c := L[i]; end if; end do; RETURN(v); end proc: ##################################################################### # fonction permettant d'évaluer la suite de Sturm à la valeur point # ##################################################################### EvalSH := proc(SH,var,point) RETURN([seq(eval(SH[i],var=point),i=1..nops(SH))]); end proc: ########################################################################### # fonction permettant de calculer le nombre de racines dans un intervalle # # !!!!! f = SH[0] ne doit pas s'annuler aux bornes de l'intervalle # ########################################################################### NbRealRoots := proc(SH,var,a,b) local va, vb; if (eval(SH[1],var=a)*eval(SH[1],var=b) = 0) then printf("Bad luck, at least one of the bound of an interval is a root of the polynomial"); # pour détecter qu'il y a une problème RETURN(-1); else va := Variation(EvalSH(SH,var,a)); vb := Variation(EvalSH(SH,var,b)); RETURN(va-vb); end if; end proc: ############################################################ # fonction retournant un intervalle d'isolation par racine # ############################################################ IsolateRealRoots := proc(SH,var,a,b) local L,n; n := NbRealRoots(SH,var,a,b); ASSERT((n <> -1),"Bounds of intervals must not be root"); if (n = 0) then RETURN([]); elif (n = 1) then RETURN([[a,b]]); else L := Split([a,b]); RETURN([ op(IsolateRealRoots(SH,var,L[1][1],L[1][2])), op(IsolateRealRoots(SH,var,L[2][1],L[2][2])) ]); end if; end proc: