计算科学 - 如何使用谱法形成泊松方程的刚度矩阵 - 吾爱随笔录

这是如何在 MATLAB 中形成切比雪夫微分矩阵的后续问题？以下代码的目标是解决泊松问题：

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

这是我认为 LHS 应该是的。

\int_{Ω} \nabla u \nabla v \approx w_{l} \sum_{i = 0}^{N} \sum_{m = 0}^{N} w_{m} D_{m i} D_{m k} u_{i l} + w_{k} \sum_{j = 0}^{N} \sum_{n = 0}^{N} w_{n} D_{n j} D_{n l} u_{k j}

$\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}$

在哪里 $D$ 是原帖中提到的微分矩阵。我通过让 $u(x,y) = \sum_{i = 0}^{N}\sum_{j = 0}^{N}u_{ij}h_i(x)h_j(y)$ 和 $v(x,y) = h_k(x)h_l(y)$

现在矩阵，从这里，应该看起来像

一些块是对角线而其他块是“完整”的原因是上述表达式中的第一项仅在以下情况下才有意义 $i=k$ , 这是因为它表示 x 方向的导数和第二项当 $j=l$ 出于同样的原因，但在 y 方向。

我想了解作者是如何实现这一点的。我在这里很虚荣，假设他是这样做的。