Prev: Cholesky decomposition
Next: Is there a manual on reading Matlab segmentation fault stack
From: Anubhab on 29 Apr 2010 03:16
Hi,

I need to do a Cholesky decomposition of a matrix, Q. The matrix arises as a weight in computing energy in problems in hydrodynamic stability. Q is negative laplacian in 2d.
Q=-(D^2-k^2); (D is a Chebyshev differentiation matrix along y and k is a wavenumber along x obtained on Fourier transforming).
Energy= y'*Q*y.
I expect Q to be positive definite and turns out all its eigenvalues are >0. I have pasted the code below

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all; clc
N=15; % number of discretization points

% Compute the Chebyshev differentiation matrix
x = cos(pi*(0:N)/N)';
c = [2; ones(N-1,1); 2].*(-1).^(0:N)';
X = repmat(x,1,N+1);
dX = X-X';
D = (c*(1./c)')./(dX+(eye(N+1))); % off-diagonal entries
D = D - diag(sum(D,2)); % diagonal entries
D2 = D^2;

% Impose homogeneous boundary conditions
D2 = D2(2:N,2:N);
I = eye(size(D2)); % Identity matrix

k=1;
Q = -(D2-k^2*I);
all(eig(Q)>0)
chol(Q);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Now when N is small (~10) it works fine. As I increase N I get the following error

??? Error using ==> chol
Matrix must be positive definite.

whereas all the eigenvalues are still positive. Any suggestions why Cholesky decomposition isn't working?
Thanks

anubhab

p.s. The eigenvalues increase in magnitude with N. By N=20 the largest eigenvalues a little less than 10^4.
 | 
Pages: 1
Prev: Cholesky decomposition
Next: Is there a manual on reading Matlab segmentation fault stack