reversible matrix operation

I am debugging my cfd code and I am checking the validity of my matrices.

I have created a finite difference laplacian matrix in cylindrical coordinates, L. I expect if I take the laplacian of a field, and then multiply by the inverse of the laplacian matrix I should be able to return the exact field that I had before. However, for some reason this is not happening. Basically I am performing the following operations.

L*A = B; C = L\B-A;

I expect that C should equal zero, but it doesn't.

I attached some code showing my problem. Can someone explain what is going wrong?

function [pdiff] = ptest()
AR = 1;
nx = 10;
ny = 20;
lx = 1;
ly = 1;
x = linspace(0,lx,nx+1); hx = lx/nx;
y = linspace(0,ly,ny+1); hy = ly/ny;
[X,Y] = meshgrid(y,x);
pcordr = avg(avg(X')');
pcordz = avg(avg(Y')');
rp = ones(ny,1)*avg(y); rp = rp';

pcenter = pcordr.^2.*pcordz.^2;
plong = reshape(pcenter,[],1);

% % Laplacian matrix
Lp = AR*kron(speye(ny),K1(nx,hx,1))+(1/AR)*kron(K1(ny,hy,1),speye(nx))+...
    (1/AR)*kron((1./rp).*Kp1(ny,hy),speye(nx));

% Test for reverse operations
pfield = Lp*plong;
pdiff = Lp\pfield - plong;

function A = K1(n,h,a11)
% a11: Neumann=1, Dirichlet=2, Dirichlet mid=3;
A = spdiags([-1 a11 0;ones(n-2,1)*[-1 2 -1];0 a11 -1],-1:1,n,n)'/h^2;

function A = Kp1(n,h)
% a11: Neumann=1, Dirichlet=2, Dirichlet mid=3;
A = spdiags([-1 1 0;ones(n-2,1)*[-1 0 1];0 -1 1],-1:1,n,n)'/(h*2);

function B = avg(A,k)
if nargin<2, k = 1; end
if size(A,1)==1, A = A'; end
if k<2, B = (A(2:end,:)+A(1:end-1,:))/2; else, B = avg(A,k-1); end
if size(A,2)==1, B = B'; end

Answer 1

[Walter Roberson]

The output that results, C: is it on the order of 1e-16 times the largest (absolute value) input? If so then you are observing round-off error.

[Peter]

no, its considerably larger

[Walter Roberson]

In my experimentation, the largest I have found has been less than the largest of max(eps(A*x)), max(eps(A)), max(eps(x))

Answer 2

[Richard Brown]

Your code doesn't run for me -- first, I don't have a function avg - did you mean to use mean?

Second, if I replace all occurences of avg with mean, the line defining the Laplacian doesn't work - in particular you can't do

(1./rp).*Kp1(ny,hy)

as your function Kp1 defines a matrix and 1./rp is a vector ...

[Peter]

right, sorry. I updated the code

Answer 3

[Richard Brown]

Your Laplacian matrix is not full rank. Because your system is consistent, you have multiple solutions, and you're just getting one of them when you do Lp \ pfield

Unfortunately because your matrix Lp is sparse, you don't see the usual warning message about the matrix being badly scaled or singular ...