4.5 Algèbre
fichier source
let make_vect = Array.make;;
let length_vect = Array.length;;
(*******************************************************************)
let make_matrix n p =
let m = make_vect n [||] in
for i = 0 to (n-1) do
m.(i) <- make_vect p 0.;
done;
m;;
(*******************************************************************)
let nb_lines m = length_vect m;;
(*******************************************************************)
let nb_columns m = length_vect m.(0);;
(*******************************************************************)
let fill_line m i f =
let p = nb_columns m in
for j = 0 to (p-1) do
m.(i).(j) <- f j;
done;;
(*******************************************************************)
let fill_matrix m f =
let n = nb_lines m in
for i = 0 to (n - 1) do
fill_line m i (function j -> f (i,j));
done;;
(*******************************************************************)
let identity n =
let m = make_matrix n n in
fill_matrix m (function i,j -> if (i = j) then 1. else 0.);
m;;
(*******************************************************************)
let copy_matrix m =
let c = make_matrix (nb_lines m) (nb_columns m) in
fill_matrix c (function i,j -> m.(i).(j));
c;;
(*******************************************************************)
let integers n =
let m = make_matrix 1 n in
fill_matrix m (function _,j -> float_of_int (j+1));
m;;
(*******************************************************************)
let str_of_matrix m =
let str = ref "" in
let n = nb_lines m in
let p = nb_columns m in
for i = 0 to (n - 1) do
for j = 0 to (p - 1) do
str := !str^" "^string_of_float m.(i).(j);
done;
str := !str^"\n"
done;
!str;;
(*******************************************************************)
let print_matrix m =
print_string (str_of_matrix m);;
(*******************************************************************)
let matrix_of_list l =
let m = make_matrix 1 (List.length l) in
let rec fill i l = match l with
[] -> ()
| h::t -> m.(0).(i) <- h ; fill (i+1) t in
fill 0 l;
m;;
(*******************************************************************)
let list_of_vector v =
let rec list_of_vector_acc acc = function
-1 -> acc
| j -> list_of_vector_acc (v.(j)::acc) (j-1) in
list_of_vector_acc [] (length_vect v - 1);;
(*******************************************************************)
let list_of_matrix m =
let rec list_of_matrix_acc acc m = function
-1 -> acc
| i -> list_of_matrix_acc ((list_of_vector m.(i))::acc) m (i-1) in
list_of_matrix_acc [] m (nb_lines m - 1);;
(*******************************************************************)
let transpose m =
let n = nb_lines m in
let p = nb_columns m in
let t = make_matrix p n in
fill_matrix t (function i,j -> m.(j).(i));
t;;
(*******************************************************************)
let rotation m =
let n = nb_lines m in
let p = nb_columns m in
let t = make_matrix n p in
fill_matrix t (function i,j -> m.(n - i - 1).(p - j - 1));
t;;
(*******************************************************************)
exception Dimensions_mismatch;;
let add_matrix a b =
let la = nb_lines a in
let lb = nb_lines b in
let ca = nb_columns a in
let cb = nb_columns b in
if (la = lb && ca = cb) then
let c = make_matrix la ca in
fill_matrix c (function i,j -> a.(i).(j) +. b.(i).(j));
c;
else
raise Dimensions_mismatch;;
(*******************************************************************)
let trace m =
let n = nb_lines m in
let p = nb_columns m in
if (n <> p) then
raise Dimensions_mismatch
else
let somme = ref 0. in
for i = 0 to (n - 1) do
somme := !somme +. m.(i).(i);
done;
!somme;;
(*******************************************************************)
let mult_matrix_const m k =
let n = nb_lines m in
let p = nb_lines m in
let c = make_matrix n p in
fill_matrix c (function (i,j) -> k *. m.(i).(j));
c;;
(*******************************************************************)
let mult_matrix a b =
let la = nb_lines a in
let lb = nb_lines b in
let ca = nb_columns a in
let cb = nb_columns b in
if (ca = lb) then
let c = make_matrix la cb in
let fill_function (i,j) =
(
let total = ref 0. in
for k = 0 to (ca - 1) do
total := !total +. a.(i).(k) *. b.(k).(j);
done;
!total;
) in
fill_matrix c fill_function;
c;
else
raise Dimensions_mismatch;;
(*******************************************************************)
let rec exp_matrix m n =
if (n = 0) then
identity (nb_columns m)
else
mult_matrix m (exp_matrix m (n - 1));;
(*******************************************************************)
let add_linear_combination m (lg, wg) (ld, wd) =
fill_line m lg (function j -> m.(lg).(j) *. wg +. m.(ld).(j) *. wd);;
(*******************************************************************)
let pivot m inv (a, b) =
let n = nb_lines m in
let valeur_pivot = m.(a).(b) in
for i = (a+1) to (n-1) do
let r = -.m.(i).(b) in
add_linear_combination inv (i, valeur_pivot) (a, r);
add_linear_combination m (i, valeur_pivot) (a, r);
done;;
(*******************************************************************)
exception Singular_matrix;;
let find_pivot m i =
let k = ref i in
let p = nb_lines m in
while(!k < p && m.(!k).(i) = 0.) do
k := !k + 1;
done;
if !k = p then
raise Singular_matrix;
!k;;
(*******************************************************************)
let swap_items m (l1, c1) (l2, c2) =
let tmp = m.(l1).(c1) in
m.(l1).(c1) <- m.(l2).(c2);
m.(l2).(c2) <- tmp;;
(*******************************************************************)
let swap_lines m i k =
for j = 0 to (nb_columns m - 1) do
swap_items m (i, j) (k, j);
done;;
(*******************************************************************)
let iteration_pivot m inv i =
let k = find_pivot m i in
if (i <> k) then
(swap_lines m i k;
swap_lines inv i k;
);
pivot m inv (i, i);;
(*******************************************************************)
let reduce_diagonal m inv =
let n = nb_lines m in
for i = 0 to (n-1) do
let coeff = m.(i).(i) in
fill_line m i (function j -> m.(i).(j) /. coeff);
fill_line inv i (function j -> inv.(i).(j) /. coeff);
done;;
(*******************************************************************)
let rec apply_function f x = function
0 -> x
| n -> apply_function f (f x) (n-1);;
(*******************************************************************)
let invert_matrix m =
let n = nb_lines m in
let reduce_and_rotate (a, b) =
(
for k = 0 to (n - 2) do
iteration_pivot a b k;
done;
(rotation a, rotation b)
)
in
let (m, inv) = apply_function reduce_and_rotate (copy_matrix m, identity n) 2 in
reduce_diagonal m inv ;
inv;;
(*******************************************************************)
let solve_system a y =
mult_matrix (invert_matrix a) y;;