1. La probabilité, ramenée à celles des voisins 

La probabilité est la somme des voisines multipliées par la probabilité du mouvement correspondant. Ceci donne une procédure récursive. Pour gagner en efficacité, il est essentiel d'utiliser l'option remember, et il est utile de mettre 0 dès que la distance à l'origine est supérieure à n. 

>	proba:=proc(i,j,n,c) option remember;  if abs(i)+abs(j)>n then 0  elif n=0 then 1.  else      proba(i-1,j,n-1,c)*(1/4+c)+proba(i+1,j,n-1,c)*(1/4-c)      +proba(i,j+1,n-1,c)*1/4+proba(i,j-1,n-1,c)*1/4  fi end:

 

Aussi, on note l'initialisation avec 1. pour n=0 et non 1, ce qui force les calculs à être tous en flottants. 

2. Premières estimations numériques 

>	Digits:=50:

 

>

N:=100:

 

>	for i to N do e[i]:=add(proba(0,0,j,1/10),j=1..2*i) od;

 

Une estimation de la probabilité de repassage à l'origine pour c=1/10 est donc : 

>	e[100]/(1+e[100]);

 

3. Accélération de convergence 

On applique la méthode de Wynn: 

>	for i to N+1 do epsilon[-1,i]:=0 od:

 

>	for i to N do epsilon[0,i]:=e[i] od:

 

>	for i to N-1 do for j to N-i do epsilon[i,j]:=epsilon[i-2,j+1]+1/(epsilon[i-1,j+1]-epsilon[i-1,j]) od; if type(i,even) then print(i,epsilon[i,N-i]/(1+epsilon[i,N-i])) fi od;

 

>	evalf(0.41274995071541979678446489712869325432505616462200,20);

 

Ces décimales semblent acquises. 

4. Valeurs exactes et récurrence 

La même procédure, mais en initialisant à 1 et non 1. pour calculer avec des rationnels. 

>	proba:=proc(i,j,n,c) option remember;  if abs(i)+abs(j)>n then 0  elif n=0 then 1  else      expand(proba(i-1,j,n-1,c)*(1/4+c)+proba(i+1,j,n-1,c)*(1/4-c)      +proba(i,j+1,n-1,c)*1/4+proba(i,j-1,n-1,c)*1/4)  fi end:

 

Il vaut mieux utiliser une version récente de gfun 

>	libname:="/Users/salvy/lib/maple/gfun/lib",libname:gfun:-version();

 

>	seq(add(proba(0,0,2*j,1/10),j=1..i),i=0..20);

 

>	gfun:-listtorec([%],u(k));

 

>	rec1:=%[1]:

 

5. Cinquante décimales de la probabilité 

>	ee:=gfun:-rectoproc(rec1,u(k)):

 

>	for i in 1000,2000,3000 do tmp:=ee(i); print(i,evalf(tmp/(1+tmp))) od:

 

Ces 50 décimales semblent acquises. 

6. Une récurrence pour les polynômes 

La procédure de la question 4 fonctionne aussi avec c symbolique : 

>	proba(0,0,4,c);

 

>	L:=[seq(add(proba(0,0,j,c),j=1..2*i),i=0..20)]:

 

>	gfun:-listtorec(L,u(k));

 

>	rec:=%[1];

 

Vérification : 

>	evalb(eval(rec,c=1/10)=rec1);

 

7. La dérive pour des nombres de pas fixés 

>	vals:=gfun:-rectoproc(rec,u(k),params=[c],evalfun=expand);

 

On sait que c est compris entre -1/4 et 1/4. Par symétrie, on peut se contenter de la racine positive : 

>	for i in 5,6,7 do 2^i,fsolve(vals(2^i,c)-1,c,0..1/4) od;

 

8. Les dérivées 

>	gfun:-listtorec(map(diff,L,c),u(k));

 

>	valdiff:=gfun:-rectoproc(op(1,%),u(k),params=[c]);

 

Vérification : 

>	seq(expand(valdiff(i,c)-diff(L[i+1],c)),i=1..5);

 

9. Une résolution par itération de Newton 

L'idée est de ne pas utiliser les polynômes développés, qui sont très gros, mais d'exploiter les procédures pour calculer leurs valeurs à des c fixés. 

Le test d'arrêt que l'on se donne sur la méthode de Newton porte sur la distance entre deux itérations consécutives. 

>	newtn:=proc(n) local v, w; w:=.1; v:=0; Digits:=Digits+5;    while abs(w-v)>Float(1,5-Digits) do v:=w; w:=v-(vals(n,v)-1)/valdiff(n,v) od; evalf(w,Digits-5) end;

 

Vérification : 

>	newtn(20);

 

>	fsolve(L[21]-1,c,0..1/4);

 

10. Les valeurs et leur limite empirique 

>	for i in 1000,2000,4000 do i,newtn(i) od;

 

>	evalf(%[2],30);

 

11. Une récurrence satisfaite par les  

>	U:=binomial(2*n,2*k)*binomial(2*k,k)*pE^k*pO^k*binomial(2*n-2*k,n-k)*pN^(n-k)*pS^(n-k);

 

>	U:=subs(pE=pEO^2/pO,pN=pNS^2/pS,U);

 

>	SumTools[Hypergeometric][Zeilberger](U,n,k,Sn);

 

La suite sommée étant nulle pour k<0 et k>n, le membre droit s'annule et la récurrence est donnée par le premier élément de cette réponse. 

>	opq:=%[1];

 

>	recq:=add(coeff(opq,Sn,i)*q(n+i),i=0..degree(opq,Sn));

 

12. Preuve de la récurrence de la question 6. 

>	collect(subs(n=k+1,pNS=1/4,pEO=(1/16-c^2)^(1/2),eval(recq,q=proc(n) u(n)-u(n-1) end)),u,expand);

 

>	op(remove(type,rec,`=`));

 

>

%%-%;

 

13. Équation différentielle satisfaite par la série génératrice 

>	deqQ:=gfun:-rectodiffeq({recq,seq(q(i)=add(eval(U,n=i),k=0..i),i=0..1)},q(n),Q(z));

 

14. Sa solution générale 

>	dsolve(op(1,deqQ),Q(z));

 

15. Sous forme hypergéométrique 

>	res:=convert(op(2,%),hypergeom);

 

16. Les conditions initiales 

On commence par isoler les arguments des séries hypergéométriques : 

>	indets(res,function);

 

>

[op(%)];

 

>	map2(op,3,%);

 

>	map(limit,%,z=0);

 

La seconde est donc divergente en z=0, alors que son facteur dans res a une limite finie. On en déduit _C2=0.  

>	res:=subs(_C2=0,res);

 

Il ne reste plus qu'a identifier _C1 : 

>	eval(res,z=0);

 

>	res:=subs(_C1=-1/I,res);

 

17. La convergence de cette série 

>	indets(res,function);

 

>	op([1,3],%);

 

Cette fraction est clairement croissante en z pour 0. Il reste à voir si la valeur à laquelle elle atteint 1 est plus grande que 1 

>	solve(%=1,z);

 

>	factor(%);

 

Les produits pEO et pNS atteignent leur maximum lorsque pE=pO=1/2 et pN=pS=1/2, où ils valent 1/2. En toutes les autres valeurs, cette fraction est donc supérieure à 1, et donc l'argument de l'hypergéométrique reste dans son disque de convergence pour  

18. La formule finale 

On part de l'expression générale 

>

res;

 

et il ne reste plus qu'à remplacer les probabilités et z par leurs valeurs 

>	result:=subs(z=1,pEO=sqrt(1/16-c^2),pNS=1/4,res);

 

>	simplify(result);

 

>	convert(%,hypergeom);

 

>	final:=1-1/%;

 

Cette formule donne finalement la probabilité de repassage par l'origine en fonction de la dérive. (On pourrait aussi utiliser l'expression en termes de fonction elliptique K ci-dessus pour l'évaluer par la moyenne arithmético-géométrique comme vu en TP). 

19. 500 décimales des valeurs cherchées  

>	Digits:=500:

 

>	evalf(subs(c=1/10,final));

 

ceci confirme notre résultat à la question 3. 

>	fsolve(final=1/2,c,0..1/4);

 

ceci confirme notre résultat à la question 10. 

Temps total pour cette session : 

>

time();

 

>