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

9 次查看（过去 30 天）
Max
Max 2024-6-9
评论： Max 2024-6-9
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.

请先登录，再回答此问题。

采纳的回答

David Goodmanson
David Goodmanson 2024-6-9
编辑：David Goodmanson 2024-6-9
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:
  1 个评论
Max
Max 2024-6-9
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!

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Surface and Mesh Plots 的更多信息

标签

产品


版本

R2023b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by