Besoin de comprendre un bout de code

**aicha84** · 08/04/2021, 22h06

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

$Formule mathématique$

avec la condition sur I:
$Formule mathématique$

(il est noté SYSTEM dans l'algorithme où $Formule mathématique$ sont des constantes.)
Je cherche à comprendre comment on définit $Formule mathématique$ à chaque itération sur t, et je vous remercie d'avance

Voici le code écrit en FF++

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
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;
}