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 * (CMux);
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.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 * (CMux);
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. 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 * (CMux);
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. surf, surface MATLAB Answers — New Questions
