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 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79
| from random import random
import numpy as np
from math import*
import matplotlib as plt
from scipy.integrate import odeint
import numpy.random as rnd
Xmax=4
Xmin=1
N=5
A=2000
B=0.08
r=0.2
K=120000
V=1
tau=0.5
m=80
X0=rnd.random((2*N))*(Xmax-Xmin)+Xmin
V0=rnd.random(2*N)
Y0=np.concatenate((X0,V0))
Z=[5,2]
temps=np.linspace(0,100,1000)
def vec(a,b):
return b-a
def ps(u,v):
return np.inner(u,v)
def norme(v):
return sqrt(ps(v,v))
def unitaire(v):
return (1/norme(v))*v
def dir (A,B):
return unitaire(vec(A,B))
def g(x):
if x>0:
return 0
else :
return x
def f(i,j):
if i==j :
return 0
else :
h=norme(vec(i,j))
D=h-2*r
return (A*exp(-D/B)+K*g(D))*dir(i,j)
def fw(i,k):
h=norme(vec(i,k))
D=h-r
return (A*exp(-D/B)+K*g(D))*dir(i,k)
def S(i):
C,E=0,0
for k in range(2*N):
C=C+fw(X0[i],k)
for j in range(2*N):
E=E+f(X0[i],X0[j])
return m*(V-V0[i])*dir(X0[i],Z)+C+E
def F(Y,t):
X=Y[0:2*N]
V=Y[2*N:4*N]
X1=V
V1=[]
for i in range (2*N):
V1=V1+[(1/m)*S(i)]
return np.concatenate((X1,V1))
Y1=odeint(F,Y0,temps)
plt.plot(temps,Y)
plt.show() |
Partager