# Leçon 1

Solution de systèmes linéaires 
Résolution d'équation algébrique par 
les méthodes de dichotomie 
          et du point fixe 

## Systèmes linéaires

### Valeurs et vecteurs propre

Décompose une matrice A comme
A = Q-1 L Q
où L est la matrice diagonalisée 

Ne traite que le cas des valeurs propres simples
(ce n'est pas une fonction très stable)

printrat:false$
Decompose(A):=block([x,U:ident(length(A)),p,sol,VP,Q,Qm],
  p:determinant(A-x*U),
  sol:solve(p,x),
  if listp(sol) then sol:map(lambda([%u],ev(x,%u)),sol)
                else sol:allroots(p),
  VP:map(lambda([lam],block([X,eq,s,xp,rn,i,res],
      X:transpose(matrix(makelist(concat('v,i),i,1,length(A)))),
      eq:(A-lam*U).X,
      s:solve(flatten(args(eq)),flatten(args(X)))[1],
      xp:ev(X,s),
      map(lambda([%u],diff(xp,%u)),%rnum_list)[1])),sol),
  Q:apply(matrix,map(lambda([%u],flatten(args(%u))),VP)),print(Q),
  if(length(Q)=length(A)) 
  then block([],Qm:invert(Q),[Qm,apply(diag_matrix,sol),Q])
  else ["il y a des multiplicités dans les valeurs propres : ",sol] 
)$

Decompose(A):=block([valps,vecps],
 [valp,vecp]:eigenvectors(A))$
A:matrix([2,1],[1,2]);
Decompose(A);

A:matrix([1,2],[2,1]);
Decompose(A);

A:matrix([1,2,3],[3,2,1],[4,5,6]);
Decompose(A);

A:matrix([2,1],[1,2])$
AA:matrix([1,0,0,0],[-1,4,1,-2],[2,1,2,-1],[1,2,1,0]);
eigenvectors(A);
[vals,vecs]:eigenvectors(AA);
for i thru length (vals[1]) do disp (val[i] = vals[1][i],
  mult[i] = vals[2][i], vec[i] = vecs[i]);

### Décomposition de Cholesky

Cholesky_decompose(m):=block
([l,dim,n,p,q,res,piv],
 dim:length(m),
 l:zeromatrix(dim,dim),
 for n:1 thru dim do
 (piv:m[n,n],
  for q:1 thru n-1 do piv:piv-l[n,q]^2,
  piv:sqrt(piv),
  l[n,n]:piv,
  for p:n+1 thru dim do 
  (res:m[p,n],
   for q:1 thru n-1 do res:res-l[p,q]*l[n,q],
   l[p,n]:res/piv
  )),l)$

A:matrix([10,7,8,7],[7,5,6,5],[8,6,10,9],[7,5,9,10]);
L:Cholesky_decompose(A);
L . transpose(L);

## Équations algébriques

### Dichotomie

code brut

f(x):=x^3-0.3333$ /*f(x):=cos(3*%pi*x)$ */ 
a:0 $ b:1 $ e:1.0e-7 $
wxplot2d(f(x),[x,a,b]);
/* algorithme */
fa:f(a) $ fb:f(b) $
if fa=0 then  (x_tilde:a) 
elseif fb=0 then (x_tilde:b)
elseif fa*fb > 0 then (x_tilde:concat("echec de l'algorithme. [a,b]=[",a,",",b,"]"))
else
(while abs(b-a)>e do
 (m:(a+b)/2,fm:f(m),
  if fm=0 then (a:m,b:m) 
  elseif fa*fm > 0 then (a:m,fa:fm)
  else (b:m,fb:fm)),
 x_tilde:float((a+b)/2)
);

Sous forme de fonction

ResolutionParDichotomie(Expression,Variable,BorneInf,BorneSup,Precision):=
 block([a:BorneInf,b:BorneSup,e:Precision,x_tilde,m,fm,xx:ev(Variable)],
  block([fa:ev(Expression,ev(Variable)=a),fb:ev(Expression,ev(Variable)=b)],
   if fa=0 then  x_tilde:a 
   elseif fb=0 then x_tilde:b
   elseif fa*fb > 0 then res:concat("echec de l'algorithme : [a,b]=[",a,",",b,"]")
   else
   (while abs(b-a)>e do
    (m:(a+b)/2,fm:ev(Expression,ev(Variable)=m),
    if fm=0 then (a:m,b:m) 
     elseif fa*fm > 0 then (a:m,fa:fm)
     else (b:m,fb:fm)),
     x_tilde:(a+b)/2.0,
     res:[x_tilde,float(ev(Expression,ev(Variable)=x_tilde))]),
   res))$
RPD:ResolutionParDichotomie$

ff:cos(3*%pi*x)$
plot2d(ff,[x,0,1])$
RPD(ff,x,0,1,1.0e-3);

### Point fixe

Fonction simple :
En entrée : 
x  : une variable
f  : une expression en x
x0 : la valeur initiale de la variable
eps: la valeur du critère d'arrêt |f-x| < eps
nblim : le nombre maximum d'itération avant l'arrêt

En sortie : [x_infini,f_evalué_pour_x_infini,nombre_d_iterations_effectuees]

point_fixe(f,x,x0,eps,nblim):=
block([fx0,n,p:0],
fx0:subst(x0,x,f),
for n:0 while (n < nblim  and abs(fx0-x0) > eps) do
(p:p+1,x0:fx0,fx0:float(subst(x0,x,f))),
[x0,fx0,p])$

ff:x^3-1/3;
plot2d(ff,[x,-2,2])$
point_fixe(x^3-1/3,x,0,1.0e-10,100);

Point fixe échoue

ff:2*x-1;
plot2d(ff,[x,-2,2])$
point_fixe(2*x-1,x,0,1.0e-10,100);

Recherche de zéro par la méthode de point fixe

zero_point_fixe(f,x,x0,eps,nblim):=point_fixe(x-f,x,x0,eps,nblim);

zero_point_fixe(x^3-1/3,x,0,1.0e-10,100);