1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41
| open Bigarray
let schrodinger1D c d n =
let f argument = 0.
and potentiel arg = arg
and h = (d -. c) /. (float_of_int (n+1))
and x = Array.create (n+2) 0.
and matriceRigidite = Array2.create float64 fortran_layout n 2
and vecteurColonne = Array.make n 0.
and z = Array.make n 0.
and u = Array.make (n-1) 0.
and d = Array.make n 0.
and v = Array.make n 0. in
for i=0 to (n+1) do
x.(i) <- c +. float_of_int i *. h;
done;
for i=1 to (n-1) do
matriceRigidite.{i,1} <- (2. /.h +. h/.3. *. (potentiel ((x.(i-1) +. x.(i)) /.2. ) +. potentiel (x.(i)) +. potentiel ((x.(i) +. x.(i+1)) /. 2.) ));
matriceRigidite.{i,2} <- 1. /.h +. h/.6. *. (potentiel (x.(i) +. x.(i+1) ) /. 2. );
done;
matriceRigidite.{n,1}<- 2. /.h +. h/.3. *. (potentiel ((x.(n-1) +. x.(n)) /.2. ) +. potentiel (x.(n)) +. potentiel ((x.(n) +. x.(n+1)) /. 2.) );
for i=0 to n-1 do
vecteurColonne.(i) <- h /. 3. *. (f (( x.(i) +. x.(i+1) ) /. 2.) +. f (x.(i+1)) +. f (( x.(i+1) +. x.(i+2) ) /. 2.) );
done;
d.(0) <- matriceRigidite.{1,1};
z.(0) <- (vecteurColonne.(0) /. d.(0));
u.(0) <- (matriceRigidite.{1,2} /. d.(0));
for i=1 to (n-2) do
d.(i) <- (matriceRigidite.{i+1,1} -. matriceRigidite.{i,2} *. u.(i-1));
u.(i) <- (matriceRigidite.{i+1,2} /. d.(i));
z.(i) <- ((vecteurColonne.(i) -. matriceRigidite.{i,2} *. z.(i-1)) /. (d.(i)));
done;
d.(n-1) <- matriceRigidite.{n,1} -. matriceRigidite.{n-1,2} *. u.(n-2);
z.(n-1) <- ((vecteurColonne.(n-1) -. matriceRigidite.{n-1,2} *. z.(n-1)) /. (d.(n-1)));
v.(n-1) <- z.(n-1);
for i=(n-2) downto 0 do
v.(i) <- z.(i) -. u.(i) *. v.(i+1);
done;
v;; |
Partager