Forum Programmation.python Equation de réaction-diffusion modélisation

-1
25
mai
2016

Bonjour à tous,
je souhaite observer des structures semblables à des zébrures avec le programme suivant, qui concerne un système de deux équations de réaction diffusion. Les seules "taches" obtenues sont celles d'un guépard, comment faire varier les paramètres pour avoir un motif de zèbre ?
Il est possible de faire varier a et b les coefficients de diffusion, la longueur de la grille, le pas… mais aussi les formules qui suivent le laplacien dans la dernière boucle for … savez vous comment faire pour obtenir des zébrures ?
Merci d'avance

import scipy as sp
from matplotlib import *
from pylab import *

Paramètres de diffusion

a =0.0000000008
b =0.0000007
size = 100 # taille de la grille en 2D
dx = 2./size # intervalle d'espace
T = 10.0 # temps total
dt = .8 * dx**2/2 # pas de la variable de temps
n = int(T/dt)

initialisation de u et v

u = np.random.rand(size, size)
v = np.random.rand(size, size)

def laplacian(Z):
Ztop = Z[0:-2,1:-1]
Zleft = Z[1:-1,0:-2]
Zbottom = Z[2:,1:-1]
Zright = Z[1:-1,2:]
Zcenter = Z[1:-1,1:-1]
return (Ztop + Zleft + Zbottom + Zright - 4 * Zcenter) / dx**2

for i in range(n):
deltaU=laplacian(u)
deltaV=laplacian(v)
ui=u[1:-1,1:-1]
#Conditions initiales
vi=v[1:-1,1:-1]
#calcul de u[i,j] à partir de ui[i,j]; pareil pour v et vi
u[1:-1,1:-1],v[1:-1,1:-1]=ui+dt*(a*deltaU+ui*(vi*2-1)-12),vi+dt*(b*deltaV+16-vi)

plt.imshow(u,cmap=plt.cm.bone, extent=[-1,1,-1,1]);
plt.xticks([]); plt.yticks([]);
plt.colorbar()

show()

  • # Espace isotrope

    Posté par . Évalué à 5.

    Si j'ai bien compris, ton plan est isotrope. Tu n'auras jamais rien d'autre que de la symétrie centrale à moins de commencer avec une rayure. Ou alors il faut rendre ton plan anisotrope en modifiant la formule de diffusion.

  • # c'est beau avec un peu de markdown

    Posté par . Évalué à 7.

    u = np.random.rand(size, size)
    v = np.random.rand(size, size)
    
    def laplacian(Z):
      Ztop = Z[0:-2,1:-1]
      Zleft = Z[1:-1,0:-2]
      Zbottom = Z[2:,1:-1]
      Zright = Z[1:-1,2:]
      Zcenter = Z[1:-1,1:-1]
      return (Ztop + Zleft + Zbottom + Zright - 4 * Zcenter) / dx**2
    
    for i in range(n):
      deltaU=laplacian(u)
      deltaV=laplacian(v)
      ui=u[1:-1,1:-1]
      #Conditions initiales
      vi=v[1:-1,1:-1]
      #calcul de u[i,j] à partir de ui[i,j]; pareil pour v et vi
      u[1:-1,1:-1],v[1:-1,1:-1]=ui+dt*(a*deltaU+ui*(vi*2-1)-12),vi+dt*(b*deltaV+16-vi)
    
    plt.imshow(u,cmap=plt.cm.bone, extent=[-1,1,-1,1]);
    plt.xticks([]); plt.yticks([]);
    plt.colorbar()
    
    show()

Suivre le flux des commentaires

Note : les commentaires appartiennent à ceux qui les ont postés. Nous n'en sommes pas responsables.