Utilisateur:RM77/Maple
Test de primalité
modifier> restart; > n:=138: > k:=2: > while (n mod k<>0) and (k^2<n) do > k:=k+1: > od: > if k^2>n > then print (`n est premier`); > else print (`n n’est pas premier, il est divisible par`, k); > fi;
Liste des n premiers nombres premiers
modifier> restart; > n:=1000: > a:=2:L:=[2]: > while nops(L)<n do > a:=a+1: > premier:=true: > x:=1: > while premier=true and L[x]^2<=a do > if a mod L[x]=0 then premier:=false:fi: > x:=x+1: > od: > if premier=true then L:=[op(L),a]:fi: > od: > L;
Tableau de Pascal
modifier> restart; > n:=10: > C:=array(1..n+1,1..n+1): > C[1,1]:=1: > for k from 2 to n+1 do > C[k,1]:=1:C[k,k]:=1:od: > for i from 3 to n+1 do > for j from 2 to i-1 do > C[i,j]:=C[i-1,j-1]+C[i-1,j]:od:od: > eval(C);
Attention : éviter des valeurs de n supérieures à 12, problème d'affichage.
Une preuve de l’existence du chaos en maths
modifier> restart; > lambda:=3.9: > f:=x->lambda*x*(1-x): > u:=1/2: > i:=[]: > for k from 0 to 1000 do > u:=evalf(f(u)): > i:=[op(i),[u,u],[u,evalf(f(u))]]: > od: > with(plots): > a:=plot([x,f(x)],x=0..2,color=[black,black]): > b:=pointplot(i,color=red,symbol=cross,connect=true): > display(a,b,view=[0..1.2,0..1.2],scaling=constrained);
Faire varier la valeur de lambda pour voir les variations de résultat.
Valeurs intéressantes :
- 3.9 (chaos)
- 3.5 (période)
- 4 (constante)
Évolution d'une population de lynx et de lapins
modifier> restart; > a:=0.02;b:=0.0001;p:=0.00001;q:=0.02; > f:=(x,t)->x*(1+a-b*t); > g:=(x,t)->t*(1+p*x-q); > u:=1000;v:=100;i:=[];n:=1000; > for k from 0 to n do > tamp1:=u;tamp2:=v; > u:=f(u,v); > v:=g(tamp1,tamp2); > i:=[op(i),[u,v]]; > od: > with(plots): > pointplot(i,color=red,symbol=cross);
Diagramme de bifurcation (ou de Feigenbaum)
modifier> restart; > with(plots): > f:=(x,l)->4*l*x*(1-x); > u0:=0.3;L:=NULL; > for k from 0 to 200 do > l:=evalf(k/200): > u:=u0: > for p from 0 to 99 do > u:=f(u,l):od: > for p from 100 to 300 do > u:=f(u,l): > L:=L,[l,u]: > od:od: > pointplot([L],symbol=point,view=[.6..1,0..1]);
Attracteur de Lorenz
modifier> restart; > s:=10;b:=evalf(8/3);r:=28;k:=0.01; > f:=(x,y,z)->x+k*(-s*x+s*y); > g:=(x,y,z)->y+k*(-x*z+r*x-y); > h:=(x,y,z)->z+k*(x*y-b*z); > x:=0;y:=1;z:=0; > L:=NULL; > for m from 1 to 10000 do > L:=L,[f(x,y,z),g(x,y,z),h(x,y,z)]: > n:=x;o:=y;x:=f(x,y,z);y:=g(n,y,z);z:=h(x,o,z): > od: > L:=[L]: > with(plots): > pointplot3d(L,connect=true);
Attracteur de Hénon
modifier> restart; > a:=1.4;b:=0.3; > f:=(x,t)->t-a*x^2+1; > g:=(x,t)->b*x; > u:=0;v:=0;i:=[];n:=10000; > for k from 0 to n do > tamp1:=u;tamp2:=v; > u:=f(u,v); > v:=g(tamp1,tamp2); > i:=[op(i),[u,v]]; > od: > with(plots): > pointplot(i,color=red,symbol=point);
2 exemples de procédures : la suite de Fibonacci
modifier- Version itérative
> fibo:=proc(n) local A,B,C,k: > A:=0: > B:=1: > for k from 2 to n do > C:=A+B: > A:=B:B:=C: > od: > if n=0 then RETURN(A); > else RETURN(B);fi: > end:
- Version récursive
> fibo2:=proc(n) local F: > if n=0 then F:=0: > else if n=1 then F:=1: > else F:=fibo2(n-1)+fibo2(n-2): > fi:fi: > RETURN(F); > end:
- Attention : cet algorithme, dans l'état actuel, a une complexité trop importante pour pouvoir être efficace. fibo2(40) utilise 2^40 opérations... On peut contourner le problème comme ceci (option remember, qui permet d’éviter de refaire les calculs déjà effectués) :
- Version récursive
> fibo3:=proc(n) option remember: > if n=0 then F:=0: > else if n=1 then F:=1: > else F:=fibo2(n-1)+fibo2(n-2): > fi:fi: > RETURN(F); > end:
Approximation : méthode de Newton
modifier> restart; > approx:=proc(b,f,e) local g,L,k,u; > g:=diff(f,x); > L:=[b,b-(subs(x=b,f)/subs(x=b,g))]; > k:=1; > u:=b-(subs(x=b,f)/subs(x=b,g)); > while evalf(abs(L[k+1]-L[k])) >=e do > u:=evalf(u-subs(x=u,f)/subs(x=u,g)); > L:=[op(L),u];k:=k+1; > od; > RETURN(evalf(L[k],5)); > end:
Utilisation :
> approx(2,(x->x^2-2)(x),10^(-4));
1.4142
Un peu de dessin
modifier- Carré de 20*20 en 2 couleurs
> with(plots): > n:=20: > L1:=[]: > for a from 0 to n do > for b from 0 to a do > L1:=[op(L1),[a,n-b]]:od:od: > A:=pointplot(L1,color=red,axes=none,scaling=constrained): > L2:=[]: > for a from 0 to n do > for b from 0 to a do > L2:=[op(L2),[n-a,b]]:od:od: > B:=pointplot(L2,color=blue,axes=none,scaling=constrained): > display(A,B);
- Cercle rouge au centre d'un carré bleu
> with(plots): > L3:=[]:n:=100:L4:=[]: > for a from 1 to n do > for b from 1 to n do > if (a-(n/2))^2+(b-(n/2))^2<200 then > L3:=[op(L3),[a,b]]: else L4:=[op(L4),[a,b]]:fi: > od:od: > C:=pointplot(L3,color=red):E:=pointplot(L4,color=blue):display(C,E,axes=none,scaling=constrained);
Méthode d'Euler
modifier> restart;
> dsolve({diff(y(x),x)+y(x)=sin(x),y(0)=0});
> assign(%):
> h:=0.5:
> L1:=[0]:
> for k from 2 to 2001 do
> L1:=[op(L1),h*(sin(h*k)-L1[k-1])+L1[k-1]]:od:
> L2:=[]:
> for k from 1 to 2000 do
> L2:=[op(L2),[h*k,L1[k]]]:od:
> with(plots):
> A:=plot(y(x),x=0..20,color=blue):
> B:=pointplot(L2,connect=true,color=red):
> display(A,B);
Penser à faire varier les valeurs de h...
Dichotomie
modifier> restart; > f:=x->x^3+x-1; > a:=0;b:=1;epsilon:=1e-8; > while b-a > epsilon do > c:=(a+b)/2; > if f(a)*f(c)>0 > then a:=c; > else b:=c; > fi; > od: > print (a,b);
Inconvénient : donne une valeur exacte... La solution consiste à placer un point . après une des valeurs des deux bornes de l'intervalle de départ (dans cet exemple, on a une fonction qui nécessite une valeur approchée pour pouvoir appliquer l'algorithme) :
> restart; > f:=x->(exp(-x))-x; > a:=0;b:=1.;epsilon:=1e-8; > while b-a > epsilon do > c:=(a+b)/2; > if f(a)*f(c)>0 > then a:=c; > else b:=c; > fi; > od: > print (a,b);