# How to form the stiffness matrix for the Poisson equation using a spectral method

This is a follow-up question to How do I form the Chebyshev differentiation matrix in MATLAB? The goal of the following code is to solve the Poisson problem:

function [PerInfError]=Poisson2Dn0dc1sev12(N)
%Solve -nabla^2u=f(x,y) on Omega=[-1,1]^2, with u(Gamma)=x^3,  leading
%to I(U) =int(0.5(nabla.u)^2-fu)dxdy
%Exact solution u(x,y) = sin(pix)sin(piy)+x^3 +y.^3
%Current example has f(x,y) = 2pi^2sin(pix)sin(piy)+6x+6y
%Using GLL nodes over 1 spectral element

%------------GLL Nodes, weights and differentiation matrix-----------------
N1=N+1;
cheb=cos(pi*(0:N)/N)';
unif=linspace(-1,1,N1)';
if N<3
x=cheb;
else
x=cheb+sin(pi*unif)./(4*N);
end
P=zeros(N1,N1); %N1xN1 zero matrix
xold=2;

while max(abs(x-xold))>eps
xold=x;
P(:,1)=1;
P(:,2)=x;
%set first collumn entries to 1 (P(:,1) = 1);
%set second collumn entries to x (P(:,2) = x);
% and so on, the collumns of this matrix are the legendre polynomials
for k=2:N
P(:,k+1)=( (2*k-1)*x.*P(:,k)-(k-1)*P(:,k-1) )/k;
%%Bonnets formula ;)
end
% Roots of (1-x^2)L'_N
x=xold-( x.*P(:,N1)-P(:,N) )./( N1*P(:,N1));
end

%---chebyshev differentiation matrix-------------------------------------------------
x=-x;
X=repmat(x,1,N1);%set every collumn of X to x
Xdiff=X-X'+eye(N1); % X MINUS ITS TRANSPOSE YOU IDIOT
L=repmat(P(:,N1),1,N1);
L(1:(N1+1):N1*N1)=1;
D=(L./(Xdiff.*L'));
D(1:(N1+1):N1*N1)=0;
D(1)=-(N1*N)/4;
D(N1*N1)=(N1*N)/4;

w=2./(N*N1*P(:,N1).^2);

%---matrix of GLL x and y coords-------------------------------------------

Xcoord = repmat(x(2:N),1,N-1)';
Ycoord = repmat(x(2:N),1,N-1);

%---defining the known exact solution uEx----------------------------------

UEx=sin(pi.*Xcoord).*sin(pi.*Ycoord);

%----- RHS contribution --- zero dirichlet boundary conditions-------------
b=zeros((N-1)^2,1);
A= zeros((N-1)^2,(N-1)^2);
for k = 1:(N-1)^2;
if mod(k,N-1)==0;
n = N-1;
else
n = mod(k,N-1);
end
m = 1+floor((k-1)./(N-1));

%----Full RHS for non zero dirichlet---------------------------------------

b(k,1)= w(m+1).*w(n+1).*(2*(pi.^2).*sin(pi.*Xcoord(m,n))...
.*sin(pi.*Ycoord(m,n)));

%--- LHS contribution---split into the three contributions-----------------
AY1=zeros(N-1,N-1);
AX1=zeros(N-1,N-1);

for ii = 1:N-1;
AX1(ii,n)=sum(w(n+1).*w(:).*D(:,m+1).*D(:,ii+1)); %for A=A1=A2
AY1(m,ii)=sum(w(m+1).*w(:).*D(:,n+1).*D(:,ii+1)); %for A=A1=A2
end

AXYSum=AX1+AY1;
A(k,:)=reshape(AXYSum,1,(N-1)^2);
end
%---Solving the resulting linear system----equating the Error--------------

UEstV= linsolve(A,b);
UEst=reshape(UEstV,N-1,N-1);
Error=abs(UEst-UEx);
PerInfError=norm(UEx-UEst,inf)./norm(UEx,inf);

end


Here is what I think the LHS should be.

$$\int_{\Omega}\nabla u \nabla v\approx w_l\sum_{i = 0}^{N}\sum_{m = 0}^{N}w_mD_{mi}D_{mk}u_{il} + w_k\sum_{j = 0}^{N}\sum_{n = 0}^{N}w_nD_{nj}D_{nl}u_{kj}$$

Where $D$ is the differentiation matrix mentioned in the original post. I got to this result by letting $u(x,y) = \sum_{i = 0}^{N}\sum_{j = 0}^{N}u_{ij}h_i(x)h_j(y)$ and $v(x,y) = h_k(x)h_l(y)$

Now the matrix, from this, should look something like

The reason for some blocks being diagonal and others being "Full" is that the first term in the above expression only makes sense when $i=k$, this is because it represents the derivative in the x direction and the second term when $j=l$ for the same reason but in the y direction.

I would like to understand how the author has implemented this. I am being quite vain here, assuming that he has done it this way.

• So that is what is called! The Stiffness matrix! – user21030 Jul 26 '16 at 11:45
• What precisely is your question? The loop (over k and ii) seems to exactly compute the LHS you give (look for the two lines with AX1 and AY1). – Christian Clason Jul 26 '16 at 11:53
• I think it was mainly a lack of matlab syntax knowledge that was holding me back. As a side note, the way the author constructs the matrix A is so elegant! Thanks Christian. – user21030 Jul 26 '16 at 15:32