I am working on a 4th order 2D Inverse Problem using the matlab OO interface in getfem++ . As a result, I calculate the bilaplacian, but I keep getting a badly scaled matrix with a close to 0 condition number. Below is the code that produces this. Is there a way to correct this.
N = 20;
h = 1/(N+1);
m = gf_mesh('cartesian', 0:h:1, 0:h:1);
mf = gf_mesh_fem(m,1);
mfd = gf_mesh_fem(m,1);
[xs, ys] = meshgrid(0:1*h:1, 0:1*h:1);
gf_mesh_fem_set(mf, 'fem', gf_fem('FEM_QK(2,4)'));
gf_mesh_fem_set(mfd, 'fem', gf_fem('FEM_QK(2,4)'));
mim = gf_mesh_im(m, gf_integ('IM_GAUSS_PARALLELEPIPED(2,10)'));
% Boundary conditions:
border = gf_mesh_get(m, 'outer_faces');
gf_mesh_set(m, 'boundary', 1, border);
% Degrees of freedom
dof = gf_mesh_fem_get(mf, 'nbdof');
dofa = gf_mesh_fem_get(mfd, 'nbdof');
u = '256.*(x.^2).*(y.^2).*((1-x).^2).*(1-y).^2';
Us = gf_mesh_fem_get(mfd, 'eval', {u});
a = '1+0.25.*sin(pi.*x).*sin(pi.*y)';
avec = gf_mesh_fem_get(mfd, 'eval', {a});
a0 = gf_mesh_fem_get(mfd, 'eval', {1.5})';
K = @(a) gf_asm('bilaplacian', mim, mf, mfd, a);
K1 = gf_asm('laplacian', mim, mfd, mfd, ones(dofa,1));
L = @(u) gf_asm('volumic',...
['u=data(#2);', ...
't=comp(Hess(#2).Hess(#2).Base(#1))(:,i,i,:,i,i,:);',...
'M(#2,#1)+=t(j,:,:).u(j)'],...
mim, ...
mfd, ...
mf, ...
u);
M = gf_asm('mass matrix', mim, mf);
Ma = gf_asm('mass matrix', mim, mfd);