Utilisateur:RM77/Maple

Test de primalité

> 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

> 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

> 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

> 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

> 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)

> 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]);

Diagramme de Feigenbaum

Attracteur de Lorenz

> 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 Lorenz

Attracteur de Hénon

> 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);

Attracteur de Hénon

2 exemples de procédures : la suite de Fibonacci

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

> 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

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

> 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

> 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);