Forum de mathématiques - Bibm@th.net

mati

Membre

Inscription : 16-05-2018

Messages : 133

[FF++] Question sur un bout de code

Bonjour,

j'ai l'algorithme suivant qui résout le système suivant de trois équations à trois inconnues U, I, V

$$\begin{cases}dU/dt&=-a UV\\ dI/dt&=-aUV - \sigma_1 I\\ dV/dt &= D d^2V/dx^2+aI(x,t-\tau) - \sigma_2 V\end{cases}$$

avec la condition sur I:
$$I(x,t-\tau) = 0 \ \mbox{si} \ t  < \tau \ \mbox{ et } \ I(x,0)=I_0$$

(il est noté SYSTEM dans l'algorithme où $a,D, \sigma_1, \sigma_2$ sont des constantes.)
Je cherche à comprendre comment on définit $I(x,t-\tau) $ à chaque itération sur t, et je vous remercie d'avance

Voici le code écrit en FF++

macro model 2
macro typeFE() P1//
//macro DBC true// by default NBC

verbosity=0;
bool delay=true; //true si tau n'est pas nul
//load "UMFPACK64"

// Parameters
// k2=a, k3=beta, sigma2=sigma, sigma1=b
real L=2.;
real T=50.;
real a=10.;
real b=0.;
real D=1.e-3;
real beta=5.;
real sigma=0.4;
int Nbx=1e4;
real dt=.01;
real tau=10.;
real u0=1.;
real i0=0.;
real v0=10.;
real x0=0.1;
int Nbt=1e2;

real Dx=L/Nbx;

load "msh3"
meshL Th=segment(Nbx,[x*L,0.]);
fespace Vh2(Th,typeFE);
Vh2 Id = (x<=x0) ? v0 : 0.;

IFMACRO(model,2)
   Vh2 V,I,U,VH,IH,UH,V0=Id,I0=i0,U0=u0;
   macro SYSTEM
   {
      solve systemV(V,VH,solver=LU)=
         int1d(Th)(V*VH/dt + D*dx(V)*dx(VH) + sigma*V*VH)
         - int1d(Th)(V0*VH/dt + beta*Itmtau*VH)
      IFMACRO(DBC)
         + on(1,V=v0)
      ENDIFMACRO
         ;
      solve systemU(U,UH,solver=LU)=
         int1d(Th)(U*UH/dt + a*U*V*UH)
         - int1d(Th)(U0*UH/dt);
      solve systemI(I,IH,solver=LU)=
         int1d(Th)(I*IH/dt + b*I*IH)
         - int1d(Th)(I0*IH/dt + a*U*V*IH);
      V0=V; I0=I; U0=U;
   }//
ENDIFMACRO
Vh2 Itmtau=I0,Itausave;
int ntau=tau/dt;
real[int,int] Itau(Vh2.ndof,ntau+2);
Itausave=I0;
Itau(:,0)=Itausave[];
int op=0,ittau=0,ittmtau=0;

int it,itsave;
int Thnv = Th.nv;
fespace Vh(Th,P1);
Vh Vs=x,Is=x,Us=x;
real[int,int] Usave(Thnv,3),Vsave(Thnv,3),Isave(Thnv,3);
Usave(:,itsave)=Us[];
Isave(:,itsave)=Is[];
Vsave(:,itsave)=Vs[];
for (real t=0;t<T;t+=dt){
   if(t<tau){
      if(delay){
         Itmtau=0;
      } else {
         Itmtau=I0;
      }
      SYSTEM;
      if(delay){
         ittau++;//l'indice de la ligne d'après dans Itau
         Itausave=I;
         Itau(:,ittau)=Itausave[];
      }
   } else {
      if(delay){
         if(ittau==(ntau+1)){
            ittau=0;
            ittmtau=0;
         }
         Itmtau[] = Itau(:,ittmtau);
         SYSTEM;
         ittmtau++;
         Itausave=I;
         Itau(:,ittau)=Itausave[];
         ittau++;
      } else {
         Itmtau=I0;
         SYSTEM;
      }
   }
 
   if(!(it%max(1.,1./dt))){
      Vs=V;
      Is=I;
      Us=U;
      itsave++;
      Usave(:,itsave)=Us[];
      Isave(:,itsave)=Is[];
      Vsave(:,itsave)=Vs[];
      Usave.resize(Thnv,itsave+3);
      Isave.resize(Thnv,itsave+3);
      Vsave.resize(Thnv,itsave+3);
   }
   if(!(it%(max(1.,1./dt/2.)))){
      real Vmax,Vmin,Umax,Umin,Imax,Imin;
      Vs=V;
      Vmax=Vs[].max; Vmin=Vs[].min;
      Vs[]=Vs[]/Vmax;
      Is=I;
      Imax=Is[].max; Imin=Is[].min;
      Is[]=Is[]/Imax;
      Us=U;
      Umax=Us[].max; Umin=Us[].min;
      Us[]=Us[]/Umax;
      meshL ThmvV=movemesh(Th,[x/L,y+Vs]);//pour normaliser on met x/L (le domaine est entre 0 et 1), x/l*2 entre 0 et 2.
      meshL ThmvI=movemesh(Th,[x/L,y+Is]);
      meshL ThmvU=movemesh(Th,[x/L,y+Us]);
      plot(ThmvV,ThmvI,ThmvU,dim=2,wait=0,cmm="time: "+t+", U min = "+Umin+", U max = "+Umax+", V min = "+Vmin+", V max = "+Vmax+", I min = "+Imin+", I max = "+Imax);
   }
   op=1;
   it++;
}

ofstream foutU("u.txt");
ofstream foutI("i.txt");
ofstream foutV("v.txt");
for(int m=0;m<Th.nv;m++) {
   for(int n=0;n<itsave;n++) {
      foutV << Vsave(m,n) << " ";
      foutI << Isave(m,n) << " ";
      foutU << Usave(m,n) << " ";
   }
   foutV<<endl;
   foutI<<endl;
   foutU<<endl;
}
   

Cordialement

Dernière modification par yoshi (09-04-2021 13:31:58)

yoshi

Modo Ferox

Inscription : 20-11-2005

Messages : 16 988

Re : [FF++] Question sur un bout de code

Bonjour,

J'ignorais l'existence de ce langage, apparemment après recherches :
- ce serait un dérivé du C++,
- Google ignore l'existence d'un forum d'entraide consacré au FF++,
- tout ce que j'ai pu trouver est ici :
  https://freefem.org/ff++
  http://www.cmap.polytechnique.fr/IMG/pd … eFem_I.pdf
  http://www.georges-sadaka.fr/FreeFem++/ … m++2D.html

Désolé...

@+

---

Arx Tarpeia Capitoli proxima...

LEG

Membre

Inscription : 19-09-2012

Messages : 694

Re : [FF++] Question sur un bout de code

Bonjour

Dans quel domaine est utilisée cet algorithme et dans quel but...? ainsi que tes formules , peut être qu'il faudrait regarder manuellement ce que cela donne puis le faire en python...
Avant de demander à quoi sert un programme que tu as récupéré quelque part...???

Car il semblerait que tu n'ais pas beaucoup de réponse ''sur les sites de mathématiques''