Mathématiques expérimentales Semestre 3


Home Projet Documents Sources

Plus d'images et le code

Images

Voici un exemple du changement de l'image à partir du changement des données initiales dans ℂ :

Donnée initiale: 2i, 3i
Donnée initiale: 2i, 4i
Donnée initiale: 2i, 5i

Autres exemples:

Donnée initiale: 1i, 3i
Donnée initiale: 1i, 4i
Exemple de l'orbite en ℝ + τℝ, données initiales: (2, 1) ; (3, 2)


Le code

Code pour les calculs:

Canopy est nécessaire pour le fonctionnement du code.

		

""" Author: Wolff Vincent, Geimer Arno, Karst Philippe """ import numpy as np import numpy.linalg as linal import math from scipy.spatial import Voronoi import scipy.spatial.distance as scp import matplotlib.pyplot as plt from split_complex_eml import Split_complex as spl import mayavi.mlab as mlab def dist(p1,p2): """Computes the euclidean distance between p1 and p2.""" return np.sqrt((p1[0]-p2[0])**2+(p1[1]-p2[1])**2) def create_words(minLength,maxLength): """Function used by complex_j_points and r_tau_r_j_points """ end_list=[] start_list=["0","1","2","3"] for word in start_list: l=0 aux_list=[] aux_list.append(word) for i in range(1,maxLength,1): l1=len(aux_list) for j in range(l,l1,1): if aux_list[j][len(aux_list[j])-1]=="0": aux_list.append(aux_list[j]+"0") aux_list.append(aux_list[j]+"2") aux_list.append(aux_list[j]+"3") elif aux_list[j][len(aux_list[j])-1]=="1": aux_list.append(aux_list[j]+"1") aux_list.append(aux_list[j]+"2") aux_list.append(aux_list[j]+"3") elif aux_list[j][len(aux_list[j])-1]=="2": aux_list.append(aux_list[j]+"0") aux_list.append(aux_list[j]+"1") aux_list.append(aux_list[j]+"2") else: aux_list.append(aux_list[j]+"0") aux_list.append(aux_list[j]+"1") aux_list.append(aux_list[j]+"3") l=l1 for line in aux_list: if len(line)>=minLength: end_list.append(line) return end_list def listOfMatrices(s,t): """Returns the list of matrices A,B such that Tr(A*B*A^(-1)*B^(-1))=-2.""" liste=[] c=1 a=-0.5*s+0.5*t+0.5*np.sqrt(2*(2*c-1)*s*t+s*s+t*t+4) l=np.sqrt((2*np.sqrt(c*s*t+1))/(c*s*t)+2/(c*s*t)+1) liste.append(np.matrix([[a,s*t],[c,a+s-t]])) liste.append(linal.inv(liste[0])) liste.append(np.matrix([[l,0],[0,1/l]])) liste.append(linal.inv(liste[2])) return liste def wordToMatrix(word,list_of_matrices): """Function used by complex_j_points andr_tau_r_j_points to convert a word to a matrix according a given list of matrices.""" result=np.matrix([[1,0],[0,1]]) for i in range(0,len(word),1): result=result*list_of_matrices[int(word[i])] return result def moebius_transformation_complex_j(A,z,t): """Returns for a given complex number z and a matrix A=[[a,b],[c,d]] with complex coefficients and t>0 a float the number (a*(z+t*j)+b)/(c*(z+t*j)+d) """ liste=[] q=abs(A[1,0]*z+A[1,1])**2+(abs(A[1,0])**2)*t*t liste.append(((A[0,0]*z+A[0,1])*(A[1,0]*z+A[1,1]).conjugate()+A[0,0]*A[1,0].conjugate()*t**2)/q) liste.append(t/q) return liste def complex_j_points(lmin=1,n=3): """Given the complex numbers s and t, the maximal length of words lmax, a complex number z and a float t>0 the function will return a list of points corresponding to the orbit of z+t*j. To represent the points of the list use draw_scatter. """ maxT=10**(-n) print("Enter s and t complex numbers.") s=complex(input()) t=complex(input()) print("Entrer the max length of the words.") lmax=int(input()) list_matrices=listOfMatrices(s,t) print("Enter the real part of z.") re=float(input()) print("Enter the imaginary part of z.") im=float(input()) z=re+1j*im print("Enter t>0.") t2=float(input()) points=[] list_of_words=create_words(lmin,lmax) for line in list_of_words: M=wordToMatrix(line,list_matrices) z1=moebius_transformation_complex_j(M,z,t2) if z1[1]<=maxT: points.append([z1[0].real,z1[0].imag]) return points def draw_scatter(points): """Given a list of points the scatter plot will be drawn.""" x=[] y=[] for p in points: x.append(p[0]) y.append(p[1]) plt.scatter(x,y,s=1) plt.show() def draw_line(points): """Draws the curve through the points of the given list.""" x=[] y=[] for p in points: x.append(p[0]) y.append(p[1]) plt.plot(x,y,'-') plt.show() def r_tau_r_j_points(lmin=7,n=5,sort=True): """Given the real numbers s1, t1, s2 and t2, the maximal length of word lmax, a split complex number z and a float t>0 the function will return a list of points corresponding to the orbit of z+t*j. The list of points will be sorted with respect to the x-coordinate if the value of the additional argument sort is True. Default value is True. The additional argument lmin is the minimal length of the words used to compute the points. To represent the points of the list use draw_scatter. To represent the curve through these points use draw_line. """ maxT=10**(-n) print("Enter the real numbers s1 and t1.") s=float(input()) t=float(input()) list_of_matrices=listOfMatrices(s,t) print("Enter the real numbers s2 and t2.") s=float(input()) t=float(input()) list_of_matrices2=listOfMatrices(s,t) print("Entrer the max length of the words.") lmax=int(input()) print("Enter the real part of z.") re=float(input()) print("Enter the split_imaginary part of z.") im=float(input()) z=spl(re,im) print("Enter t>0.") t2=float(input()) z1=[z,t2] print("Enter the endpoint xmax of the interval [-xmax,xmax].") xmax=float(input()) list_of_words=create_words(lmin,lmax) points=[] for line in list_of_words: M1=wordToMatrix(line,list_of_matrices) M2=wordToMatrix(line,list_of_matrices2) M=spl.sl_2_to_spl(M1,M2) z2=spl.moebius_transformation_r_tau_r_j(M,z1,None) if z2!=None: if z2[1]<=maxT and abs(z2[0].real)<=xmax and abs(z2[0].spl_im)<=xmax: points.append([z2[0].real,z2[0].spl_im]) if sort: points.sort(key=lambda pair:pair[0]) return points def eliminate_points(points,slope=1): """Given a list of points sorted with respect to the x-coordinate the function will remove all the points such that the slope of the line through 2 consecutive points is in absolute value bigger than the parameter. Default value of the parameter slope is 1. """ indices=[] for i in range(1,len(points),1): if abs((points[i][1]-points[i-1][1])/(points[i][0]-points[i-1][0]))>slope: indices.append(i) l=[] for i in range(0,len(points),1): if not(i in indices): l.append(points[i]) return l def max_Circles(points,rmin=10**(-3)): """Given a list of points a list of maximal circles will be returned. The list is of the form [c_x,c_y,c_r] where c_x is the x-coordinate and c_y is the y-coordinate of the center and c_r is the radius. The last element of the list is the value xmax used to draw the circles. The additional argument rmin is the smallest possible radius. To get all circles set rmin to 0. To draw the circles use draw_circles. """ vor=Voronoi(points) points.sort(key=lambda pair:pair[0]) if abs(points[0][0])>abs(points[len(points)-1][0]): xmax=abs(points[0][0]) else: xmax=abs(points[len(points)-1][0]) liste=[x for x in vor.vertices if dist(x,[0,0])<=xmax] end_list=[] for p in liste: Y=scp.cdist([p],vor.points) end_list.append([p[0],p[1],Y.min()]) end_list.append(xmax) return end_list def draw_circles(circles): """Draws the circles of the list.""" ax = plt.gca() ax.cla() for i in range(0,len(circles)-1,1): ax.add_artist(plt.Circle((circles[i][0],circles[i][1]),circles[i][2],fill=False)) ax.set_xlim(-circles[len(circles)-1],circles[len(circles)-1]) ax.set_ylim(-circles[len(circles)-1],circles[len(circles)-1]) plt.show() def center_radius_of_hyperbola(x1,y1,x2,y2,x3,y3): """Computes the center of the hyperbola passing through the points (x1,y1), (x2,y2) and (x3,y3) if it exists. Returns None if such a hyperbola does not exist. """ q=((x2 - x3)*y1 - (x1 - x3)*y2 + x1*y3 - x2*y3) if q==0: return None else: den_x=-1/2*((y1 - y3)*(y2**2) - (x1**2)*y3 + (x2**2)*y3 + (y1**2)*y3 - (x2**2 - x3**2 + y3**2)*y1 + (x1**2 - x3**2 - y1**2 + y3**2)*y2) den_y=1/2*((x1 - x3)*(x2**2) + (x1**2)*x3 + (x2 - x3)*(y1**2) - (x1 - x3)*(y2**2) - (x3**2 - y3**2)*x1 - (x1**2 - x3**2 + y3**2)*x2) den_x=den_x/q den_y=den_y/q r=(y1-den_y)**2-(x1-den_x)**2 return [den_x,den_y,r] def list_hyperbolas(points,hyperbolas): """Function used by max_Hyperbolas.""" list_center=[] for hyp in hyperbolas: center=center_radius_of_hyperbola(points[hyp[0]][0],points[hyp[0]][1],points[hyp[1]][0], points[hyp[1]][1],points[hyp[2]][0],points[hyp[2]][1]) if center!=None: list_center.append(center) return list_center def eq_hyperbola(cx,cy,x,y): """Returns the number (y-cy)**2-(x-cx)**2.""" return (y-cy)**2-(x-cx)**2 def rearange(points,hyp,i): """Function used by max_Hyperbolas.""" hyp1=[hyp[0],hyp[1],i] hyp2=[hyp[1],hyp[2],i] hyp3=[hyp[0],hyp[2],i] hyp1=sort_hyperbola(points,hyp1) hyp2=sort_hyperbola(points,hyp2) hyp3=sort_hyperbola(points,hyp3) liste=[] if not(inside(points,hyp1,hyp[2])): liste.append(hyp1) if not(inside(points,hyp2,hyp[0])): liste.append(hyp2) if not(inside(points,hyp3,hyp[1])): liste.append(hyp3) return liste def inside(points,hyp,i): """Function used by max_Hyperbolas.""" center=center_radius_of_hyperbola(points[hyp[0]][0],points[hyp[0]][1],points[hyp[1]][0], points[hyp[1]][1],points[hyp[2]][0],points[hyp[2]][1]) if center==None: return True return eq_hyperbola(center[0],center[1],points[i][0],points[i][1])-center[2]>0 def check(points,hyperbolas,hyp): """Function used by max_Hyperbolas.""" b2=True for k in range(0,len(points),1): if b2: if not(k in hyp): if inside(points,hyp,k): t=rearange(points,hyp,k) if len(t)!=0: hyperbolas.append(t[0]) if len(t)>=2: hyperbolas.append(t[1]) if len(t)==3: hyperbolas.append(t[2]) if hyp in hyperbolas: hyperbolas.remove(hyp) check(points,hyperbolas,t[0]) if len(t)>=2: check(points,hyperbolas,t[1]) if len(t)==3: check(points,hyperbolas,t[2]) b2=False else: if hyp in hyperbolas: hyperbolas.remove(hyp) def sort_hyperbola(points,hyp): """Sorts the list hyp such that hyp[0] has the lowest x-coordinate.""" list_pts=[] for i in range(0,3,1): list_pts.append([points[hyp[i]][0],hyp[0]]) list_pts.sort(key=lambda pair:pair[0]) return [list_pts[0][1],list_pts[1][1],list_pts[2][1]] def max_Hyperbolas(points): """Given a list of points a list of maximal hyperbolas will be returned. The list is of the form [c_x,c_y,r]. The equation of the corresponding hyperbola is (y-c_y)^2-(x-c_x)^2=r. """ hyperbolas=[] hyperbolas.append([0,1,2]) for i in range(3,len(points),1): b=True for j in range(0,len(hyperbolas),1): if j=2: hyperbolas.append(t[1]) if len(t)==3: hyperbolas.append(t[2]) hyperbolas.remove(hyperbolas[j]) b=False else: hyperbolas.remove(hyperbolas[j]) if b: hyperbolas.append([i-2,i-1,i]) check(points,hyperbolas,hyperbolas[len(hyperbolas)-1]) for hyp in hyperbolas: check(points,hyperbolas,hyp) t2=[hyperbolas[0]] hyperbolas.remove(hyperbolas[0]) for hyp in hyperbolas: if hyp in t2 or [hyp[0],hyp[2],hyp[1]] in t2 or [hyp[1],hyp[0],hyp[2]] in t2 or [hyp[1],hyp[2],hyp[0]] in t2 or [hyp[2],hyp[0],hyp[1]] in t2 or [hyp[2],hyp[1],hyp[0]] in t2: pass else: t2.append(hyp) return list_hyperbolas(points,t2) def plot_sphere(R,a,b,half=True): """Plots a sphere with radius R and center (a,b). If the additional argument half is set to False, then spheres will be drawn, else half-spheres.""" [phi,theta] = np.mgrid[0:2*np.pi:15j,0:np.pi:15j] x = R*np.cos(phi)*np.sin(theta) + a y = R*np.sin(phi)*np.sin(theta) + b if half: z = abs(R*np.cos(theta)) else: z=R*np.cos(theta) return mlab.mesh(x, y, z) def draw_half_spheres(circles): """Draws spheres for a given list of circles and s,t and maxLength represents the parameters used to create the original set of points. """ for i in range(0,len(circles)-1,1): plot_sphere(circles[i][2],circles[i][0],circles[i][1]) def hyperbola_points(center,n=2,prec=100): """Returns two lists of points on the hyperbola of given center and radius.""" l1=[] l2=[] cx=center[0] cy=center[1] radius=center[2] if radius>0: interval=np.linspace(-n,n,prec) for number in interval: l1.append([math.sqrt(radius)*number+cx,math.sqrt(radius)*math.sqrt(1+number*number)+cy]) l2.append([math.sqrt(radius)*number+cx,-math.sqrt(radius)*math.sqrt(1+number*number)+cy]) return [l1,l2] else: return None def draw_hyperbola(center,n=2,prec=100,show=True): """Given the center and radius of an hyperbola, the hyperbola will be drawn.""" hyp=hyperbola_points(center,n=n,prec=prec) if hyp!=None: l1=hyp[0] l2=hyp[1] x=[p[0] for p in l1] y=[p[1] for p in l1] plt.plot(x,y,'-b') x=[p[0] for p in l2] y=[p[1] for p in l2] plt.plot(x,y,'-b') if show: plt.show() def draw_listOfHyperbolas(list_Hyperbola): """Given a list of hyperbolas, the corresponding picture will be drawn.""" for hyp in list_Hyperbola: draw_hyperbola(hyp,show=False) plt.show() def draw_hyperboloid(center,half=False): """Draws an hyperboloid of given center and radius.""" cx=center[0] cy=center[1] if center[2]>0: r=np.sqrt(center[2]) [s,t] = np.mgrid[-2:2:15j,0:2*np.pi:15j] x = r*s + cx y = r*np.sqrt(1+s**2)*np.sin(t) + cy if half: z = abs(r*np.sqrt(1+s**2)*np.cos(t)) else: z = r*np.sqrt(1+s**2)*np.cos(t) return mlab.mesh(x, y, z) def draw_listOfHyperboloid(list_of_hyperbolas): """Draws a list of hyperboloids.""" for hyp in list_of_hyperbolas: draw_hyperboloid(hyp,half=True) def plot_3dline(points): """Plots a 3d line passing through the given points.""" x=[p[0] for p in points] y=[p[1] for p in points] z=[0 for i in range(0,len(x))] mlab.plot3d(x,y,z)

Code pour les opérations (addition, multiplication,...) en ℜ+τℜ:

		

""" Author Wolff Vincent, Geimer Arno, Karst Philippe Class provides implementation of most of the standard operations for split complex numbers except /. """ from math import sqrt class Split_complex(object): def __init__(self,real,split_imag): self.real=real self.spl_im=split_imag def __str__(self): return(str(self.real)+"+tau*"+str(self.spl_im)) __repr__=__str__ def __add__(self,other): real=self.real+other.real im=self.spl_im+other.spl_im return Split_complex(real,im) def __sub__(self,other): return Split_complex(self.real-other.real,self.spl_im-other.spl_im) def __mul__(self,other): real=self.real*other.real+self.spl_im*other.spl_im im=self.spl_im*other.real+self.real*other.spl_im return Split_complex(real,im) def __rmul__(self,other): return Split_complex(self.real*other,self.spl_im*other) def __pow__(self,other): res=Split_complex(1,0) for i in range(1,other+1,1): res=res*self return res def inverse(self): """For a split complex number z returns z divided by z.norm_2() if defined.""" return Split_complex.div_relle(self.conjuguate(),self.norm_2()) def norm_2(self): """For a split complex number x+tau*y returns x^2-y^2.""" return self.real*self.real-self.spl_im*self.spl_im def norm(self): """For a split complex number x+tau*y returns sqrt(x^2-y^2) if x^2-y^2>=0. """ return sqrt(self.norm_2()) def conjuguate(self): """For a split complex number x+tau*y returns x-tau*y.""" return Split_complex(self.real,-self.spl_im) @staticmethod def mult_relle(z,q): """Function used by div_relle. For multiplication of split complex number z with float a write a*z.""" return Split_complex(z.real*q,z.spl_im*q) @staticmethod def div_relle(z,q): """ For a split complex number z=x+tau*y returns z2=x/q+tau*y/q if q!=0.""" return Split_complex.mult_relle(z,1/q) @staticmethod def moebius_transf(A,z): """ For a split complex number z and a matrix A=[[a,b],[c,d]] with split complex coefficients returns (a*z+b)/(c*z+d) if defined.""" a=A[0][0] b=A[0][1] c=A[1][0] d=A[1][1] z1=a*z+b z2=c*z+d res=z1*z2.conjuguate() q=z2.norm_2() return Split_complex.div_relle(res,q) @staticmethod def moebius_transformation_r_tau_r_j(A,z,*args): """ For a number z=z1+t*j where t>0, z1 split complex and a matrix A=[[a,b],[c,d]] with split complex coefficients returns (a*z+b)/(c*z+d) if defined. If an additional argument is given (e.g. an split complex number) this will be returned if A*z is not defined.""" z1=z[0] t=z[1] a=A[0][0] b=A[0][1] c=A[1][0] d=A[1][1] z2=a*z1+b z3=c*z1+d q=z3.norm_2() c=c.conjuguate() z3=z3.conjuguate() q=q+c.norm_2()*t*t res=z2*z3+t*t*a*c if q==0 and len(args)!=0: return args[0] else: return [Split_complex.div_relle(res,q),t/q] @staticmethod def sl_2_to_spl(A,B): """ For A,B square matrices of dimension 2 the function returns the matrix A*omega1+B*omega2 where omega1=(1+tau)/2 and omega2=(1-tau)/2.""" C=[] omega1=Split_complex(0.5,0.5) omega2=Split_complex(0.5,-0.5) C.append([A[0,0]*omega1+B[0,0]*omega2,A[0,1]*omega1+B[0,1]*omega2]) C.append([A[1,0]*omega1+B[1,0]*omega2,A[1,1]*omega1+B[1,1]*omega2]) return C @staticmethod def birapport(z1,z2,z3,z4): """ Returns ((z1-z3)*(z2-z4))/((z1-z4)*(z2-z3)) if defined.""" z5=(z1-z3)*(z2-z4) z6=(z1-z4)*(z2-z4) q=z6.norm_2() z6=z6.conjuguate() res=Split_complex.div_relle(z5*z6,q) return res.spl_im