lapcomplexeigenvalue.edp

FreeFem++可以生成高质量的有限元网格。可以用于流体力学

//  laplace with matrix  
verbosity=1;
mesh Th=square(20,20,[pi*x,pi*y]);
fespace Vh(Th,P2);
Vh<complex> u1,u2;
Vh ur,ui;
int n=u1.n;
complex[int] Bu1(n),Bu2(n);
complex[int] Au1(n),Au2(n);

complex  sigma = 0; //1.000+1i;                     
complex nu=1+1i;

varf  op(u1,u2)= int2d(Th)(  nu*(dx(u1)*dx(u2) + dy(u1)*dy(u2))  - sigma*u1*u2  )
                    +  on(1,2,3,4,u1=0)
;
varf  a(u1,u2)= int2d(Th)(  nu*(dx(u1)*dx(u2) + dy(u1)*dy(u2))   )
                    +  on(1,2,3,4,u1=0)
;
                   
varf b([u1],[u2]) = int2d(Th)(  u1*u2 ) ;//+  on(1,2,3,4,u1=0);


matrix<complex> OP= op(Vh,Vh,solver=UMFPACK); 
matrix<complex> A= a(Vh,Vh,solver=GMRES); 
matrix<complex> B= b(Vh,Vh,solver=GMRES,eps=1e-20); 
int nev=10;

complex[int] ev(nev); // to store 10 eigen value 
Vh<complex>[int] eV(nev);   // to store 10 eigen vector  

int k=EigenValue(OP,B,sigma=sigma,value=ev,vector=eV,
	         tol=1e-10,maxit=90000,ncv=100);

k=min(k,nev); //  some time the number of converged eigen value 
              // can be greater than nev;
for (int kk=0;kk<k;kk++)
{ 
  int i=kk;	
  u1=eV[i];	
  complex v= ev[i];
  Bu1=B*u1[];
  Au1=A*u1[];
    //  The Rayleigh quotient lambda = x'Ax/x'Bx 
    //  given the eigen value 
  complex xAx = u1[]'*Au1 ;
  complex xBx = u1[]'*Bu1 ;
  //  A u = l * B u
  Bu1 =   v*Bu1;
  u1[] = Au1 -Bu1; 
  if(norm(u1[].sum)>1e-5)
        cout << "BUG :::   zero ==  " <<u1[].sum << endl;
  cout << " ---- " <<  i<< " " <<  v << endl;
  ur=real(eV[i]);
  ui=imag(eV[i]);
//  plot(ur,cmm="Eigen  Vector (real)  "+i+" valeur =" + v  ,wait=1,value=1);
//  plot(ui,cmm="Eigen  Vector (imag)  "+i+" valeur =" + v  ,wait=1,value=1);
}