Ajout du super programme interpolant une liste de points avec une application polynomiale de degré 3 par morceaux et C^2

.gitignore

@@ -1 +1,3 @@
 __pycache__
+/Spline/spline
+*.cmi

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

Spline/generate.sh

@@ -0,0 +1,6 @@
+for i in {0..1000}
+do
+	echo $i
+	./spline
+	mv plot_out.png plot_$i.png
+done

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;;