Approximation of Solution to Poisson-Equation via RBFs but with a change of Basis

Hello! I am currently trying to implement some idea into Matlab but I can’t do it somehow. Probably because my brain is not big enough or because I don’t have that much experience in Matlab… This might be a long text but maybe I will find someone that is interested in this topic aswell and could help me…
First the Theory behind all that:
Lets say we have an open and bounded Domain and we want to solve the Equation

where is known and we have Dirichlet-Boundarycondition

for some known . One approach to do so is to choose

where holds some collocationpoints ( in the Interior and on the Boundary of ) and is an RBF-Kernel, i.e. for some , e.g. is the Gaussian-RBF. The approximand then should meet the conditions
1) ,
2)
wich leads us to solve the linear system of equations

to get the Coefficients where
1) ,
2) ,
3) ,
4) .
Now this one I was already able to implement for the example

with the analytical solution

where I used the inverse Multiquadric
.
This is the corresponding script:
% IMQ-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);

% Generate Evaluation Points
Neval = 100; % Sqrt of Number of Evaluationpoints
a=0;b=4;c=a;d=b; % Bounds for rectangle
[~,~,xeval]=GetRectGrid(a,b,c,d,Neval);
% Reshape for Plot
X = reshape(xeval(:,1),Neval,Neval);
Y = reshape(xeval(:,2),Neval,Neval);

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

% Generate Collocatoinpoints
N = 25; % Sqrt of number of collocationpoints
[xint,xbdy,x]=GetRectGrid(a,b,c,d,N);
NI = size(xint,1);
NB = size(xbdy,1);

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

% Build Collocationmatrix for the system of equations
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 Collocatoinpoints
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 Approximand
warning(‘off’,’MATLAB:nearlySingularMatrix’)
approx = EM * (CMux);
warning(‘on’,’MATLAB:nearlySingularMatrix’)

% Plot Approximand and analytical solution
figure(1)
approx2 = reshape(approx,Neval,Neval);
surf(X,Y,approx2,’FaceAlpha’,0.5,’EdgeColor’,’interp’);
colormap cool; camlight; lighting gouraud; material dull;
title([‘Approximation of solution to Poisson-Equation for N=’,num2str(N^2)])

function D = DistMatr(A,B)
% Helping function for Evaluating RBFs.

% See if matrices are compatible
A_width = width(A);
B_width = width(B);
if A_width~=B_width
error(‘The matrices are not compatible.’);
end

%Berechnung der Abstandsmatrix
A2 = sum(A.^2,2);
B2 = sum(B.^2,2);
D = sqrt(max(bsxfun(@plus,A2,B2′) – 2*A*B’,0));

end

function [xint,xbdy,x] = GetRectGrid(a,b,c,d,N)
% Generated Grid on choosen rectangle and sorts for interior and boundary
% points.

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

% Sort for interior and boundary points
k = 1;l=1; % Counter (There are probably better solutions to this but idc)
xbdy = zeros(4*N-4,2); % Boundary Points
xint = zeros(N^2-(4*N-4),2); % Interior Points

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
end
Okay so now I want to do basically the same but first I want to do a change of basis. Currently we are working with the Standardbasis that spans some finite dimensional functionspace. The reason to change the basis is that it appears that you gain better evaluation stability if you want to achieve higher accuracy by e.g. choosing very small values for . There’s a lot of theory behind that but lets just lay that by the side for now.
The Idea is to make a simple change of basis to gain some new basis that spans the same space. Therefore we chose
.
Note that we now have instead of for simplicity. We can rewrite this approach to

where is the Coefficientmatrix that holds the Coefficients . Now theres a Theorem that I will not proof here but I will refer to this work by Maryam Pazouki and Robert Schaback in Theorem 3.1 and Theorem 6.1 but the Theorem basically says that every possible basis we can get by this is uniquley characterizable by a factorization of the Matrix

into some Matrix and the inverse of the Coefficientmatrix , i.e.. . For example if we take the Cholesky-Decomposition of we end up with as the coefficientmatrix wich leads us to the so called Newton-Basis , described here in chapter 3.7. Since the RBF-Kernel is because we take the double Laplacian of it, the Newtonbasisfunktions are too since they are just linear combinations of the Kernel and the laplacian is a linear differential operator. Therefore
and .
In summary I want do the approximation of the solution of the Poisson-Equation but not use the standard rbf basis but do a change of basis into the newton basis and then chose the approach

for the approximand to solve the linear system of equations that I describe above but now it has Submatrices and . Obviously after that I’d like to Plot the resulting approximand to see how it looks. I already tried it, but it didnt work like that… I just have alot of trouble to implement that… Here is my Approach:
% Inverse Multiquadric RBF and its Laplacians
rbf = @(e,r) 1./sqrt(1+(e*r).^2); ep = 2;
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 its 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);

% Evaluationpoints
Neval = 100;
a=0;b=4;c=a;d=b;
[~,~,xeval]=GetRectGrid(a,b,c,d,Neval);
X = reshape(xeval(:,1),Neval,Neval);
Y = reshape(xeval(:,2),Neval,Neval);

% Collocationpoints
N = 20;
[xint,xbdy,x]=GetRectGrid(a,b,c,d,N);
x=[xint;xbdy]; % Sorting x… hmm…
NI = size(xint,1);
NB = size(xbdy,1);

% 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));
% .. and on the Evaluationpoints
ueval = u(xeval);

% Evaluate needed Kernel-Matrices
% % % This part works fine. B is looking good.
A1 = rbf(ep,DistMatr(x,x));
A2 = rbf(ep,DistMatr(xeval,x));
% Cholesky decomposition
L = chol(A1,’lower’);
% Each column is one newtonbasisfunction evaluated on the evaluationpoints
newton = A2/L’;
% This is the Coefficent Matrix
B = inv(L)’;

% Collocationmatrix
% % % This is probably not the way to do it… I dont know how to use B…
DM_II = DistMatr(xint,xint);
LLCM = L2rbf(ep,DM_II)*B(1:NI,1:NI);
DM_IB = DistMatr(xint,xbdy);
LBCM = Lrbf(ep,DM_IB).*B(1:NI,NI+1:end);
DM_BI = DistMatr(xbdy,xint);
BLCM = Lrbf(ep,DM_BI).*B(NI+1:end,1:NI);
DM_BB = DistMatr(xbdy,xbdy)*B(NI+1:end,NI+1:end);
BBCM = rbf(ep,DM_BB);
CM = [LLCM LBCM; BLCM BBCM];

% Evaluate Approximand
warning(‘off’,’MATLAB:nearlySingularMatrix’)
approx = newton * (CMux);
warning(‘on’,’MATLAB:nearlySingularMatrix’)

figure(1) % Here you can see how its just not working how I want it to.
approx2 = reshape(approx,Neval,Neval);
surf(X,Y,approx2,’FaceAlpha’,0.5,’EdgeColor’,’interp’);
colormap cool; camlight; lighting gouraud; material dull;
title([‘Approximation of the solution to the Poisson-Equation for N=’,num2str(N^2)])

figure(2);
ueval2 = reshape(ueval,Neval,Neval);
surf(X,Y,ueval2,’FaceAlpha’,0.5,’EdgeColor’,’interp’)
colormap cool; camlight; lighting gouraud; material dull;
title(‘Solution to the Poisson-Equation’)
I hope someone of you could help me manage this problem. Thank you very much! Kind regards, Max.Hello! I am currently trying to implement some idea into Matlab but I can’t do it somehow. Probably because my brain is not big enough or because I don’t have that much experience in Matlab… This might be a long text but maybe I will find someone that is interested in this topic aswell and could help me…
First the Theory behind all that:
Lets say we have an open and bounded Domain and we want to solve the Equation

where is known and we have Dirichlet-Boundarycondition

for some known . One approach to do so is to choose

to get the Coefficients where
1) ,
2) ,
3) ,
4) .
Now this one I was already able to implement for the example

with the analytical solution

% 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);

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

% Generate Collocatoinpoints
N = 25; % Sqrt of number of collocationpoints
[xint,xbdy,x]=GetRectGrid(a,b,c,d,N);
NI = size(xint,1);
NB = size(xbdy,1);

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

% Evaluate Datafunction on Collocatoinpoints
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 Approximand
warning(‘off’,’MATLAB:nearlySingularMatrix’)
approx = EM * (CMux);
warning(‘on’,’MATLAB:nearlySingularMatrix’)

function D = DistMatr(A,B)
% Helping function for Evaluating RBFs.

% See if matrices are compatible
A_width = width(A);
B_width = width(B);
if A_width~=B_width
error(‘The matrices are not compatible.’);
end

%Berechnung der Abstandsmatrix
A2 = sum(A.^2,2);
B2 = sum(B.^2,2);
D = sqrt(max(bsxfun(@plus,A2,B2′) – 2*A*B’,0));

end

function [xint,xbdy,x] = GetRectGrid(a,b,c,d,N)
% Generated Grid on choosen rectangle and sorts for interior and boundary
% points.

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

% Datafunction and its 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);

% Evaluationpoints
Neval = 100;
a=0;b=4;c=a;d=b;
[~,~,xeval]=GetRectGrid(a,b,c,d,Neval);
X = reshape(xeval(:,1),Neval,Neval);
Y = reshape(xeval(:,2),Neval,Neval);

% Collocationpoints
N = 20;
[xint,xbdy,x]=GetRectGrid(a,b,c,d,N);
x=[xint;xbdy]; % Sorting x… hmm…
NI = size(xint,1);
NB = size(xbdy,1);

% Evaluate Approximand
warning(‘off’,’MATLAB:nearlySingularMatrix’)
approx = newton * (CMux);
warning(‘on’,’MATLAB:nearlySingularMatrix’)

figure(2);
ueval2 = reshape(ueval,Neval,Neval);
surf(X,Y,ueval2,’FaceAlpha’,0.5,’EdgeColor’,’interp’)
colormap cool; camlight; lighting gouraud; material dull;
title(‘Solution to the Poisson-Equation’)
I hope someone of you could help me manage this problem. Thank you very much! Kind regards, Max. Hello! I am currently trying to implement some idea into Matlab but I can’t do it somehow. Probably because my brain is not big enough or because I don’t have that much experience in Matlab… This might be a long text but maybe I will find someone that is interested in this topic aswell and could help me…
First the Theory behind all that:
Lets say we have an open and bounded Domain and we want to solve the Equation

where is known and we have Dirichlet-Boundarycondition

for some known . One approach to do so is to choose

to get the Coefficients where
1) ,
2) ,
3) ,
4) .
Now this one I was already able to implement for the example

with the analytical solution

% 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);

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

% Generate Collocatoinpoints
N = 25; % Sqrt of number of collocationpoints
[xint,xbdy,x]=GetRectGrid(a,b,c,d,N);
NI = size(xint,1);
NB = size(xbdy,1);

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

% Evaluate Datafunction on Collocatoinpoints
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 Approximand
warning(‘off’,’MATLAB:nearlySingularMatrix’)
approx = EM * (CMux);
warning(‘on’,’MATLAB:nearlySingularMatrix’)

function D = DistMatr(A,B)
% Helping function for Evaluating RBFs.

% See if matrices are compatible
A_width = width(A);
B_width = width(B);
if A_width~=B_width
error(‘The matrices are not compatible.’);
end

%Berechnung der Abstandsmatrix
A2 = sum(A.^2,2);
B2 = sum(B.^2,2);
D = sqrt(max(bsxfun(@plus,A2,B2′) – 2*A*B’,0));

end

function [xint,xbdy,x] = GetRectGrid(a,b,c,d,N)
% Generated Grid on choosen rectangle and sorts for interior and boundary
% points.

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

% Datafunction and its 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);

% Evaluationpoints
Neval = 100;
a=0;b=4;c=a;d=b;
[~,~,xeval]=GetRectGrid(a,b,c,d,Neval);
X = reshape(xeval(:,1),Neval,Neval);
Y = reshape(xeval(:,2),Neval,Neval);

% Collocationpoints
N = 20;
[xint,xbdy,x]=GetRectGrid(a,b,c,d,N);
x=[xint;xbdy]; % Sorting x… hmm…
NI = size(xint,1);
NB = size(xbdy,1);

% Evaluate Approximand
warning(‘off’,’MATLAB:nearlySingularMatrix’)
approx = newton * (CMux);
warning(‘on’,’MATLAB:nearlySingularMatrix’)

Cart

Cart