diff --git a/.gitignore b/.gitignore index bee8a64..54c3833 100644 --- a/.gitignore +++ b/.gitignore @@ -1 +1,3 @@ __pycache__ +/Spline/spline +*.cmi diff --git a/Spline/Makefile b/Spline/Makefile new file mode 100644 index 0000000..aac3d33 --- /dev/null +++ b/Spline/Makefile @@ -0,0 +1,8 @@ +%.cmx: %.ml + ocamlopt -linkall -c $< + +spline: spline.cmx + ocamlopt -o spline -linkall str.cmxa spline.cmx + +clean: + rm -rf spline *.cmx *.cmi *.o diff --git a/Spline/generate.sh b/Spline/generate.sh new file mode 100644 index 0000000..0b060d6 --- /dev/null +++ b/Spline/generate.sh @@ -0,0 +1,6 @@ +for i in {0..1000} +do + echo $i + ./spline + mv plot_out.png plot_$i.png +done diff --git a/Spline/spline.ml b/Spline/spline.ml new file mode 100644 index 0000000..902ae1b --- /dev/null +++ b/Spline/spline.ml @@ -0,0 +1,226 @@ +(** CRÉATION DES TYPES **) + +(* Type représentant un polynôme de degré trois à coefficients réels *) +type poly3 = Poly3 of (float*float*float*float) | NotAPoly3;; + +(* Type représentant un polynôme de degré un, à coéfficients réels *) +type poly1 = Poly1 of (float*float) | NotAPoly1;; + +(* Type représentant un polynôme de degré trois à coefficients des polynomes de degré 1 *) +(* Par exemple, on a P(X,Y) = (a+bY)X^3 + (c+dY)X^2 + (e+fY)X^1 + (g+hY) *) +type poly3poly1 = Poly3Poly1 of (poly1*poly1*poly1*poly1) | NotAPolyPoly;; + +(* Représente une spline, sans les coordonées de la subdivisision «portant» les polynomes de degré trois *) +type spline = poly3 array;; + + +(* Librairie pour l'affichage *) +open Owl;; + +let x0 = Mat.uniform 8 8;; + +(* On définit les opérations +,-,.(proportion réelle) sur les polynômes de degré 1*) +let (++) = fun x y -> match (x,y) with + | (Poly1(a,b),Poly1(c,d)) -> Poly1(a+.c,b+.d) + | _ -> NotAPoly1;; + +let (--) = fun x y -> match (x,y) with + | (Poly1(a,b),Poly1(c,d)) -> Poly1(a-.c,b-.d) + | _ -> NotAPoly1;; + +let ( ** ) = fun x y -> match y with + | Poly1(a,b) -> Poly1(x*.a,x*.b) + | _ -> NotAPoly1;; + +(* Cette fonction renvoie la fonction polynomiale associée à un polynome de degré trois (pour l'affichage) *) +let ffromPoly3 p = match p with + | Poly3(a,b,c,d) -> fun x -> a*.x*.x*.x +. b*.x*.x +. c*.x +. d + | NotAPoly3 -> fun x -> nan;; +(* Cette fonction prend en argument une spline (tableau de polynômes de degré trois) et les +coordonées de la subdivision la portant et renvoie la fonction réelle correspondant *) +let ffromPoly3arr arr xs x= let n = Array.length xs in + let rec aux i = + if i=n-1 then 0. (* Hors des champs *) + else + if xs.(i)<=x && x<=xs.(i+1) (* On suppose la continuité *) + then ffromPoly3 arr.(i) x + else aux (i+1) + in aux 0;; + + +(* GÉNÉRATION DE LA SPLINE *) +(* L'idée est ici de considérer des polynômes d'ordre trois «presque déterminés». Ils correspondent au +type poly3poly1 défini plus haut. On peux considérer que quand une fonction renvoie un tel objet, il +renvoie une représentation de l'ensemble des objets respectant les contraintes spécifiées, cet ensemble +étant {P(X)(y) pour y réel}*) + +(* Récupère le polynôme «presque déterminé» correspondant à la première section. Il y a trois contraintes pour +quatres inconnues (les coefficients du polynôme). On renvoie donc le polynome presque déterminé P(X)(Y) +tel que si l'on choisit le coefficient en X^3 du polynôme pour être égal à "a", alors le polynôme respectant +les contraintes sera P(X)(a) *) +let getFirst (x:float) (y:float) (alpha:float) (beta:float) = + let a = Poly1(1.,0.) in + let b = Poly1((-3.) *. x , 0.) in + let c = Poly1( (2.*.x*.x) +. (2.*.x*.y) -. (y*.y) , (beta -. alpha) /. (y -. x)) in + let d = Poly1( (y*.y*.x)-.(2.*.x*.x*.y) , alpha -. (((beta -. alpha) /. (y -. x)) *. x)) in + Poly3Poly1(a,b,c,d);; + +(* On fait de la même manière que ci-dessus, avec cette fois ci quatre contraintes et quatre inconnues. +Par contre, on calcule avec les coefficients qui sont des polynômes de degré un, car ils sont encore «indéterminés»*) +let getMiddle (x:float) (y:float) (alpha:poly1) (beta:poly1) (gamma:poly1) (delta:poly1) = + let a = ((1./.(y-.x))/.(y-.x)) ** (((1./.(y-.x))**(delta--alpha)) -- beta -- (((y-.x)/.2.)**gamma)) in + let b = (1./.2.)**gamma -- ((3.*.x)**a) in + let c = (beta -- (x**gamma)) ++ ((3.*.x*.x)**a) in + let d = ((alpha -- ((x*.x*.x)**a)) -- ((x*.x)**b)) -- x**c in + Poly3Poly1(a,b,c,d);; + +(* Cette fonction renvoie l'unique valeur de y tel que le polynôme pp3(X)(y) vérifie la contrainte à gauche (dérivée seconde nulle) *) +let getEndLambda pp3 y = match pp3 with + | Poly3Poly1(Poly1(aL,aC),Poly1(bL,bC),Poly1(cL,cC),Poly1(dL,dC)) -> -1. *. (bC +. 3.*.aC*.y)/.(3.*.aL*.y +. bL) + | _ -> nan;; + +(* Cette fonction est (P,x) -> P(x)(X) *) +let applyAndFriendsPoly1Poly3 pp3 x = match pp3 with + | Poly3Poly1(a,b,c,d) -> ( + (x*.x*.x)**a ++ (x*.x)**b ++ x**c ++ d, + (3.*.x*.x)**a ++ (2.*.x)**b ++ c, + (6.*.x)**a ++ 2.**b ) + | NotAPolyPoly -> (NotAPoly1,NotAPoly1,NotAPoly1);; + +(* Cette fonction est (l,P) -> P(l) *) +let solvePoly1 lambda poly = match poly with + | Poly1(a,b) -> a*.lambda+.b + | NotAPoly1 -> nan;; + +(* Cette fonction est (l,P) -> P(X)(l) *) +let solvePoly1Poly3 lambda pp3 = match pp3 with + | Poly3Poly1(a,b,c,d) -> Poly3(solvePoly1 lambda a,solvePoly1 lambda b,solvePoly1 lambda c, solvePoly1 lambda d) + | NotAPolyPoly -> NotAPoly3;; + +(* Cette fonction applique la fonction ci-dessus à un tableau *) +let solvePolyArray (lambda:float) (arr:poly3poly1 array) = + let n = Array.length arr in + let out = Array.make n NotAPoly3 in + for i=0 to (n-1) do + out.(i) <- solvePoly1Poly3 lambda arr.(i) + done; + out;; + +(* ENFIN on fusionne toutes les fonctions faites précedemment *) +let interpolate xs ys = + let n = Array.length xs in + let pp3 = Array.make (n-1) NotAPolyPoly in + pp3.(0) <- getFirst xs.(0) xs.(1) ys.(0) ys.(1); (* Le premier polynomes «indéterminé» est renvoyé par getFirst *) + for i=1 to (n-2) do (* Ensuite, pour chacun des segments de la subdivision *) + let alpha,beta,gamma = applyAndFriendsPoly1Poly3 pp3.(i-1) xs.(i) and delta = Poly1(0.,ys.(i+1)) in (* On explicite nos trois contraintes *) + (* On pourrait mettre alpha = Poly1(0.,ys.(i)) *) + pp3.(i) <- getMiddle xs.(i) xs.(i+1) alpha beta gamma delta (* Puis on récupère le polynôme indéterminé imposé sur ce segment*) + done; + let lambda = getEndLambda pp3.(n-2) xs.(n-1) in (* Enfin, on récupère la valeur de lambda=le coefficient en X^3 de notre premier polynôme *) + solvePolyArray lambda pp3;;(* Avant de l'appliquer à tous *) + + +(** En dessous, des tests d'utilisation des fonctions **) +let xs = [|0.;1.;2.;3.;4.;5.;6.;7.;8.;9.|];; +let ys = [|-7.;14.;-19.;15.;-2.;-11.;-9.;-16.;10.;3.|];; +interpolate xs ys;; +interpolate [|-1.;0.;1.;2.|] [|1.;-1.;0.;2.|];; +getFirst (-1.) 0. 1. (-1.);; + +getMiddle 1. 2. (Poly1(0.,1.)) (Poly1(2.,1.)) (Poly1(6.,0.)) (Poly1(0.,0.));; + +getEndLambda (Poly3Poly1(Poly1(-5., -2.),Poly1(18., 6.), Poly1(-19., -5.), Poly1(6., 2.))) 2.;; + +let ppp1 x = match x with + | Poly1(a,b) -> + print_string (string_of_float a); + print_string "X + "; + print_endline (string_of_float b) + | NotAPoly1 -> ();; + +open Owl;; + +let linterpol xs ys x = let n = Array.length xs in + let rec aux i = + if i=n then 0. (* Hors des champs *) + else + if xs.(i)<=x && x<=xs.(i+1) (* On suppose la continuité *) + then ys.(i) +. (ys.(i+1) -. ys.(i)) *. (x -. xs.(i)) /. (xs.(i+1)-.xs.(i)) + else aux (i+1) + in aux 0;; + + +let getLI xs i = + let n = Array.length xs in + let out = Array.make n 0. in + out.(0)<-1.; + for j=0 to n-1 do + if i<>j then + let xixj = 1. /. (xs.(i) -. xs.(j)) in + for k=n-1 downto 1 do + out.(k) <- xixj *. (out.(k-1) -. xs.(j)*.out.(k)) + done; + out.(0) <- (-1.) *.xixj *. xs.(j) *. out.(0) + done; + out;; + +let polyCLToFirst p lambda q = + let n = Array.length p in + for i=0 to (n-1) do + p.(i) <- p.(i) +. (lambda *. q.(i)) + done;; + +let lagrange xs ys = + let n = Array.length xs in + let out = Array.make n 0. in + for i=0 to n-1 do + polyCLToFirst out ys.(i) (getLI xs i) + done; + out;; + +let polyfun poly x = let out = ref 0. and xx = ref 1. in + for i=0 to ((Array.length poly)-1) do + out := !out +. (poly.(i) *. !xx); + xx := !xx *. x + done; + !out;; + +lagrange [|0.;1.;2.|] [|0.;1.;3.|];; + +(* Génération de points aléatoires à interpoler *) +Random.self_init ();; +let n = 10;; +let xs = Array.make n 0.;; +for i=1 to (n-1) do + xs.(i) <- (xs.(i-1) +. (Random.float 10.) +. 1.) +done;; +let ys = Array.make n 0.;; +for i=0 to (n-1) do + ys.(i) <- (-10.) +. (Random.float 20.) +done;; + +let toPlot1 = ffromPoly3arr (interpolate xs ys) xs;; +let toPlot2 = polyfun (lagrange xs ys);; +let toPlot3 = linterpol xs ys;; +let minX = xs.(0) and maxX = xs.((Array.length xs)-1);; + + +let pplot h f = + let x = Mat.linspace minX maxX 1200 in + let y = Mat.map f x in + Plot.plot ~h ~spec:[] x y;; + +let h = Plot.create "plot_out.png" in +Plot.set_background_color h 12 12 12; +Plot.set_foreground_color h 0 0 255; +Plot.set_pen_size h 2.; +Plot.plot_fun ~h toPlot2 minX maxX; +Plot.set_foreground_color h 0 255 240; +Plot.set_pen_size h 4.; +Plot.plot_fun ~h toPlot3 minX maxX; +Plot.set_foreground_color h 42 255 0; +Plot.set_pen_size h 2.; +Plot.plot_fun ~h toPlot1 minX maxX; +Plot.set_foreground_color h 200 200 200; +Plot.set_yrange h (-20.) 20.; +Plot.output h;;