0
$\begingroup$

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;
%---GLL Quad Weights-------------------------------------------------------

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

enter image description here

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.

$\endgroup$
  • $\begingroup$ So that is what is called! The Stiffness matrix! $\endgroup$ – user21030 Jul 26 '16 at 11:45
  • 1
    $\begingroup$ 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). $\endgroup$ – Christian Clason Jul 26 '16 at 11:53
  • $\begingroup$ 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. $\endgroup$ – user21030 Jul 26 '16 at 15:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy