> restart:

MAO_TP - E1

Évaluation du 10/01/2007

>

On charge quelques bibliothèques

Vérifier les numéros de version des bibliothèmes disponibles

> with(pldx): # il faut la version 29

> with(simul): # il faut la version 21

Warning, the name ever_loaded has been redefined

> with(LinearAlgebra):

> xima:= proc(pl); op(1,pl): convert(%,list): map(op,%): map2(op,2,%); (max@op)(%); end:

> kernelopts(ASSERT=true); with(plots, display);

> xmacro(); macro(ptplo=plots[pointplot]);

>

Homographies A

Construction du sujet

1.7, -2.3; (z-%[1])/(z-%[2]); phi:= unapply(%,z);

-0.8*phi(z)=phi(Z); (normal@expand)(solve(%, Z)); f:= unapply(%, z);

>

Points fixes

> f := unapply((-7.*z+351.9)/(90.*z+47.), z);

> alpha, beta:= solve(f(z)=z);

> ka, kb:= D(f)(alpha), D(f)(beta);

Escargot

> plot(f(t), t=-10..10, view=[-10..10, -10..10],
color=black, linestyle=16, discont=true),
plot(t, t=-10..10): display(%); pl1:=%:

>

> x0:= -2.7; li:= seq( (f@@k)(x0), k=0..40);

> abs(li[-1]-li[-2]);

>

> pl4:= plots[pointplot]([li[1]$3, seq(X$4, X=li[2..-2]), li[-1]$3 ], style=line, color=magenta):

repérage de coordonnées sur le graphe

> display(pl1, pl4, view=[-8..8, -8..8], tickmarks=[[-2.7,-10,4],[-10,4]], labels=[``,``]);

>

>

Homographies B

Construction du sujet

1.7, -2.4; (z-%[1])/(z-%[2]); phi:= unapply(%,z);

-0.8*phi(z)=phi(Z); (normal@expand)(solve(%, Z)); f:= unapply(%, z);

>

Points fixes

> f:= unapply( (-5.5*z+183.6)/(45.*z+26.) , z);

> alpha, beta:= solve(f(z)=z);

> ka, kb:= D(f)(alpha), D(f)(beta);

Escargot

> plot(f(t), t=-10..10, view=[-10..10, -10..10],
color=black, linestyle=16, discont=true),
plot(t, t=-10..10): display(%); pl1:=%:

>

> x0:= -2.7; li:= seq( (f@@k)(x0), k=0..40);

> abs(li[-1]-li[-2]);

>

> pl4:= plots[pointplot]([li[1]$3, seq(X$4, X=li[2..-2]), li[-1]$3 ], style=line, color=magenta):

repérage de coordonnées sur le graphe

> display(pl1, pl4, view=[-8..8, -8..8], tickmarks=[[-2.7,-10,4],[-10,4]], labels=[``,``]);

>

>

Matrices

> f(z): ``(numer(%))/denom(%); qq:= [numer, denom]( f(z)); ma:= Matrix(2,2,(j,k)->coeff(qq[j],z,2-k));

> mu, mp:= Eigenvectors(ma);

> Column(mp,1); %/%[2];

> alpha,beta; ma.<alpha,1>; -%/82;

>

>

Création des fichiers

fich:= "/home/douillet/docs/Ensait/maotp/maotp_ds17a/dat_maotp_ds17a.txt";
N:=nextprime(410); p:= 11; fu:= unapply(23-5*x+x^2/2,x); plot(fu(x), x=0..10);
evp:= evalf(Pi): g:= proc(z); fu(z/33.)+6*cos(2*evp*z/p)+4*ra(); end;
fclose(fich); hdl:= fopen(fich, WRITE);
for i to N do fprintf(hdl,"%a\n", g(i)); od:
fclose(fich);

fich:= "/home/douillet/docs/Ensait/maotp/maotp_ds17b/dat_maotp_ds17b.txt";
N:=nextprime(410); p:= 11; fu:= unapply(23-5*x+x^2/2,x); plot(fu(x), x=0..10);
evp:= evalf(Pi): g:= proc(z); fu(z/33.)+6*cos(2*evp*z/p)+4*ra(); end;
fclose(fich); hdl:= fopen(fich, WRITE);
for i to N do fprintf(hdl,"%a\n", g(i)); od:
fclose(fich);

>

Série Temporelle A

Lecture du fichier

Déclarer le nom de la copie locale du fichier

> fich:= "F:/docs/Ensait/maotp/maotp_ds17a/dat_maotp_ds17a.txt";
fich:= "/home/douillet/docs/Ensait/maotp/maotp_ds17/dat_maotp_ds17a.txt";

On lit tous les nombres

> lire:= table(): tmp:= 'tmp':
hdl:= fopen(fich, READ) :
for i to 5000 while tmp <> [ ] do
tmp:= fscanf(hdl,"%a");
lire[i]:= op (tmp);
od:
fclose(fich); n:= i-2 ; i:='i':

On en fait une liste

> li:= [seq(lire[k],k=1..n)]:

>

Paramètres de dispersion

> m:= moy(li);

> s, s2:= sqrt(var(li)), var(li);

> vals:= xab(li); nbr(li);

>

Visualisation

> ptplo([seq([k, li[k]], k=1..n)], style= line, view=[DEFAULT, 0..rhs(vals)], color=magenta); pl_data:= %:

> pl2:= plot([m-s, m, m+s], t=1..n, color=red, linestyle=[8,24,8]) :

> plots[display](pl_data, pl2);

>

>

Transformée de Fourier

nn:=20: ft:= table([]);

ft:= [seq( add(li[k]*exp(2*I*evalf(Pi)*k*j/nn), k=1..nn)/nn, j=0..nn)];

ptplo( [seq([k, abs(ft[k+1])], k=0..nn)], style=line, view=[DEFAULT, 0..20]);

>

Attention : cela dure "un certain temps"

> nn:=n : ft:= table([]);

> pif:= evalf(2*I*Pi/nn): expp:= proc(kj) option remember; exp(pif*kj); end;

> stamp:= time() :
ft:= evalf([seq( add(li[k]*expp(modp(k*j,nn)), k=1..nn)/nn, j=0..nn)]) :
elapsed_time = time() - stamp ;

> ptplo( [seq([k, abs(ft[k+1])], k=0..nn)], style=line, view=[DEFAULT, 0..20]);

> decal:= 2: [seq(abs(ft[k+1]), k=1+decal..n/2)]: member((max@op)(%), %, 'ouca'); ouca:= ouca+decal;

> ici:=ouca-4..ouca+4;
plot([seq([k, abs(ft[k+1])], k=ici)]); pic:= add(k * abs(ft[k+1]), k=ici)/add(abs(ft[k+1]), k=ici);

> n/pic; per:= round(%); lon:= floor(n/per);

> perg := floor(per/2); perd := per-%;

>

Moyenne mobile (trend)

> moymob:= proc(ma::list, per) global perg, perd; local k, kx, mb, tmp;
kx := nops(ma):
for k from 1 to perg-1 do mb[k] := missing od:;
for k from kx-perd+1 to kx do mb[k] := missing od:
tmp := add(ma[k], k=1..per); mb[perg] := tmp/per;
for k from perg+1 to kx-perd do
tmp := tmp-ma[k-perg]+ma[k+perd]: mb[k] := %/per
od: [seq(mb[k],k=1..kx)];
end:

> trend:= moymob(li, per):

> pl_trend:= plot([seq([k, trend[k]],k=1..n)], color=blue):

> plots[display](pl_data, pl_trend);

>

>

Correction des Variations saisonnières

> for u from 1 to per do meu[u]:= moy([seq(li[u+per*j],j=0..lon-1)]); od:
meu_:= moy([seq(meu[u],u=1..per)]);
for u from 1 to per do meu[u]:= meu[u]-meu_: od:

> evalf([seq(meu[m], m=1..per)], 4);
plot( [seq([k, meu[1+irem(k-1,per)]], k=1..30)]);

>

> licor:= [seq(li[k] -meu[1+irem(k-1,per)], k=1..n)] :

> [nbr, moy, var](licor[1..lon*per]); [nbr, moy, var](li[1..lon*per]);

>

> pl_sai_sim:= plot([seq( [k, licor[k] ], k=1..n)], linestyle=16):

> plots[display]( pl_data, pl_sai_sim);

>

Résiduels

> resid := [missing$(perg-1), seq( licor[k] - trend[k], k=perg .. n-perd), missing$perd] : nops(%), n;

> plot([seq( [k, resid[k]], k=1..n)], linestyle=8) ;

> xhisto( subs(missing=NULL, resid), area=1, numbars=15);

>

>

Série Temporelle B

Lecture du fichier

Déclarer le nom de la copie locale du fichier

> fich:= "F:/docs/Ensait/maotp/maotp_ds17a/dat_maotp_ds17a.txt";
fich:= "/home/douillet/docs/Ensait/maotp/maotp_ds17/dat_maotp_ds17b.txt";

On lit tous les nombres

> lire:= table(): tmp:= 'tmp':
hdl:= fopen(fich, READ) :
for i to 5000 while tmp <> [ ] do
tmp:= fscanf(hdl,"%a");
lire[i]:= op (tmp);
od:
fclose(fich); n:= i-2 ; i:='i':

On en fait une liste

> li:= [seq(lire[k],k=1..n)]:

>

Paramètres de dispersion

> m:= moy(li);

> s, s2:= sqrt(var(li)), var(li);

> vals:= xab(li); nbr(li);

>

Visualisation

> ptplo([seq([k, li[k]], k=1..n)], style= line, view=[DEFAULT, 0..rhs(vals)], color=magenta); pl_data:= %:

> pl2:= plot([m-s, m, m+s], t=1..n, color=red, linestyle=[8,24,8]) :

> plots[display](pl_data, pl2);

>

>

Transformée de Fourier

nn:=20: ft:= table([]);

ft:= [seq( add(li[k]*exp(2*I*evalf(Pi)*k*j/nn), k=1..nn)/nn, j=0..nn)];

ptplo( [seq([k, abs(ft[k+1])], k=0..nn)], style=line, view=[DEFAULT, 0..20]);

>

Attention : cela dure "un certain temps"

> nn:=n : ft:= table([]);

> pif:= evalf(2*I*Pi/nn): expp:= proc(kj) option remember; exp(pif*kj); end;

> stamp:= time() :
ft:= evalf([seq( add(li[k]*expp(modp(k*j,nn)), k=1..nn)/nn, j=0..nn)]) :
elapsed_time = time() - stamp ;

> ptplo( [seq([k, abs(ft[k+1])], k=0..nn)], style=line, view=[DEFAULT, 0..20]);

> decal:= 2: [seq(abs(ft[k+1]), k=1+decal..n/2)]: member((max@op)(%), %, 'ouca'); ouca:= ouca+decal;

> ici:=ouca-4..ouca+4;
plot([seq([k, abs(ft[k+1])], k=ici)]); pic:= add(k * abs(ft[k+1]), k=ici)/add(abs(ft[k+1]), k=ici);

> n/pic; per:= round(%); lon:= floor(n/per);

> perg := floor(per/2); perd := per-%;

>

Moyenne mobile (trend)

> moymob:= proc(ma::list, per) global perg, perd; local k, kx, mb, tmp;
kx := nops(ma):
for k from 1 to perg-1 do mb[k] := missing od:;
for k from kx-perd+1 to kx do mb[k] := missing od:
tmp := add(ma[k], k=1..per); mb[perg] := tmp/per;
for k from perg+1 to kx-perd do
tmp := tmp-ma[k-perg]+ma[k+perd]: mb[k] := %/per
od: [seq(mb[k],k=1..kx)];
end:

> trend:= moymob(li, per):

> pl_trend:= plot([seq([k, trend[k]],k=1..n)], color=blue):

> plots[display](pl_data, pl_trend);

>

>

Correction des Variations saisonnières

> for u from 1 to per do meu[u]:= moy([seq(li[u+per*j],j=0..lon-1)]); od:
meu_:= moy([seq(meu[u],u=1..per)]);
for u from 1 to per do meu[u]:= meu[u]-meu_: od:

> evalf([seq(meu[m], m=1..per)], 4);
plot( [seq([k, meu[1+irem(k-1,per)]], k=1..30)]);

>

> licor:= [seq(li[k] -meu[1+irem(k-1,per)], k=1..n)] :

> [nbr, moy, var](licor[1..lon*per]); [nbr, moy, var](li[1..lon*per]);

>

> pl_sai_sim:= plot([seq( [k, licor[k] ], k=1..n)], linestyle=16):

> plots[display]( pl_data, pl_sai_sim);

>

Résiduels

> resid := [missing$(perg-1), seq( licor[k] - trend[k], k=perg .. n-perd), missing$perd] : nops(%), n;

> plot([seq( [k, resid[k]], k=1..n)], linestyle=8) ;

> xhisto( subs(missing=NULL, resid), area=1, numbars=15);

>

>

>