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

#1

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();

load "iovtk"

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;

  }
(Simon Garnotel) #2

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

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

load "iovtk"

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;
}
#5

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


load "UMFPACK64"
defaulttoUMFPACK64();

load "iovtk"

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;

  }
(Simon Garnotel) #6

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

#8

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

max_0

(Simon Garnotel) #9

Ç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”