Système d'équations de diffusion non linéaires
Bonjour à tous (premier message sur le forum),
je cherche à résoudre un problème physique par des méthodes d'analyse numérique :
Je dispose au départ de N équations aux dérivées partielles sur mes inconnues Ci(x,t) (concentrations de l'espèce i à la position x à l'instant t) qui ressemblent fortement à des équations de diffusion, auxquelles on ajoute des termes non linéaires "en plus" pour modéliser certains phénomènes.
J'ai un peu de mal à écrire les équations sur ordi mais je tente quand meme :
d(Ci)/dt = Di * d^2(Ci)/dx^2 + f(x) - Di*Si(x)*Ci + [Somme sur (n+p = i)]k(n,p)*Cn*Cp + k(i,i)*Ci+1
- [Somme sur (n+i = p)]k(n,i)*Ci*Cn - 2k(i,i)*(Ci)^2 - k(i,i)*Ci
la fonction f(x) et Si(x), les coefficients Di et k(i,i) sont connus et fixés pour tous les i.
Bref ca ressemble pas mal à une équation de diffusion, le sommes non linéaires à droite transformant ca en "gros bordel"...
J'ai vu les conseils et liens vers les cours d'algorythmique et je suis en train d'essayer de résoudre ca en langage Fortran (peu importe le langage il me semble non ?), mais j'aurais aimé votre avis sur ma méthode de résolution (j'ai fait des choix de méthodes mais un peu subjectifs, peut être que d'autres méthodes sont plus adaptées) :
J'ai d'abord choisi d'utiliser des méthodes de différence finies.
en gros je commence par discrétiser l'espace en L points (longueur du pas non constante pour des raisons physiques - pas plus fin au bord de l'intervalle- mais pour simplifier je note le pas spatial h) , puis j'utilise la méthode "différence centrale" pour écrire
d^2(Ci)/dx^2 sous la forme (Ci,p+1 + Ci,p-1 - 2Ci,p)/h^2 , p désignant l'espace. En intégrant je transforme les autres concentrations en Ci,p*h. J'obtiens donc N*L équations différentielles du premier ordre (ODE)
Ensuite pour la variable temps, j'utilise la méthode de Crank- Nicholson qui semble adaptée aux procédés de diffusion (l'est elle ici ?)
Je mets alors, pour chaque i, l'ensemble des Ci,p sous la forme d'un vecteur spatial Ci, ce qui me donne N équations :
A*Ci(t+dt) = F1(Ci(t)) + F2(C0(t+dt),...,CN(t+dt))
F1 étant relativement simple, F2 exprimant sous forme vectorielle le gros bordel de somme vu précédemment...
A est une matrice tridiagonale.
Cette équation ressemblant fortement à une équation de neutronique, j'utilise une méthode appelée "source itération" : Cela consiste à résoudre ceci pour tous les especes i avec des termes explicites (en t) dans F2, à l'aide d'un algorythme de résolution de systèmes linéaires (algorythme de Thomas pour moi, je ne sais pas si c'est le mieux...), ce qui me donne les nouvelles concentrations Ci(t+dt), puis comme c'est inexact je re-résoud le système en remplacant les concentrations dans F2 par mes nouvelles concentrations en t +dt, et ainsi de suite jusqu'à ce que l'erreur entre le terme de gauche A*Ci et le terme de droite F1 + F2 soit très faible.
Je passe ensuite au pas de temps suivant, et ainsi de suite...
Malheureusement en lisant divers cours d'analyse numérique sur le net je ne trouve cette méthode "source itération" nulle part, pourtant je sais qu'elle est viable pour les équations multigroupes en neutronique qui ressemblent beaucoup à ce problème, pensez vous qu'elle est viable ici ?
Je pourrais aussi utiliser la méthode de Newton-Raphson pour résoudre ce système (il me semble que c'est ce que font les solveurs non ?) , est ce mieux adapté ?
Voila je suis ouvert à toute remarque et conseil sur n'importe quelle partie de la résolution de ce système, voire à des propositions de méthodes complètement différentes pour améliorer ce code. A priori je suis limité par la puissance de calcul (N peut atteindre des milliers !) pour l'instant, et je ne sais même pas si les solutions obtenues sont précises...
Merci pour votre aide !
Répondre avec citation
Partager