# Calculer le max sur (x,y) d'une fonction (x,y,t)

Bonjour
j’ai le code suivant et je souhaite calculer
eh(t)= max_{x,y} |uh(x,y,t)-wh(x,y,t)|
où uh et wh sont deux solutions de deux problèmes différents calculer dans le code. Le voici joint et je vous remercie par avance pour toute aide.

``````load "UMFPACK64"
defaulttoUMFPACK64();

verbosity=1;

func real qdelta (real t, real s, real w) {
real t1 = t - rint(t);
real s1 = s - rint(s);
real delta=1.0/sqrt(w*pi);
if ( square(s1)+square(t1) <= square(delta))
return w;
else
return 0.0;
};

func real q (real t, real s, real w) {
real t1 = t - rint(t);
real s1 = s - rint(s);
real delta=1.0/sqrt(w*pi);
if ( square(s1)+square(t1) <= square(delta))
return 1;
else
return 0.0;
};

func real f (real x)
{
return 0.25*sin(2*pi*x);
};

// initial values for the spatial domain
int N=200;
real eps=0.5, w=3.0;

real t=0.00;
real dt=0.1;

int ni=5;       // number of values of h=eps^2 w
real[int] h(ni);     // value of h at each step
real[int] error(ni); // error in norm H1 t each step
real[int] alpha(ni); // alpha =1/log(2). log(error(h))/log(error(h/2))
ofstream fout("error.txt");
fout << "# step h error alpha" << endl;

for (int i=0;i<=ni-1;++i) {

mesh Th=square(N,N,[x,y]);
fespace Vh(Th,P1,periodic=[[2,y],[4,y],[1,x],[3,x]]);
fespace Wh(Th,P2,periodic=[[2,y],[4,y],[1,x],[3,x]]);
Vh uh=0.0;
Vh wh=0.0;
Vh vh;
Vh oldU=0.0;
Vh oldW=0.0;
real T=1.0;// arbitraire
int[int] Order= [1];

string DataName;

Vh fh=f;

Vh nith = 0.0; // for nonlinear iterations
Vh g0,g1,g2,g3;
Wh uplot;

// Nonlinear part is calculated in previous iteration (nith) and becomes a part of the right hand side
problem het(uh,vh,init=t) =
int2d(Th)((dx(uh)*dx(vh)+dy(uh)*dy(vh))*dt)
+int2d(Th)(uh*vh)
-int2d(Th)(oldU*vh)
+int2d(Th)(qdelta(x/eps,y/eps,w)*g0*vh*dt) +int2d(Th)(qdelta(x/eps,y/eps,w)*g1*uh*vh*dt)
-int2d(Th)(fh*vh*dt)
;

problem hom(wh,vh,init=t) =
int2d(Th)((dx(wh)*dx(vh)+dy(wh)*dy(vh))*dt)
+int2d(Th)(wh*vh)
-int2d(Th)(oldW*vh)
+int2d(Th)(g2*vh*dt) +int2d(Th)(g3*wh*vh*dt)
-int2d(Th)(fh*vh*dt)
;

uplot=q(x/eps,y/eps,w);

savevtk("q_"+i+".vtu",Th,uplot,dataname="qdelta", order=Order);
cout << "mean value of qdelta: " << int2d(Th)(uplot);

real max = -1e100; //initialisation de max
ofstream fouti("max.txt");//créer le fichiers où seront écrites les valeurs

for (t=0;t<=T;t+=dt) {
cout << "Time = " << t << endl;
fh=f(t);
oldU=uh;
oldW=wh;
int k=0;
int noit = 10; // maximum no of nonlinear iterations
real errnlhet = 0;  // error of nonlinear iteration of het
real normuh = 0;
real errnlhom =0; // error of nonlinear iteration of hom
real normwh=0;

// Nonlinear iterations for heterogeneous solution
nith = uh; // start nonlinear iterations  with last solution for het
g0=uh*uh/(1.0+uh)^2;
g1=1.0/(1.0+uh)^2;
cout << "    Solving nonlinear heterogeneous problem.\n";
for(k=1;k<noit; ++k) {
het;
normuh = sqrt(int2d(Th)(square(uh))); // L^2 nor of the solution of het
errnlhet = sqrt(int2d(Th)(square(uh - nith))); // L^2 norm of the error in nonlinear iteration of het
cout << "        " << k << ": error = " << errnlhet << endl;
if(errnlhet <= 1E-5*normuh) // stopping criteria. If relative error
break;        // is less than 1E-5 stop nonlinear iterations of het
nith = uh; // prepare new iteration
g0=uh*uh/(1.0+uh)^2;
g1=1.0/(1.0+uh)^2;
}
// Check convergence for het
if(k >= noit) {
cout << "Nonlinear solver for het did not converged.\n";
cout << "Error = " << errnlhet << endl;
cout << "|| uh || = " << normuh << endl;
exit(1); // exit the program with error code 1
}

nith= wh; // start nonlinear iteration with last solution for hom
g2=wh*wh/(1.0+wh)^2;
g3=1.0/(1.0+wh)^2;
cout << "    Solving nonlinear homogeneous problem.\n";
for(k=1;k<noit; ++k) {
hom;
normwh = sqrt(int2d(Th)(square(wh))); // L^2 nor of the solution of hom
errnlhom = sqrt(int2d(Th)(square(wh - nith))); // L^2 norm of the error in nonlinear iteration of hom
cout << "        " << k << ": error = " << errnlhom << endl;
if (errnlhom <= 1E-5*normwh) // stopping criteria. If relative error
break;        // is less than 1E-5 stop nonlinear iterations of het
nith = wh; // prepare new iteration
g2=wh*wh/(1.0+wh)^2;
g3=1.0/(1.0+wh)^2;
}
// Check convergence for het
if(k >= noit) {
cout << "Nonlinear solver for hom did not converged.\n";
cout << "Error = " << errnlhom << endl;
cout << "|| wh || = " << normwh << endl;
exit(1); // exit the program with error code 1
}

Vh eh=abs(uh-wh); // Définition de la fonction eh pour laquelle on veut calculer son max en (x,y)
real umaxt = eh[].max; //  max au temps temp
umax = max( max, umaxt);// max jusqu'au temps t
fouti << umaxt << " "<< umax << endl;
}

savevtk("u_"+i+".vtu",Th,uh,dataname="uh", order=Order);
savevtk("w_"+i+".vtu",Th,wh,dataname="vh", order=Order);

real errH1 = sqrt(int2d(Th)(square(uh - wh)) + int2d(Th)(square(dx(uh)-dx(wh)) + square(dy(uh)-dy(wh))));

cout << " ==== errH1 ( eps = " << eps <<", omega = " << w << ", eps^2*omega = " << eps^2*w << ") = "<< errH1 <<endl;

h(i)=eps^2*w;
error(i)=errH1;
alpha(i)=i>0?1.0/log(2.0)*log(error(i))/log(error(i-1)):1.0;
fout << i <<" "<< h(i) <<" "<< error(i) <<" "<< alpha(i) << endl;
fout.flush;

eps=eps/2;
w=2*w;

}
``````

See below the code with the max for each time step and the max of all time steps:

``````//load "UMFPACK64"
//defaulttoUMFPACK64();

verbosity=1;

func real qdelta (real t, real s, real w) {
real t1 = t - rint(t);
real s1 = s - rint(s);
real delta=1.0/sqrt(w*pi);
if ( square(s1)+square(t1) <= square(delta))
return w;
else
return 0.0;
}

func real q (real t, real s, real w) {
real t1 = t - rint(t);
real s1 = s - rint(s);
real delta=1.0/sqrt(w*pi);
if ( square(s1)+square(t1) <= square(delta))
return 1;
else
return 0.0;
}

func real f (real x) {
return 0.25*sin(2*pi*x);
}

// initial values for the spatial domain
int N=20;
real eps=0.5, w=3.0;

real t=0.00;
real dt=0.1;

int ni=5;       // number of values of h=eps^2 w
real[int] h(ni);     // value of h at each step
real[int] error(ni); // error in norm H1 t each step
real[int] alpha(ni); // alpha =1/log(2). log(error(h))/log(error(h/2))
ofstream fout("error.txt");
fout << "# step h error alpha" << endl;

for (int i=0;i<=ni-1;++i) {
mesh Th=square(N,N,[x,y]);
fespace Vh(Th,P1,periodic=[[2,y],[4,y],[1,x],[3,x]]);
fespace Wh(Th,P2,periodic=[[2,y],[4,y],[1,x],[3,x]]);
Vh uh=0.0;
Vh wh=0.0;
Vh vh;
Vh oldU=0.0;
Vh oldW=0.0;
real T=1.0;// arbitraire
int[int] Order= [1];

string DataName;

Vh fh=f;

Vh nith = 0.0; // for nonlinear iterations
Vh g0,g1,g2,g3;
Wh uplot;

// Nonlinear part is calculated in previous iteration (nith) and becomes a part of the right hand side
problem het(uh,vh,init=t) =
int2d(Th)((dx(uh)*dx(vh)+dy(uh)*dy(vh))*dt)
+int2d(Th)(uh*vh)
-int2d(Th)(oldU*vh)
+int2d(Th)(qdelta(x/eps,y/eps,w)*g0*vh*dt) +int2d(Th)(qdelta(x/eps,y/eps,w)*g1*uh*vh*dt)
-int2d(Th)(fh*vh*dt)
;

problem hom(wh,vh,init=t) =
int2d(Th)((dx(wh)*dx(vh)+dy(wh)*dy(vh))*dt)
+int2d(Th)(wh*vh)
-int2d(Th)(oldW*vh)
+int2d(Th)(g2*vh*dt) +int2d(Th)(g3*wh*vh*dt)
-int2d(Th)(fh*vh*dt)
;

uplot=q(x/eps,y/eps,w);

savevtk("q_"+i+".vtu",Th,uplot,dataname="qdelta", order=Order);
cout << "mean value of qdelta: " << int2d(Th)(uplot);

real umax = -1e100; //initialisation de max
//ofstream fouti("max.txt");//créer le fichiers où seront écrites les valeurs
// je retire se fichier car il est écrasé à chaque boucle sur i

for (t=0;t<=T;t+=dt) {
cout << "Time = " << t << endl;
fh=f(t);
oldU=uh;
oldW=wh;
int k=0;
int noit = 10; // maximum no of nonlinear iterations
real errnlhet = 0;  // error of nonlinear iteration of het
real normuh = 0;
real errnlhom =0; // error of nonlinear iteration of hom
real normwh=0;

// Nonlinear iterations for heterogeneous solution
nith = uh; // start nonlinear iterations  with last solution for het
g0=uh*uh/(1.0+uh)^2;
g1=1.0/(1.0+uh)^2;
cout << "    Solving nonlinear heterogeneous problem.\n";
for(k=1;k<noit; ++k) {
het;
normuh = sqrt(int2d(Th)(square(uh))); // L^2 nor of the solution of het
errnlhet = sqrt(int2d(Th)(square(uh - nith))); // L^2 norm of the error in nonlinear iteration of het
cout << "        " << k << ": error = " << errnlhet << endl;
if(errnlhet <= 1E-5*normuh) // stopping criteria. If relative error
break;        // is less than 1E-5 stop nonlinear iterations of het
nith = uh; // prepare new iteration
g0=uh*uh/(1.0+uh)^2;
g1=1.0/(1.0+uh)^2;
}
// Check convergence for het
if(k >= noit) {
cout << "Nonlinear solver for het did not converged.\n";
cout << "Error = " << errnlhet << endl;
cout << "|| uh || = " << normuh << endl;
exit(1); // exit the program with error code 1
}

nith= wh; // start nonlinear iteration with last solution for hom
g2=wh*wh/(1.0+wh)^2;
g3=1.0/(1.0+wh)^2;
cout << "    Solving nonlinear homogeneous problem.\n";
for(k=1;k<noit; ++k) {
hom;
normwh = sqrt(int2d(Th)(square(wh))); // L^2 nor of the solution of hom
errnlhom = sqrt(int2d(Th)(square(wh - nith))); // L^2 norm of the error in nonlinear iteration of hom
cout << "        " << k << ": error = " << errnlhom << endl;
if (errnlhom <= 1E-5*normwh) // stopping criteria. If relative error
break;        // is less than 1E-5 stop nonlinear iterations of het
nith = wh; // prepare new iteration
g2=wh*wh/(1.0+wh)^2;
g3=1.0/(1.0+wh)^2;
}
// Check convergence for het
if(k >= noit) {
cout << "Nonlinear solver for hom did not converged.\n";
cout << "Error = " << errnlhom << endl;
cout << "|| wh || = " << normwh << endl;
exit(1); // exit the program with error code 1
}

Vh eh=abs(uh-wh); // Définition de la fonction eh pour laquelle on veut calculer son max en (x,y)
real umaxt = eh[].max; //  max au temps temp
cout << "Umax for time t = " << umaxt << endl;
//fouti << umaxt << " "<< umax << endl;

umax = max(umax, umaxt);// max jusqu'au temps t
cout << "Global max = " << umax << endl;
}

cout << "Final max = " << umax << endl;

savevtk("u_"+i+".vtu",Th,uh,dataname="uh", order=Order);
savevtk("w_"+i+".vtu",Th,wh,dataname="vh", order=Order);

real errH1 = sqrt(int2d(Th)(square(uh - wh)) + int2d(Th)(square(dx(uh)-dx(wh)) + square(dy(uh)-dy(wh))));

cout << " ==== errH1 ( eps = " << eps <<", omega = " << w << ", eps^2*omega = " << eps^2*w << ") = "<< errH1 <<endl;

h(i)=eps^2*w;
error(i)=errH1;
alpha(i)=i>0?1.0/log(2.0)*log(error(i))/log(error(i-1)):1.0;
fout << i <<" "<< h(i) <<" "<< error(i) <<" "<< alpha(i) << endl;
fout.flush;

eps=eps/2;
w=2*w;
}
``````

Bonjour
merci beaucoup pour la réponse. En fait moi j’ai besoin de calculer
eh(t_k)= max_{x,y} |uh(x,y,t_k)-wh(x,y,t_k)| pour t_k de l’intervalle de temps. Puis mettre les data dans un fichier afin de dessiner eh en fonction du temps t. On refait tout ça pour chaque i: à chaque i un nouveau fichier où il y a les data de eh(t_k). Voici joint comment j’ai modifié le code. C’est correcte? S’il vous plaît

``````
defaulttoUMFPACK64();

verbosity=1;

func real qdelta (real t, real s, real w) {
real t1 = t - rint(t);
real s1 = s - rint(s);
real delta=1.0/sqrt(w*pi);
if ( square(s1)+square(t1) <= square(delta))
return w;
else
return 0.0;
};

func real q (real t, real s, real w) {
real t1 = t - rint(t);
real s1 = s - rint(s);
real delta=1.0/sqrt(w*pi);
if ( square(s1)+square(t1) <= square(delta))
return 1;
else
return 0.0;
};

func real f (real x)
{
return 0.25*sin(2*pi*x);
};

// initial values for the spatial domain
int N=200;
real eps=0.5, w=3.0;

real t=0.00;
real dt=0.1;

int ni=5;       // number of values of h=eps^2 w
real[int] h(ni);     // value of h at each step
real[int] error(ni); // error in norm H1 t each step
real[int] alpha(ni); // alpha =1/log(2). log(error(h))/log(error(h/2))
ofstream fout("error.txt");
fout << "# step h error alpha" << endl;

for (int i=0;i<=ni-1;++i) {

ofstream fouti("max_"+i+".txt");
fouti << "# t max" << endl;

mesh Th=square(N,N,[x,y]);
fespace Vh(Th,P1,periodic=[[2,y],[4,y],[1,x],[3,x]]);
fespace Wh(Th,P2,periodic=[[2,y],[4,y],[1,x],[3,x]]);
Vh uh=0.0;
Vh wh=0.0;
Vh vh;
Vh oldU=0.0;
Vh oldW=0.0;
real T=1.0;// arbitraire
int[int] Order= [1];

string DataName;

Vh fh=f;

Vh nith = 0.0; // for nonlinear iterations
Vh g0,g1,g2,g3;
Wh uplot;

// Nonlinear part is calculated in previous iteration (nith) and becomes a part of the right hand side
problem het(uh,vh,init=t) =
int2d(Th)((dx(uh)*dx(vh)+dy(uh)*dy(vh))*dt)
+int2d(Th)(uh*vh)
-int2d(Th)(oldU*vh)
+int2d(Th)(qdelta(x/eps,y/eps,w)*g0*vh*dt) +int2d(Th)(qdelta(x/eps,y/eps,w)*g1*uh*vh*dt)
-int2d(Th)(fh*vh*dt)
;

problem hom(wh,vh,init=t) =
int2d(Th)((dx(wh)*dx(vh)+dy(wh)*dy(vh))*dt)
+int2d(Th)(wh*vh)
-int2d(Th)(oldW*vh)
+int2d(Th)(g2*vh*dt) +int2d(Th)(g3*wh*vh*dt)
-int2d(Th)(fh*vh*dt)
;

uplot=q(x/eps,y/eps,w);
//plot(uplot,fill=true,value=true,ps="q".eps");
savevtk("q_"+i+".vtu",Th,uplot,dataname="qdelta", order=Order);
cout << "mean value of qdelta: " << int2d(Th)(uplot);

real umax = -1e100; //initialisation de max

for (t=0;t<=T;t+=dt) {
cout << "Time = " << t << endl;
fh=f(t);
oldU=uh;
oldW=wh;
int k=0;
int noit = 10; // maximum no of nonlinear iterations
real errnlhet = 0;  // error of nonlinear iteration of het
real normuh = 0;
real errnlhom =0; // error of nonlinear iteration of hom
real normwh=0;

// Nonlinear iterations for heterogeneous solution
nith = uh; // start nonlinear iterations  with last solution for het
g0=uh*uh/(1.0+uh)^2;
g1=1.0/(1.0+uh)^2;
cout << "    Solving nonlinear heterogeneous problem.\n";
for(k=1;k<noit; ++k) {
het;
normuh = sqrt(int2d(Th)(square(uh))); // L^2 nor of the solution of het
errnlhet = sqrt(int2d(Th)(square(uh - nith))); // L^2 norm of the error in nonlinear iteration of het
cout << "        " << k << ": error = " << errnlhet << endl;
if(errnlhet <= 1E-5*normuh) // stopping criteria. If relative error
break;        // is less than 1E-5 stop nonlinear iterations of het
nith = uh; // prepare new iteration
g0=uh*uh/(1.0+uh)^2;
g1=1.0/(1.0+uh)^2;
}
// Check convergence for het
if(k >= noit) {
cout << "Nonlinear solver for het did not converged.\n";
cout << "Error = " << errnlhet << endl;
cout << "|| uh || = " << normuh << endl;
exit(1); // exit the program with error code 1
}

nith= wh; // start nonlinear iteration with last solution for hom
g2=wh*wh/(1.0+wh)^2;
g3=1.0/(1.0+wh)^2;
cout << "    Solving nonlinear homogeneous problem.\n";
for(k=1;k<noit; ++k) {
hom;
normwh = sqrt(int2d(Th)(square(wh))); // L^2 nor of the solution of hom
errnlhom = sqrt(int2d(Th)(square(wh - nith))); // L^2 norm of the error in nonlinear iteration of hom
cout << "        " << k << ": error = " << errnlhom << endl;
if (errnlhom <= 1E-5*normwh) // stopping criteria. If relative error
break;        // is less than 1E-5 stop nonlinear iterations of het
nith = wh; // prepare new iteration
g2=wh*wh/(1.0+wh)^2;
g3=1.0/(1.0+wh)^2;
}
// Check convergence for het
if(k >= noit) {
cout << "Nonlinear solver for hom did not converged.\n";
cout << "Error = " << errnlhom << endl;
cout << "|| wh || = " << normwh << endl;
exit(1); // exit the program with error code 1
}
Vh eh=abs(uh-wh); // Définition de la fonction eh pour laquelle on veut calculer son max en (x,y)
real umaxt = eh[].max; //  max au temps temp
cout << "Umax for time t = " << umaxt << endl;
fouti << t <<" "<< umaxt<< endl;
//umax = max(umax, umaxt);// max jusqu'au temps t
//cout << "Global max = " << umax << endl;

}

//cout << "Final max = " << umax << endl;

savevtk("u_"+i+".vtu",Th,uh,dataname="uh", order=Order);
savevtk("w_"+i+".vtu",Th,wh,dataname="vh", order=Order);

real errH1 = sqrt(int2d(Th)(square(uh - wh)) + int2d(Th)(square(dx(uh)-dx(wh)) + square(dy(uh)-dy(wh))));

cout << " ==== errH1 ( eps = " << eps <<", omega = " << w << ", eps^2*omega = " << eps^2*w << ") = "<< errH1 <<endl;

h(i)=eps^2*w;
error(i)=errH1;
alpha(i)=i>0?1.0/log(2.0)*log(error(i))/log(error(i-1)):1.0;
fout << i <<" "<< h(i) <<" "<< error(i) <<" "<< alpha(i) << endl;
fout.flush;

eps=eps/2;
w=2*w;

}``````

Oui c’est bien ça qu’il faut faire

Bonjour
merci beaucoup pour l’aide. J’ai une dérnière question qui ne concerne pas FF++ directement, je m’excuse. On ne sait jamais peut être que vous pourrez m’aider. Le but de ce calcul du max est de faire une courbe qui montre que uh et wh sont proches l’une de l’autre. Avec ce qu’on a fait, on obtient la courbe ci jointe (le max temps a temps en fonction du temps). Ma question est est ce que cette courbe montre que uh et wh sont proches? Ou bien que fau-il dessiner comme courbe pour voire que les deux sont proches?

Bien cordialement

Ça dépend de ce que vous appelez proches.

Si les résultats sont sensés être égaux, alors l’erreur devrait plutôt être de l’ordre de la précision du solveur (10^-10 - 10^-6).

Si les résultats sont sensés être différents, il faut une base de comparaison pour savoir si il sont effectivement “proches”