Code works in Octave, but not in Matlab...

I am working on learning how to solve eigenvalue PDEs for a class, and I have written a bit of code that first discretizes the problem, then finds the eigenvalues for the resulting matrix, and then calculates the error in the second eigenvalue using the (known) solution.

When I run my code in Octave, the error estimate is exactly what I would expect (the code plots the error for varying h = [1/4 1/8 1/16 1/32 1/64 hf]; where hf = 0.01, set by the classwork). Also in Octave, the plot of the solution looks just like I would expect.

However when I run the same code in Matlab, the eigenvalue blows up, and I get a rapid oscillation inside of the oscillation that I would expect. Below is the solution for Octave (this is the correct solution.)

The following is the solution that Matlab produces.

This is really, really driving me crazy. I have run the code on several different machines, using several different versions of Matlab (perhaps the same Octave, I'm not sure) and I get the same result.

I have pasted the code below. Much appreciation, and many thanks.

clear 
clf

h = [1/4 1/8 1/16 1/32 1/64 0.01]';
N = 1./h;

xi = 0; 
xf=1;
yi = 0; 
yf = 0;

lamda = -4*pi*pi;
error = zeros(length(h),1);

B = cell(1,length(h));
V = cell(1,length(h));
D = cell(1,length(h));

for k = 1:length(h)

  x = xi+h(k):h(k):xf-h(k);

  A = sparse(N(k)-1,N(k)-1);

  A(1,1) = 2*( 1/((1+x(1))^2) -1/h(k)/h(k));
  A(1,2) = 1/h(k)/h(k) - 1/h(k)/(1+x(1));

  for j=2:N(k)-2
    A(j,j-1) = 1/h(k)/h(k) + 1/h(k)/(x(j)+1);
    A(j,j) = 2*( (1+x(j))^(-2) -1/h(k)/h(k));
    A(j,j+1) = 1/h(k)/h(k) - 1/h(k)/(1+x(j));
  end

  A(end,end-1) = 1/h(k)/h(k) + 1/h(k)/(1+x(end));
  A(end,end) = 2*( (1+x(end))^(-2) -1/h(k)/h(k));

  B{k} = full(A);
  
  [V{k}, D{k}] = eig(B{k});
  
  error(k) = abs(lamda - D{k}(2,2));
end

ErrorDiv4 = error./4;
ErrorDiv4 = [NaN; ErrorDiv4(1:end-1)];
%~ T = table(h,error,ErrorDiv4)
T = [h error ErrorDiv4]

l2 = D{end}(2,2);

C = B{end}-l2*eye(N(end)-1); 
b = zeros(N(end)-1,1);
y = null(C);
y = [yi; y; yf];
x = [xi:h(end):xf];

Q = plot(x,y.')

Answer 1

[John D'Errico]

This looks like you are expecting something that need not happen. Without even looking that deeply at it, the plots are suggestive that you are grabbing a wrong eigenvalue. So a rapidly varying result, rather than a smooth, slowly varying result.

I'd look into the help for eig in Octave, versus MATLAB. There is no presumption that the eigenvalues are returned in any specific order from eig in MATLAB, but you are always taking a specific element from D.

[Zack Fifer]

Yes, in retrospect it was silly to assume that the eigenvalues would be sorted in the order that I expected.

Setting:

EE = eig(B{k});
EE = sort(EE,'descending');
E{k} = EE;

fixed the whole thing.

[John D'Errico]

Ah, good. It was just a guess on my part, but then I've seen many people have problems stemming from the order of eigenvalues as returned from eig.