equations de reaction/diffusion

Damienrobert

equations de reaction/diffusion

Message par Damienrobert » 27 août 2016 16:40

Bonjour à tous. Dans le cadre de mon TIPE je cherche à coder un algorithme python pour modéliser la création des motifs sur les pelages des animaux. J’ai essayé de coder par moi-même, soit avec des fonctions python, soit avec un schéma numérique… mais cela ne donne pas grand-chose (mon « pelage » reste uni sans motif). J’ai donc trouvé cet algorithme en MAPLE (je n’y connais rien du tout) qui normalement fonctionne mais je n’ai pas réussi à le traduire correctement. Quelqu’un saurait-il me le transcrire en python ? merci beaucoup

> with(linalg) ;
Les fonctions de reaction :
> f :=(u,v,a) -> a-u+u^2*v ;
> g :=(u,v,b)->b-u^2*v ;
La solution homogene stationnaire (u_eq,v_eq) :
> s :=solve({f(u,v,a)=0,g(u,v,b)=0},{u,v}) ;
> u_eq :=(a,b)->a+b ;
> v_eq :=(a,b)->b/((a+b)^2) ;
La matrice jacobienne au point d’équilibre : (on la définit d’abord comme une liste de listes)
> AA :=(a,b,gamma,d,lambda)->[[gamma*D[1](f)(u_eq(a,b),v_eq(a,b),a,b) –lambda, gamma*D[2](f)(u_eq(a,b),v_eq(a,b),a,b)], [gamma*D[1](g)(u_eq(a,b),v_eq(a,b)),gamma*D[2](g)(u_eq(a,b),v_eq(a,b)-lambda*d]] ;
Pour obtenir une matrice on utilise convert
> A :=(a,b,gamma,d,lambda)->convert(AA(a,b,gamma,d,lambda),matrix) ;
La valeur critique d_c est solution de l’equation quadratique
> q :=(a,b,d)->d^2*(D[1](f)(u_eq(a,b),a,b))^2+2*(2*D[2](f)(u_eq(a,b),v_eq(a,b),a,b)*D[1](g)(u_eq(a,b),v_eq(a,b))-D[1](f)(u_eq(a,b),v_eq(a,b),a,b)*D[2](g)(u_eq(a,b),v_eq(a,b)))*d+(D[2](g)(u_eq(a,b),v_eq(a,b)))^2 ;
> solve(q(0.25,0.75,d) ;
> eigenvects(A(0.25,0.75,8,20,(1*Pi^2)/1)) ;
On choisit des valeurs des parameters a, b, gamma et d, puis on fait le graphe du determinant en function de lambda
> plot(det(A(0.25,0.75,8,20,lambda)),lambda=0.4..2.5) ;
On trouve les valeurs pour les different modes sur un domaine rectangulaire
> tableau :=proc(m_min,m_max,n_min,n_max,p,q)
local h,i,j,k ;
h :=i->[seq(evalf(Pi^2*(k^2/(p^2)+i^2/(q^2)),3),k=m_min..m_max)] ;
array([seq(h(j),j=n_min..n_max)]) ;
end :
> tableau(0,10,0,10,5,2) ;

Pour mieux comprendre, voice les equations en jeu:
équations de recation-diffusion.png
équations de recation-diffusion.png (11.89 Kio) Consulté 2678 fois

Messages : 9679

Inscription : 30 juil. 2008 16:59

Profil de l'utilisateur : Élève de lycée

Re: equations de reaction/diffusion

Message par fakbill » 27 août 2016 17:47

Non.
Ce qu'il te faut c'est comprendre ce qui ne va pas dans ce que tu fais.
Tenter de reproduire des algos au pif ne vas pas t'aider.
Tu cherches à répondre à quelle question? Quels paramètres donnent un état intéressant (pas tout plat)? Si oui, commence par faire résoudre tes équations par un truc dont on sait qu'il marche PUIS, si tu as fini ton étude paramétrique, essaye de coder toi même le schéma.
Tu as codé des schémas. ok. Les as tu testés sur des exemples simples? ils fonctionnent?

"soit avec des fonctions python, soit avec un schéma numérique" : gni??? un schéma numérique est un algo qu'on peut coder en ce qu'on veut.
Pas prof.
Prépa, école, M2, thèse (optique/images) ->ingé dans le privé.

Damienrobert

Re: equations de reaction/diffusion

Message par Damienrobert » 28 août 2016 16:58

J'ai bien compris comment fonctionne cet algo, qui j'ai l'impression ressemble à peu près au mien, le but de ma question était donc d'avoir d'abord un algo qui marche afin de résoudre le mien.
Le fait est que mon algo ne cherche pas résoudre un problème, mais à modéliser l'apparition de motifs sur un pelage à partir de conditions initiales données... il n'y a donc pas de solution basique pour vérifier qu'il fonctionne (à part peut-être tout nul sur toute la surface comme conditions initiales). Le problème est de plus qu'il y a plusieurs facteurs pour aboutir à quelque chose (les concentrations de u et de v au départ, le temps pendant lequel on le fait tourner, la taille de la surface...) et je ne sais pas lesquels prendre, pour ça que j'essaye de prendre ceux de cet algorithme.
Pour résoudre les équations, il faut connaître les fonctions python qui trouvent les solutions de système d'équations différentielles partielles non linéaires couplées... que je n'ai pas trouvé même en cherchant dans la bibliothèque numpy.
J'ai donc essayé d'emboîter trois boucles, de discrétiser les dérivées et d'utiliser des matrices... mais les lignes de calculs renvoient des résultats nan (not a number) ce que je ne comprends pas. (j'ai rajouté des print pour déceler cette erreur)
C'est donc parce qu'il me semblait difficile d'arriver à un algorithme exploitable vu le nombre de paramètres que je demandais la traduction cet algorithme maple (qui est normalement correct). ;-)
Merci en tout cas pour cette première réponse. Voici ci-dessous mon essai d'algo :-)

import numpy as np
import matplotlib.pyplot as plt

a = 0.25
b=0.75
d=20
Lx = 10
Ly = 10
T = 5.0
NX = 100
NY = 100
NT = 3
dx = Lx/(NX-2)
dy = Ly/(NY-2)
dt = T/NT

U = np.zeros((NX,NY))
V = np.zeros((NX,NY))

for i in range(0,NX):
for j in range(0,NY):
U[i,j] = 0
V[i,j] = 0

DU = np.zeros((NX,NY))
DV = np.zeros((NX,NY))

for t in range(0,NT):
for i in range(2,NX-1):
for j in range(2,NY-1):
DU[i,j] , DV[i,j] = dt*(a - U[i,j] + U[i,j]*U[i,j]*V[i,j]) + dt*(U[i-1,j] -2*U[i,j] + U[i+1,j])/(dx**2) + dt*(U[i,j-1] -2*U[i,j] + U[i,j+1])/(dy**2) , dt*(b - U[i,j]*U[i,j]*V[i,j]) + d*dt*(V[i-1,j] -2*V[i,j] + V[i+1,j])/(dx**2) + d*dt*(V[i,j-1] -2*V[i,j] + V[i,j+1])/(dy**2)
U[i,j] += DU[i,j]
V[i,j] += DV[i,j]

def couleur(u,v):
if u>v:
return('.r')
elif v>u:
return('.g')
else:
return('.b')

for x in range(0,NX):
for y in range(0,NY):
plt.plot([x],[y],couleur(U[x][y],V[x][y]))

plt.show()

Messages : 9679

Inscription : 30 juil. 2008 16:59

Profil de l'utilisateur : Élève de lycée

Re: equations de reaction/diffusion

Message par fakbill » 28 août 2016 17:26

Heu si ça existe c'est peut être dans scipy mais pas dans numpy.
dx = Lx/(NX-2) : boum! en python2 ou en python3? Dans tous les cas, c'est une mauvaise idée de diviser deux entiers et de laisser implicitement le langage faire la conversion en float ou pas.

"C'est donc parce qu'il me semblait difficile d'arriver à un algorithme exploitable vu le nombre de paramètres que je demandais la traduction cet algorithme maple (qui est normalement correct). ;-)"
QUand tu dis ça, tu confonds l'algo et ses paramètres.

Le maple est assez simple.
tu as besoin de https://www.maplesoft.com/support/help/ ... th=convert
et de solve. Le reste est clair non?

Par contre, ton code maple fait je ne sais trop quoi (je n'ai pas réfléchi) mais ce n'est pas un schéma numérique de résolution d'EDP comme ce que tu essayes de faire en python/scipy/numpy.
Pas prof.
Prépa, école, M2, thèse (optique/images) ->ingé dans le privé.

Messages : 9679

Inscription : 30 juil. 2008 16:59

Profil de l'utilisateur : Élève de lycée

Re: equations de reaction/diffusion

Message par fakbill » 28 août 2016 17:31

Tu as regardé un peu avec google??
J'ai trouvé ça :
https://github.com/thomasdeniau/pyfauxfur
Si tu veux le faire marcher il va te falloir installer qt, opengl et weave (sous windows ça risque d'être difficile...je ne sais pas).
Cependant, tout est là :
https://github.com/thomasdeniau/pyfauxf ... ageData.py entre la ligne 133 et 177 (c'est du code C que le module weave va compiler et lier ça au code python à la volée. Ca permet d'avoir de meilleures perfo qu'un code en pur python....cela dit on peut l'écrire en numpy et, si c'est bien écrit, il est probable que les perfs ne soient pas si mauvaises. Bref tu t'en fiche...le code est facile à lire).
Pas prof.
Prépa, école, M2, thèse (optique/images) ->ingé dans le privé.

Messages : 9679

Inscription : 30 juil. 2008 16:59

Profil de l'utilisateur : Élève de lycée

Re: equations de reaction/diffusion

Message par fakbill » 28 août 2016 17:43

Pas prof.
Prépa, école, M2, thèse (optique/images) ->ingé dans le privé.

Damienrobert

Re: equations de reaction/diffusion

Message par Damienrobert » 30 août 2016 21:43

Merci pour tous ces liens.
Effectivement certains ne fonctionnent pas sur mon ordinateur car les modules (notamment fipy.random.uniform) ne sont pas présents.
J'ai aussi essayé de faire un nouveau code avec les fonctions fipy... mais sans grand succès, peut-être que mon prof de maths saura l'améliorer.
Merci encore
Damien :)

Répondre