Consider the P1 finite element matrix, denoted by A, obtained from the problem
\begin{array}{rcl} -\Delta u & =& f \quad\mbox{on } \Omega,\\ \partial_n u&=& 0 \quad \mbox{on } \partial \Omega. \end{array}
We know A is singular. Then consider the P1 finite element matrix, denoted by A0, obtained from the problem
\begin{array}{rcl} -\Delta u & =& f \quad\mbox{on } \Omega,\\ u&=& 0 \quad \mbox{on } \partial \Omega. \end{array}
We know A0 is not singular and could be obtained from A by removing rows and columns associated to boundary dofs.
FreeFem++ generates A0 from A by multiplying by tgv=1e30 diagonal entries of A associated to boundary dofs.
Removing rows and columns for boundary dofs in A is made by the order matrix A00 = A(I^-1,I^-1); where I has -1 at these boundary dofs. We also check that cond(A0)=O(h^-2)
The following code (thanks to Frédéric Hecht) computes the L^2 condition numbers of A and of A0.
download example: CondNumb.edp or return to Matrices and Arrays
verbosity=0; // // Mesh generation // int nsize=3; mesh Th=square(nsize,nsize,[x*pi,y*pi]); real hsize=pi/nsize; fespace Vh(Th,P1); // // identity matrix of n=Vh.ndof size // int n = Vh.ndof; real[int] one(n); one=1.; matrix Id=one; // set(Id,solver=sparsesolver); macro grad(u) [dx(u),dy(u)]// // // To generate finite element matrix with tgv at boundary dofs // varf a0(u,v) = int2d(Th)(grad(u)'*grad(v)) + on(1,2,3,4,u=0); // // To generate finite element matrix with no tgv at boundary dofs // It amounts to a singular matrix because is the fem of // // -Lap(u)=f con \partial_n u=0 on boundary // varf a(u,v) = int2d(Th)(grad(u)'*grad(v)) ; // // to generate a vector with 1 on boundary dofs // varf von(u,v)= on(1,2,3,4,u=1); real[int] on1=von(0,Vh, tgv=1); // one on boundary dofs int[int] I(n); // to get indices of boundary dofs for (int k=0,i=0;i<n;++i) if( on1[i]) I[i]=-1; // set -1 at indices on boundary dofs else I[i]=k++; //numbers all remaining dofs starting at 0 // // get matrix for problem with tgv at boundary dofs // matrix A0 = a0(Vh,Vh,solver=sparsesolver); // // get matrix for singular problem (no tgv at boundary dofs) // matrix A = a(Vh,Vh,solver=sparsesolver); // // Computational trick to get the cond(A) in L2 norm using the // eigenvalues // int nev=3; // the number of eigenvalues searched real[int] evmin(nev),evmax(nev); // // closest eig values to zero for A. Saved in evmin[0] // int k0=EigenValue(A,Id,sym=true,value=evmin); // // closest eig values to zero for A^-1. Saved in evmax[0] // As a consequence: // // 1/evmax[0] is the largest eig values for A // int k=EigenValue(Id,A,sym=true,value=evmax); cout << " # of Converged eig values for A is = "<<k0 << " # of Converged eig values for A^{-1} is = " << k << endl; cout << 1/evmin[0] << " " << evmax[0] << endl; cout << " Pure Neumann problem (singular) " <<endl; // // L^2 condition number= (largest eig value)/(min eig value)=(1/evmax[0]) / evmin[0] // cout << " cond number of A = " << 1/(evmax[0] * evmin[0]) << endl; // // Remove from A the rows and columns for boundary dofs // i.e., with tgv=1 on diagonal entry // cout<<"---------------------------------------------------------------------- "<<endl; cout<<"Vector I sets -1 at boundary dofs and numbers the rest starting from 0 "<<endl; cout<<"---------------------------------------------------------------------- "<<endl; cout<<I<<endl; matrix A00 = A(I^-1,I^-1); // // Remove from A0 the rows and columns for boundary dofs, // i.e., with tgv=1e30 on diagonal entry // It serves to check result // matrix A01 = A0(I^-1,I^-1); matrix Id0 = Id(I^-1,I^-1); // // set solvers // set(A00,solver=sparsesolver); set(A01,solver=sparsesolver); set(Id0,solver=sparsesolver); //same trick as before for A00 k0=EigenValue(A00,Id0,sym=true,value=evmin); k =EigenValue(Id0,A00,sym=true,value=evmax); cout << " # of Converged eig values for A00 is = "<<k0 << " # of Converged eig values for A00^{-1} is = " << k << endl; cout << k0 << " " << k << endl; cout << 1/evmin[0] << " " << evmax[0] << endl; cout << " Homogeneous Dirichlet problem (no singular) " <<endl; cout << " Matrix deduced from that of pure Neumann problem " <<endl; real trucond= 1/(evmax[0] * evmin[0]); cout << " True cond number = " << trucond << endl; cout << " h = " <<hsize << endl; cout<<" cond= O(h^-2) ==> cond*h^2 approx constant "<<endl; cout << " True cond number*h^2 = " << trucond*hsize^2 << endl; //same trick as before for A01 k0=EigenValue(A01,Id0,sym=true,value=evmin); k =EigenValue(Id0,A01,sym=true,value=evmax); cout << "# of Converged eig values for A01 is = "<<k0 << "# of Converged eig values for A01^{-1} is = " << k << endl; cout << k0 << " " << k << endl; cout << 1/evmin[0] << " " << evmax[0] << endl; cout << " Homogeneous Dirichlet problem (no singular) " <<endl; cout << " Matrix generated from that of Dirichlet problem " <<endl; cout << " True cond number = " << 1/(evmax[0] * evmin[0]) << endl;
return to Matrices and Arrays