Besoin de comprendre la logique de cet algo — Les-mathematiques.net The most powerful custom community solution in the world

Besoin de comprendre la logique de cet algo

Bonjour,
j'ai l'algorithme suivant qui définit la fonction $I$ qui véfifie la condition
$$I(x,t-\tau) = 0 \quad \mbox{si} \quad t < \tau \qquad \mbox{ et } \qquad I(x,0)=I_0,
$$ puis qui résout le système de trois équations à trois inconnues $U, I, V$ suivant :
$$
dU/dt=-a UV, \quad dI/dt=-aUV - \sigma_1 I, \quad dV/dt = D d^2V/dx^2+aI(x,t-\tau) - \sigma_2 V.

$$ (il est noté SYSTEM dans l'algorithme où $a,D, \sigma_1, \sigma_2$ sont des constantes.
L'algorithme est le suivant.
int ntau=tau/dt;
Itausave=I0;
Itau(:,0)=Itausave[];
int ittau=0,ittmtau=0;

for (real t=0;t<T;t+=dt){
   if(t<tau){
      if(tau>0){
         Itmtau=0;
      } else {
         Itmtau=I0;
      }
      SYSTEM;
      if(tau>0){
         ittau++;//l'indice de la ligne d'après dans Itau
         Itausave=I;
         Itau(:,ittau)=Itausave[];
      }
   } else {
      if(tau>0){
         if(ittau==(ntau+1)){
            ittau=0;
            ittmtau=0;
         }
         Itmtau[] = Itau(:,ittmtau);
         SYSTEM;
         ittmtau++;
         Itausave=I;
         Itau(:,ittau)=Itausave[];
         ittau++;
      } else {
         Itmtau=I0;
         SYSTEM;
      }
   }
Est-ce que quelqu'un peut m'expliquer ce que ce bout de programme fait ?
Cordialement.

Réponses

  • Bonjour,
    on ne devrait pas avoir à poser ce genre de question car un tel programme ne devrait pas exister. Pas un seul commentaire. Peut-être y a-t-il un en-tête mais ce n'est pas sûr.
    Cordialement.
  • Bonjour.

    Si, il y en a un, par contre pas de chance, il explique une chose qui n'a pas lieu d'être expliquée comparativement à d'autres opérations qui sont fournies sans explications.

    Effectivement, c'est un très mauvais exemple de code :


    _ Quelque déclarations de variables, mais loin d'être suffisant.
    _ Nommage non homogène des variables
    _ Beaucoup d'opérations cachées dans des boucles implicites. Au passage, je me demande si cette machinerie tourne vraiment. J'aimerais avoir un exemple d'instanciation.
    _ Utilisation outrancière de 'if then', peut-être légitime, mais cela ne participe pas à la lisibilité et cela peut très certainement être amélioré, seulement, en l'état ce serait compliqué sans plus d'explication sur la méthode d'origine (je suppose du Euler, sans certitude).

    À bientôt.

    Cherche livres et objets du domaine mathématique :

    Intégraphes, règles log et calculateurs électromécaniques.

  • Bonjour,

    Déjà, contrairement à ce qu'écrit Mati, ce n'est pas un algorithme, c'est un bout de programme, écrit dans on ne sait quel langage, isolé de son contexte et de ses racines, avec des variables non déclarées dès la première ligne, etc...
    Bref, c'est illisible.
    Si un étudiant m'avait rendu une telle cochonnerie, c'était $0$ sans plus rentrer dedans.

    Cordialement,

    Rescassol
  • Bonjour,
    je m'excuse je pensais bien faire en n'envoyant que le bout de code qui concerne ma question bien précise à savoir comment programmer I(x,t-\tau) à chaque itératrion sur t.
    Le programme est écrit en FF++ et le voici au complet. J'espère que quelqu'un peut m'aider
    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;
    }
    
  • C'est un peu mieux, mais tu es toujours aussi avare d'explications ?

    Il y a juste trois lignes commentées, dont deux qui sont franchement accessoires.

    Cordialement.

    Cherche livres et objets du domaine mathématique :

    Intégraphes, règles log et calculateurs électromécaniques.

  • Mais si je pose la question c'est parce que justement je cherche ces explication. Notamment à propos de la boucle sur if (delay) c'est à dire si tau > 0
Connectez-vous ou Inscrivez-vous pour répondre.
Success message!