1// Computation of the eigen value and eigen vector of the 2// Dirichlet problem on square $]0,\pi[^2$ 3// Functionnal interface 4// ---------------------------------------- 5// we use the inverse shift mode 6// the shift is given with sigma real 7// ------------------------------------- 8// find $\lamda$ such that: 9// $$ \int_{\omega} \nabla u_ \nabla v = \lamba \int_{\omega} u \nabla v $$ 10verbosity=1; 11mesh Th=square(20,20,[pi*y,pi*x]); 12fespace Vh(Th,P2); 13Vh u1,u2; 14int n=Vh.ndof; 15 16real sigma = 00; // value of the shift 17 18varf a(u1,u2)= int2d(Th)( dx(u1)*dx(u2) + dy(u1)*dy(u2) - sigma* u1*u2 ) 19 + int1d(Th)(u1*u2) ; // Boundary condition 20 21varf b([u1],[u2]) = int2d(Th)( u1*u2 ) ; // no Boundary condition 22 23matrix A= a(Vh,Vh,solver=sparsesolver); 24matrix B= b(Vh,Vh,solver=sparsesolver); 25 26func real[int] FA1(real[int] & u) { real[int] Au=A^-1*u;return Au;} 27func real[int] FB(real[int] & u) { real[int] Au=B*u;return Au;} 28func real[int] FA(real[int] & u) { real[int] Au=A*u;return Au;} 29func real[int] FB1(real[int] & u) { real[int] Au=B^-1*u;return Au;} 30// important remark: 31// the boundary condition is make with exact penalisation: 32// we put 1e30=tgv on the diagonal term of the lock degre of freedom. 33// So take dirichlet boundary condition just on $a$ variationnal form 34// and not on $b$ variationnanl form. 35// because we solve 36// $$ w=A^-1*B*v $$ 37 38int nev=20; // number of computed eigen valeu close to sigma 39 40real[int] ev(nev); // to store nev eigein value 41Vh[int] eV(nev); // to store nev eigen vector 42 43 44int k=EigenValue(n,A1=FA1,B=FB,A=FA,B1=FB1,sym=true,sigma=sigma,value=ev,vector=eV,tol=1e-10,maxit=1000,ncv=200 45 ,mode=3,which="LM"); 46// tol= the tolerace 47// maxit= the maximal iteration see arpack doc. 48// ncv see arpack doc. 49// the return value is number of converged eigen value. 50k=min(k,nev); // some time the number of converged eigen value 51 // can be greater than nev; 52int nerr=0; 53for (int i=0;i<k;i++) 54{ 55 u1=eV[i]; 56 real gg = int2d(Th)(dx(u1)*dx(u1) + dy(u1)*dy(u1)) + int1d(Th)(u1*u1); 57 real mm= int2d(Th)(u1*u1) ; 58 real err = int2d(Th)(dx(u1)*dx(u1) + dy(u1)*dy(u1) - (ev[i])*u1*u1) + int1d(Th)(u1*u1) ; 59 if(abs(err) > 1e-6) nerr++; 60 cout << " ---- " << i<< " " << ev[i] << " err= " << err << " --- "<<endl; 61 62 // FFCS: add 3D view capabilities 63 plot(eV[i],cmm="Eigen Vector "+i+" valeur =" + ev[i] ,wait=0,value=1,dim=3,fill=1,CutPlane=0,ShowAxes=0); 64} 65 66