Legendre 谱方法及其算例

谱方法按照正交多项式的不同可分为 Legendre 谱方法，Chebshev 谱方法，Fourier 谱方法，Laguerre 谱方法，Hermite 谱方法。这里介绍一维问题的 Legendre 谱方法，然后 MATLAB 编程算出数值解以及误差。

两点边值问题

考虑两点边值问题：

$$\left\{\begin{aligned} &-u''(x)+\alpha u(x)=f(x),\quad x \in I=(-1,1), \\ &u(-1)=0,\ u(1)=0. \end{aligned} \right.$$

Legendre-Galerkin 方法

弱形式：

$$\left\{\begin{aligned} &\text{Find } u \in H_{0}^{1}(I) \text{ such that} \\ &(u',v')+\alpha(u,v)=(f,v),\quad v \in H^{1}_{0}(I). \end{aligned} \right.$$

令 满足边界条件，[Shen2011, P146] 可得 , ，所以

$$\phi_{k}(x)=L_{k}(x) - L_{k+2}(x).$$

设 是 的一组基函数，我们将逼近解展开为

$$u_{N}(x)=\sum_{k=0}^{N-2} \hat{u}_{k} \phi_{k}(x) \in X_{N}.$$

谱格式：

$$\left\{\begin{aligned} &\text{Find } u_N \in X_N \text{ such that} \\ &(u_N',v_N')+(u_N,v_N)=(f,v_N),\quad v_N \in X_N. \end{aligned} \right.$$

将 的展开式代入谱格式方程，取 。令

$$\begin{aligned} &f_{k}=\int_{I} f_{N} \phi_{k} d x, \quad \boldsymbol{f}=\left(f_{0}, f_{1}, \ldots, f_{N-2}\right)^{T}, \\ &u_{N}=\sum_{j=0}^{N-2} \hat{u}_{j} \phi_{j}, \quad \boldsymbol{u}=\left(\hat{u}_{0}, \hat{u}_{1}, \ldots, \hat{u}_{N-2}\right)^{T}, \\ &s_{k j}=-\int_{I} \phi_{j}^{\prime \prime} \phi_{k} d x, \quad m_{k j}=\int_{I} \phi_{j} \phi_{k} d x, \\ \end{aligned}$$ $$S=\left(s_{k j}\right)_{0 \leq k, j \leq N-2}, \quad M=\left(m_{k j}\right)_{0 \leq k, j \leq N-2}.$$

可得以下方程组

$$(S+\alpha M) \boldsymbol{u}=\boldsymbol{f}.$$

刚度矩阵 是一个对角矩阵 [P146-4.22]，它的非零元素为

$$s_{kk}=-(4 k+6)b_{k}=4k+6,\quad k=0,1,\cdots,N-2.$$

质量矩阵 是一个对称的五对角矩阵 [P146-4.23]，它的非零元素为

$$m_{kj}=\left\{ \begin{aligned} &\frac{2}{2k+1}+\frac{2}{(2k+5)} , && j=k, \\ &-\frac{2}{(2k+5)}, && j=k+2, \\ &-\frac{2}{(2k+1)}, && j=k-2. \end{aligned}\right.$$

如果想让生成的刚度矩阵 为单位对角阵，即 (相应的质量矩阵也会改变)， 可选择基函数

$$\tilde{\phi}_{k}(x) :=\frac{1}{\sqrt{-b_{k}(4 k+6)}} \phi_{k}(x).$$

数值例子

对于前面提到的两点边值问题，取 ，如果取精确解为 ，则右端项为 , .

首先选取 作为基函数.

% LegenSM.m
% Legendre spectral-Galerkin method for the model equation
% -u_xx+u=f in (-1,1) with boundary condition: u(-1)=u(1)=0;
% exact solution: u=sin(kw*pi*x);
% RHS: f=kw*kw*pi^2*sin(kw*pi*x)+sin(kw*pi*x);
% Rmk: Use routines lepoly(); legs(); lepolym();
clear all; close all;
kw=1;
Nvec=[4:2:28];
L2_Err=[]; Max_Err=[];  % initialization error
% Loop for various modes N to calculate numerical errors
for N=Nvec
    [xv,wv]=legs(N+1);      % Legendre-Gauss points and weights
    Lm=lepolym(N,xv);       % matrix of Legendre polynomals
    u=sin(kw*pi*xv);        % exact function
    f=kw*kw*pi^2*sin(kw*pi*xv)+sin(kw*pi*xv);  % Right-hand-side(RHS)
    % Calculting coefficient matrix
    S=diag(4*(0:N-2)+6);       % stiffness matrix
    M=diag(2./(2*(0:N-2)+1)+2./(2*(0:N-2)+5))...
        -diag(2./(2*(0:N-4)+5),2)...
        -diag(2./(2*(2:N-2)+1),-2);  % mass matrix
    A=S+M;
    % Solving the linear system
    Pm=Lm(1:end-2,:)-Lm(3:end,:);  % matrix of Phi(x)
    b=Pm*diag(wv)*f;               % solving RHS
    uh=A\b;                        % expansion coefficients of u_N(x)
    un=Pm'*uh;                     % compositing the numerical solution
    
    L2_err=sqrt(((un-u).^2)'*wv);  % L^2 error
    Max_err=max(abs(un-u));        % maximum pointwise error
    L2_Err=[L2_Err;L2_err];
    Max_Err=[Max_Err;Max_err];
end
% Plot the L^2 and maximum pointwise error
plot(Nvec,log10(L2_Err),'bo-','MarkerFaceColor','w','LineWidth',1)
hold on
plot(Nvec,log10(Max_Err),'rd-','MarkerFaceColor','w','LineWidth',1)
grid on,
legend('L^2 error','L^{\infty} error')
% title('Convergence of Legendre-Galerkin method','fontsize',12)
set(gca,'fontsize',12)
xlabel('N','fontsize', 14), ylabel('log_{10}Error','fontsize',14)

% print -dpng -r600  LegenSM.png
% print -depsc2 LegenSM.eps

进一步，取 ，如果取精确解为 ，则右端项为 , .

选取 作为基函数.

% LegenSM1.m
% Legendre spectral-Galerkin method for the model equation
% -u_xx+u=f in (-1,1) with boundary condition: u(-1)=u(1)=0;
% exact solution: u=sin(kw*pi*x); 
% RHS: f=kw*kw*pi^2*sin(kw*pi*x)+sin(kw*pi*x); 
% Rmk: Use routines lepoly(); legs(); lepolym();
clear all; close all;
kw=10;
Nvec=[32:2:76];
% Initialization for error
L2_Err=[]; Max_Err=[];    
% Loop for various modes N to calculate numerical errors
for N=Nvec
    [xv,wv]=legs(N+1);        % Legendre-Gauss points and weights
    Lm=lepolym(N,xv);         % matrix of Legendre polynomals
    u=sin(kw*pi*xv);          % exact solution
    f=kw*kw*pi^2*sin(kw*pi*xv)+sin(kw*pi*xv);  % RHS
    % Calculting coefficient matrix
    S=eye(N-1);               % stiff matrix
    M=diag(1./(4*(0:N-2)+6))*diag(2./(2*(0:N-2)+1)+2./(2*(0:N-2)+5))...
        -diag(2./(sqrt(4*(0:N-4)+6).*sqrt(4*(0:N-4)+14).*(2*(0:N-4)+5)),2)...
        -diag(2./(sqrt(4*(2:N-2)-2).*sqrt(4*(2:N-2)+6).*(2*(2:N-2)+1)),-2);  % mass matrix
    A=S+M;
    % Solving the linear system
    Pm=diag(1./sqrt(4*(0:N-2)+6))*(Lm(1:end-2,:)-Lm(3:end,:));  % matrix of Phi(x)
    b=Pm*diag(wv)*f;          % solving RHS
    uh=A\b;                   % expansion coefficients of u_N(x)
    un=Pm'*uh;                % compositing the numerical solution
    
    L2_error=sqrt(((un-u).^2)'*wv);  % L^2 error
    Max_error=norm(abs(un-u),inf);   % maximum pointwise error 
    L2_Err=[L2_Err;L2_error];
    Max_Err=[Max_Err;Max_error];
end
% Plot L^2 and maximum pointwise error
plot(Nvec,log10(L2_Err),'bo-','MarkerFaceColor','w','LineWidth',1)
hold on
plot(Nvec,log10(Max_Err),'rd-','MarkerFaceColor','w','LineWidth',1)
grid on
legend('L^2 error','L^{\infty} error','location','NorthEast')
set(gca,'fontsize',12)
xlabel('N','fontsize', 14), ylabel('log_{10}Error','fontsize',14)

% sets axis tick and axis limits
xticks(30:10:80)
yticks(-15:3:0)
xlim([30 80])
ylim([-15 0])

% print -dpng -r600 LegenSM1.png
% print -depsc2 LegenSM1.eps

参考书籍

• Spectral Methods: Algorithms, Analysis and Applications.