Commit bb3677b7 by Benjamin Gallois

Initial

parent 159484e8
%{*******************************************************************************
* @file fonctions.m
* @author Benjamin GALLOIS
* @date xx/05/2016
* @version final
* @resume Computes the distance between two points.
******************************************************************************}%
function [d] = distance(x,y,theta0)
d = sqrt((x*cos(y)-cos(theta0))^2+(x*sin(y)-sin(theta0))^2);
end
'''/*******************************************************************************
* @file divergence.py
* @author Benjamin GALLOIS
* @date xx/05/2016
* @version final
* @resume Plot the angle of divergence according to the parameter of the
problem G
*******************************************************************************/'''
from fonctions import *
#T=True arrete la simulation quand il y a convergence
def simulation(G,N,T):
#initialisation
T=0.1
R0=1
V0=G/T
E=np.zeros((N,2))
angle=[]
B=ini(1)
E[0,0]=B[0]
E[0,1]=B[1]
#evolution
for i in range(N-1):
E=evolution(T,E,V0,R0,N)
b=placement(1,E)
E[i+1,0]=1
E[i+1,1]=b
angle.append(np.abs(E[i+1,1]-E[i,1]))
meanangle=0
if angle[i]>(np.pi):
angle[i]=np.pi*2-angle[i]
print 'En cours itération',i,'/',N-1
angle[i]=(angle[i]*180)/np.pi
if i>1 and T=True:
if round(angle[i-1])==round(angle[i]) and round(angle[i-2])==round(angle[i]) and angle[i]!=0:
meanangle=angle[i]
break
else:
continue
if i>1 and T=False:
meanangle=angle[i]
return E[:,0],E[:,1],meanangle
G=np.linspace(0.01,0.4,50)
G0=[]
angle=[]
T=True
for i in G:
A=simulation(i,50,T)
angle.append(A[0])
plt.figure()
plt.subplot()
plt.plot(G,angle,'x')
plt.xlabel('G')
plt.ylabel('Angle de divergence')
plt.subplot()
plt.plot(A[0]*np.cos(A[1]),A[0]*np.sin(A[1]),'o')
theta=np.linspace(0,2*np.pi,100)
plt.plot(np.cos(theta),np.sin(theta))
plt.xlabel('x')
plt.ylabel('y')
plt.savefig('Divergence', format='png')
plt.show()
%{*******************************************************************************
* @file energy.m
* @author Benjamin GALLOIS
* @date xx/05/2016
* @version final
* @resume Computes the energy at the distance d.
******************************************************************************}%
function [E] = energy(d)
if d ~= 1
E = exp(-d/0.01);
else
E = 0;
end
end
%{*******************************************************************************
* @file fonction.m
* @author Benjamin GALLOIS
* @date xx/05/2016
* @version final
* @resume Computes the function of G accordinf to time.
******************************************************************************}%
function G=fonction(N,n)
G1 = 0.8*exp(-15*(n/N));
G2 = 0.5*log((N-n+1)/(N-n));
if G1>G2
G = G1;
end
if G2>=G1
G=G2;
end
end
%{*******************************************************************************
* @file fonction2.m
* @author Benjamin GALLOIS
* @date xx/05/2016
* @version final
* @resume Computes the function of G according to time.
******************************************************************************}%
function G=fonction2(N,n)
step = 1.5/N;
G=1.5-step*n;
end
'''/*******************************************************************************
* @file fonctions.py
* @author Benjamin GALLOIS
* @date xx/05/2016
* @version final
* @resume Contains all functions needed to run the simulations.
*******************************************************************************/'''
import numpy as np
import matplotlib.pyplot as plt
import scipy as sc
from matplotlib import animation
#initialisation, place un premier point aléatoirement sur le cercle de rayon R0
def ini(R0):
i=np.random.randint(0,360)
A=np.zeros(2)
A[0]=R0
A[1]=(i*np.pi)/180
return A
#energie de repulsion d'une seule particule en A(ra,thetaa) au point B(rb,thetab)
def energy(ra,thetaa,rb,thetab):
xa=ra*np.cos(thetaa)
xb=rb*np.cos(thetab)
ya=ra*np.sin(thetaa)
yb=rb*np.sin(thetab)
d=np.sqrt((xa-xb)**2+(ya-yb)**2)
return np.exp(-d/0.1)
def placement(R0,E):
theta=np.linspace(0,2*np.pi,2000)
En=np.zeros(len(E[:,0]))
ener=[]
def somme(E,En,R0,theta):
for j in range(len(E[:,0])):
if E[j,0]!=0:
En[j]=energy(E[j,0],E[j,1],R0,theta)
else:
En[j]=0
return np.sum(En)
for i in range(len(theta)):
ener.append(somme(E,En,R0,theta[i]))
A=np.min(ener)
B=np.argwhere(ener==A)
if len(B)>1:
minu=B[np.random.randint(0,len(B))]
return theta[minu]
else:
return theta[B]
#fait se déplacer chaque point avec une vitesse radiale v(r)
def evolution(T,E,V0,R0,N):
for i in range(len(E[:,0])):
if E[i,0]!=0:
E[i,0]=E[i,0]+((V0*E[i,0]*T)/R0)
return E
%{*******************************************************************************
* @file main.m
* @author Benjamin GALLOIS
* @date xx/05/2016
* @version final
* @resume Simulation of the growing.
******************************************************************************}%
clear()
N = 800;
data = zeros(N,2);
angle = zeros(N-1,1);
circle = linspace(0,2*pi,3600);
G = zeros(1,N);
for k = 1:N
if k == 1
data(1,1) = 1;
data(1,2) = 0;
G(k,1) = fonction(N,k);
data(1,1) = data(1,1)+G(k,1)*data(1,1);
end
if k >= 2
G(k,1) = fonction(N,k);
for i = 1:N;
data(i,1) = data(i,1)+G(k,1)*data(i,1);
end
ener = zeros(1,3600);
for i = 1:3600;
E = zeros(1,N);
for j = 1:N
d = distance(data(j,1),data(j,2),circle(i));
E(1,j) = energy(d);
if d == 1
break
end
end
ener(1,i) = sum(E(1,:));
end
[A,B] = min(ener(1,:));
data(k,1) = 1;
data(k,2) = circle(B);
angle(k-1,1) = abs(data(k,2)-data(k-1,2));
if angle(k-1,1)>pi;
angle(k-1,1)=2*pi-angle(k-1,1);
end
subplot(2,2,1);
scatter(data(:,1).*cos(data(:,2)),data(:,1).*sin(data(:,2)),10,'filled');
axis square
axis off
axis([-25,25,-25,25])
subplot(2,2,2);
p1=plot((angle(1:k-1)*180)/pi,'-bo','MarkerSize',5,'MarkerFaceColor','b');
axis([0,N,0,200])
ylabel('Divergence angle')
xlabel('Time')
moy=num2str((angle(k-1,1)*180)/pi);
text(k-1+10,(angle(k-1,1)*180)/pi, moy);
subplot(2,2,3)
p2 = plot(G(1:k),'-bo','MarkerSize',5,'MarkerFaceColor','b');
axis([0,N,0,0.8])
ylabel('$G=\frac{T V_0}{R_0}$','Interpreter','LaTex')
xlabel('Time')
text2=num2str(k);
pause(0.00001);
print(text2,'-dpng','-r600')
end
end
%{*******************************************************************************
* @file main2.m
* @author Benjamin GALLOIS
* @date xx/05/2016
* @version final
* @resume Simulation of the growing.
******************************************************************************}%
clear()
N = 200;
data = zeros(N,2);
angle = zeros(N-1,1);
circle = linspace(0,2*pi,3600);
G = zeros(1,N);
for k = 1:N
if k == 1
data(1,1) = 1;
data(1,2) = 0;
H = fonction2(N,k);
G(k,1) = H ;
data(1,1) = data(1,1)+H*data(1,1);
end
if k >= 2
H = fonction2(N,k);
G(1,k) = H ;
for i = 1:N;
data(i,1) = data(i,1)+H*data(i,1);
end
ener = zeros(1,3600);
for i = 1:3600;
E = zeros(1,N);
for j = 1:N
d = distance(data(j,1),data(j,2),circle(i));
E(j) = energy(d);
if d == 1
break
end
end
ener(i) = sum(E);
end
[A,B] = min(ener);
data(k,1) = 1;
data(k,2) = circle(B);
angle(k-1,1) = abs(data(k,2)-data(k-1,2));
if angle(k-1,1)>pi;
angle(k-1,1)=2*pi-angle(k-1,1);
end
subplot(2,2,1);
plot(data(:,1).*cos(data(:,2)),data(:,1).*sin(data(:,2)),'bo');
subplot(2,2,2);
p1=plot((angle(1:k-1)*180)/pi);
p1.Marker = 'o'
axis([0,N,0,200])
moy=num2str((angle(k-1,1)*180)/pi)
text(k-1,(angle(k-1,1)*180)/pi, moy)
ylabel('Divergence angle')
xlabel('Time')
subplot(2,2,3)
p2 = plot(G(1:k))
p2.Marker = 'o'
ylabel('$G=\frac{T V_0}{R_0}$','Interpreter','LaTex')
xlabel('Time')
axis([0,200,0,1.5])
pause(0.001)
end
end
'''/*******************************************************************************
* @file patern.py
* @author Benjamin GALLOIS
* @date xx/05/2016
* @version final
* @resume Plot the bifurcation diagram according to the parameter of the
simulation G.
*******************************************************************************/'''
from fonctions import *
def simulationaveccondition(N,theta):
E=np.zeros((N,2))
angle=[]
B=ini(1)
E[0,0]=B[0]
E[0,1]=B[1]
R0=1.
V0=1.
G=[]
G0=0.001
pas=1./N
for i in range(N-1):
if i<(N/4):
G0=G0
T=G0/V0
E=evolution(T,E,V0,R0,N)
E[i+1,0]=1
E[i+1,1]=E[i,1]+theta
angle.append(np.abs(E[i+1,1]-E[i,1]))
angle[i]=(angle[i]*180)/np.pi
G.append(G0)
else:
G0=G0+pas
T=G0/V0
E=evolution(T,E,V0,R0,N)
b=placement(1,E)
E[i+1,0]=1
E[i+1,1]=b
angle.append(np.abs(E[i+1,1]-E[i,1]))
if angle[i]>(np.pi):
angle[i]=np.pi*2-angle[i]
angle[i]=(angle[i]*180)/np.pi
G.append(G0)
if i>(N/4+2):
if abs(angle[i]-angle[i-1])<5.:
continue
if abs(angle[i]-angle[i-1])>5.:
break
else:
continue
if len(G)>(N/3):
T=True
else :
T=False
return G,angle,T
def simulationlibre(N,H):
E=np.zeros((N,2))
angle=[]
B=ini(1)
E[0,0]=B[0]
E[0,1]=B[1]
R0=1.
V0=2
G=[]
pas=H/N
G0=H+pas
for i in range(N-1):
G0=G0-pas
print G0
T=G0/V0
E=evolution(T,E,V0,R0,N)
b=placement(1,E)
E[i+1,0]=1
E[i+1,1]=b
angle.append(np.abs(E[i+1,1]-E[i,1]))
if angle[i]>(np.pi):
angle[i]=np.pi*2-angle[i]
angle[i]=(angle[i]*180)/np.pi
G.append(G0)
return G,angle
for i in np.linspace(60,120,60):
print i
A=simulationaveccondition(300,(i*np.pi)/180.)
if A[2]==True:
plt.plot(A[0],A[1],'r.')
B=simulationlibre(300,1.2)
plt.plot(B[0],B[1],'r.')
plt.xlabel('G')
plt.ylabel('Angle de divergence')
plt.axis([0, 1.2, 60, 180])
plt.title("Diagramme de bifurcation de l'angle de divergence.")
plt.savefig('Bifurcation2.png')
plt.show()
'''/*******************************************************************************
* @file simulation.py
* @author Benjamin GALLOIS
* @date xx/05/2016
* @version final
* @resume Simulates the growing and plot the divergence angle according to
the parameter of the problem G.
*******************************************************************************/'''
from fonctions import *
import time
#prend (le rayon, la vitesse, le nombre de points et la periode entre apparition de chaque point), renvoit pour chaque point les coordonnées polaires (r, theta, G, l'angle de divergence avec le voisin precedant, la valeur de G, l'angle de divergence quand il y a convergence)
def simulation(R0,V0,N,T,H):
#initialisation
G=(V0*T/R0)
E=np.zeros((N,2))
angle=[]
B=ini(1)
E[0,0]=B[0]
E[0,1]=B[1]
#evolution
for i in range(N-1):
E=evolution(T,E,V0,R0,N)
b=placement(1,E)
E[i+1,0]=1
E[i+1,1]=b
angle.append(np.abs(E[i+1,1]-E[i,1]))
meanangle=0
if angle[i]>(np.pi):
angle[i]=np.pi*2-angle[i]
print 'En cours itération',i,'/',N-1
angle[i]=(angle[i]*180)/np.pi
if i>1:
if round(angle[i-1])==round(angle[i]) and round(angle[i-2])==round(angle[i]) and angle[i]!=0:
meanangle=angle[i]
if H==True:
break
else:
continue
return E[:,0],E[:,1],angle,G,meanangle
N=input('Nombre de points ?')
G=input('G ?')
T=0.1
V0=G/T
start_time = time.time()
A=simulation(1,V0,N,T,False)
print("--- %s seconds ---" % (time.time() - start_time))
plt.figure()
plt.plot(A[0]*np.cos(A[1]),A[0]*np.sin(A[1]),'o')
theta=np.linspace(0,2*np.pi,100)
plt.plot(np.cos(theta),np.sin(theta))
plt.xlabel('x')
plt.ylabel('y')
plt.figure()
plt.plot(range(len(A[2])),A[2],'o')
plt.xlabel('(i+1)-i')
plt.ylabel('Angle de divergence')
plt.show()
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment