How do I get rid of these weird lines in the surf plot?

Hello! I was coding something in Matlab (R2023b) and my surf plots all seem to look nice but then they suddenly changed but I didnt even change anything in the code. Even scripts that I didn't touch at all for day have that... I then updated to R2024a but it didnt help at all... Here is one of the codes:

%SymPoisson.m
close all;
clear varibles;
tic

% Inverse Multiquadric RBF and Laplacians
rbf = @(e,r) 1./sqrt(1+(e*r).^2); ep = 1;
Lrbf = @(e,r) e^2*((e*r).^2-2)./(1+(e*r).^2).^(5/2);
L2rbf = @(e,r) 3*e^4*(3*(e*r).^4-24*(e*r).^2+8)./(1+(e*r).^2).^(9/2);

% Datafunction and Laplacian
u = @(x) sin(pi*x(:,1)).*cos(pi*x(:,2)/2);
Lu = @(x) -1.25*pi^2*sin(pi*x(:,1)).*cos(pi*x(:,2)/2);

% Evaluation Points in Rectangle [a,b]x[c,d]
Neval = 100;
a=0;b=4;c=0;d=4;
[~,~,xeval]=GetRectGrid(a,b,c,d,Neval); % xeval is  (Neval^2)x2 Matrix with each line beeing a Point in R^2
% Reshape for Plot
X = reshape(xeval(:,1),Neval,Neval);
Y = reshape(xeval(:,2),Neval,Neval);

% Initialize Arrays with Errors and Conditionnumbers
m=2;M=50; % Bounds for Iteration
rms_fehler = zeros(size(m:M));
max_fehler = zeros(size(m:M));
avg_fehler = zeros(size(m:M));
cond2 = zeros(size(m:M));
anzahlen  = zeros(size(m:M));

% Evaluate Datafunction on Rectangle
ueval = u(xeval);



k=1; %Counter
for N = m:M
disp(['Anzahl Stützstellen: ', num2str(N^2)])

% Get Collocationpoints on Boundary and Interior
[xint,xbdy,x]=GetRectGrid(a,b,c,d,N);
NI = size(xint,1);
NB = size(xbdy,1);

% Evaluation-Matrix
DM_int = DistMatr(xeval,xint);
LEM = Lrbf(ep,DM_int);
DM_bdy = DistMatr(xeval,xbdy);
BEM = rbf(ep,DM_bdy);
EM = [LEM BEM];

% Collocation-Matrix
DM_II = DistMatr(xint,xint);
LLCM = L2rbf(ep,DM_II);
DM_IB = DistMatr(xint,xbdy);
LBCM = Lrbf(ep,DM_IB);
BLCM = LBCM';
DM_BB = DistMatr(xbdy,xbdy);
BBCM = rbf(ep,DM_BB);
CM = [LLCM LBCM; BLCM BBCM];

% Evaluate Datafunction on Collocationpoints
ux = zeros(N^2,1);
ux(1:NI) = Lu(xint);
% Boundary-Conditions
indx = find(xbdy(:,2)==0 | xbdy(:,2)==4);
ux(NI+indx) = sin(pi*xbdy(indx,1));


% Evaluate the Approx. function
warning('off','MATLAB:nearlySingularMatrix')
approx = EM * (CM\ux);
warning('on','MATLAB:nearlySingularMatrix')

% Collect errors and conditionnumbers
rms_fehler(k) = norm(approx-ueval)/sqrt(Neval);
max_fehler(k) = max(abs(approx-ueval));
avg_fehler(k) = mean(abs(approx-ueval));
cond2(k) = cond(CM);
anzahlen(k) = N^2;
k=k+1;

% Plot of the last Approx. function and Datafunction and Error-Function
if N == M 
   
    % These Plots are WEIRD
    figure(N)
    approx2 = reshape(approx,Neval,Neval);
    surf(X,Y,approx2,'FaceAlpha',0.5,'EdgeColor','interp');
    colormap cool; camlight; lighting gouraud; material dull;
    title(['Approximation to the Solution of the Poisson-Equation for N=',num2str(N^2)])
    
    figure(1);
    ueval2 = reshape(ueval,Neval,Neval);
    surf(X,Y,ueval2,'FaceAlpha',0.5,'EdgeColor','interp')
    colormap cool; camlight; lighting gouraud; material dull;
    title('Solution of the Poisson-Equation')

    figure(N+1);
    fehler = abs(approx2-ueval2);
    surf(X,Y,fehler,'FaceAlpha',0.5,'EdgeColor','interp');
    colormap cool; camlight; lighting gouraud; material dull;
    title(['Errorfunction for N=',num2str(N^2)])
    
end

end


% Plot Error-Decay
figure(2)
semilogy(anzahlen,rms_fehler','-r')
hold on;
semilogy(anzahlen,max_fehler','-b')
hold on;
semilogy(anzahlen,avg_fehler','-g')
hold off;
title('Error-Decay')
legend('RMS-Error','Max-Error','Average Error')
xlabel('Number of Collocationpoints')
ylabel('Error')

% Plot Collocation-Points
figure(4);
scatter(xint(:,1),xint(:,2),'ob');
hold on;
scatter(xbdy(:,1),xbdy(:,2),'or');
hold off;
title(['Collocationpoints N=',num2str(M^2)])
legend('Interior-Points', 'Boundary-Points')
grid on;

% Plot Conditionnumbers 
figure(5)
semilogy(anzahlen,cond2)
title('Cond2 of Collocation-Matrix')
xlabel('Number of Collocation-Points')
ylabel('Cond2')

toc

function [xint,xbdy,x] = GetRectGrid(a,b,c,d,N)

% Assertions
if(a>=b)
 error('a muss echt kleiner als b sein!');
end
if(c>=b)
 error('c muss echt kleiner als d sein!');
end
if(N<=1)
    error('N muss echt größer als 1 sein!');
end

% Get the Grid
[X1,X2] = meshgrid(linspace(a,b,N),linspace(c,d,N));
x = [X1(:),X2(:)];

% Sort for inner points and boundary points

k = 1;l=1; % Zähler
xbdy = zeros(4*N-4,2); % Hier kommen die Randpunkte rein
xint = zeros(N^2-(4*N-4),2); % und hier die inneren Punkte

for j = 1:N^2
    if x(j,1) == a || x(j,1) == b || x(j,2) == c || x(j,2) == d
        xbdy(k,:) = x(j,:);
        k = k+1;
    else
        xint(l,:) = x(j,:);
        l = l+1;
    end
end

x = [xint;xbdy];

end

function D = DistMatr(A,B)

A_width = width(A);
B_width = width(B);

if A_width~=B_width
    error('Die Matrizen sind nicht kompatibel.');
end

% Evaluate DistanceMatrix
A2 = sum(A.^2,2);
B2 = sum(B.^2,2);
D = sqrt(max(bsxfun(@plus,A2,B2') - 2*A*B',0));

end

Now look at the Surf-Plots:

I really don't know why it does that... It looks like it connects Points with lines that shouldn't be connected... Weird is, that the follow code works:

%NewtonBasisPlot.m
tic
close all;
clear variables;

% Gaussian RBF
K = @(epsilon,x) exp(-(epsilon*x).^2); epsilon = 1;

n = 25; 
N = 50;

% Get Evaluation-Grid
Omega = GetPoints('grid',N,0,1);

% Get Support-Points by P-Greedy Algorithm
X=GetPoints('P-Greedy',n,0,1,K,epsilon);

% Choose wich Newtonbasisfunction to Plot
NewtonIndex = [7 25]; %Bsp: N_7(x) und N_25(x)

% Kernel-Matrices
A_X = K(epsilon,DistMatr(X,X));
A_Omega = K(epsilon,DistMatr(Omega,X));

% Newton Basis via Cholesky
L = chol(A_X,'lower');
newton = A_Omega/L';

%  Plot Newtonbasisfunctions
x=reshape(Omega(:,1),N,N);
y=reshape(Omega(:,2),N,N);
Nmat = cellfun(@(nvec) reshape(nvec,N,N), ...
    num2cell(newton,1),'UniformOutput',0);
h_surf = zeros(length(NewtonIndex));
for k=1:length(NewtonIndex)
    h_surf(k) = figure;
    surf(x,y,Nmat{NewtonIndex(k)},'FaceAlpha',0.5,'EdgeColor','interp')
    colormap cool; camlight; lighting gouraud; material dull;
    hold on;
end

% Plot Support-Points
figure('Name','Points')
scatter(X(:,1),X(:,2),'xb','LineWidth',1);
grid on;
toc


function X = GetPoints(type,n,a,b,K,epsilon)

assert(a<b)

if (strcmp(type,'grid'))

    [X1,X2] = meshgrid(linspace(a,b,n));
    X = [X1(:),X2(:)];
    
% Via Powerfunction
elseif(strcmp(type,'P-Greedy'))

    % Number of Points to choose from
    N = 50;
    assert(n<N*N)
    
    Omega = GetPoints('grid',N,a,b);
    
    % Start with "random" Point
    X = Omega(floor((N*N)/2),:);

    for j=1:n

        % Kernel-Matrices
        A_X = K(epsilon,DistMatr(X,X));
        A_Omega = K(epsilon,DistMatr(Omega,X));
        
        % Evaluate Powerfunction via P^2(x) = K(x,x)-T(x)*A_{K,X}*T(x)^T,
        Kxx = K(epsilon,0)*ones(N*N,1); %K(x,x)=phi(0) for RBF-Kernel
        warning('off','MATLAB:nearlySingularMatrix')
        B = A_Omega / A_X; %T(x)*A_{K,X}
        C = (B*A_Omega'); %T(x)*A_{K,X}*T(x)^T
        warning('on','MATLAB:nearlySingularMatrix')
        Psquare = Kxx - diag(C); %P^2(x)

        % Add Point if its max of powerfunction
        [~,index]=max(abs(Psquare));
        X=[X;Omega(index(1),:)];  
        
        % Console
        disp(['P-Greedy Iteration: ', num2str(j)])
    end
    X=X(1:end-1,:);
    
else

    Error('Punkteauswahl-Typ nicht bekannt!');
  
end

end

And its giving me the usual and wanted surface Plot of the Functions:

Maybe someone of you guys has a clue whats going on here... Thank you very much! Kind Regards, Max.

Answer 1

[David Goodmanson]

Hello Max,

It appears that the GetRectGrid function creates points whose order is imcompatible with the reshape functions that follow on the next two lines and create X and Y. If you immediately follow

[~,~,xeval]=GetRectGrid(a,b,c,d,Neval);
with
xeval = sortrows(xeval); 

then it works to use reshape on the new xeval to obtain X and Y. I thought it might be necessary after the command

[xint,xbdy,x]=GetRectGrid(a,b,c,d,N)

to do sortrows on those three as well, but it seems not to be the case. Here is a sample plot after using sortrows:

[Max]

Thank you very much for your Answer! It works now! Apparently I shouldn't do x=[xint,xbdy] at the end of the function GetRectGrid. I did that because I needed the points to be sortet like that in some other script where I want to use this function. I wasn't aware that this isn't working like that.

Thank you so much!