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

%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

%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