Chebyshev 谱方法及其算例

谱方法的正交多项式为 Chebyshev 多项式，为 Chebyshev 谱方法。这里介绍一维问题的 Chebyshev 谱方法，然后 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.$$

Chebyshev-Galerkin 方法

Chebyshev 多项式的权函数 .

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

$$\phi_{k}(x)=T_{k}(x)-T_{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)_{\omega}+\alpha (u_N,v_N)_{\omega}=(f,v_N)_{\omega},\quad v_N \in X_N. \end{aligned} \right.$$

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

$$\begin{aligned} &f_{k}=\int_{I} f_{N} \phi_{k} \omega 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} \omega d x, \quad m_{k j}=\int_{I} \phi_{j} \phi_{k} \omega 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}.$$

质量矩阵 是一个对称正定的五对角矩阵 [P149-4.29]，它的非零元素为

$$m_{j k}=m_{k j}=\left\{ \begin{aligned} &\frac{\pi}{2}\left(c_{k}+a_{k}^{2}+b_{k}^{2}\right), && j=k, \\ &\frac{\pi}{2}\left(a_{k}+a_{k+1} b_{k}\right), && j=k+1, \\ &\frac{\pi}{2} b_{k}, && j=k+2. \end{aligned}\right.$$

其中 ，当 时，.

刚度矩阵 是一个上三角矩阵 [P149-4.30]，其元素如下：

$$s_{k j}=\left\{ \begin{aligned} &2 \pi(k+1)(k+2), && j=k, \\ &4 \pi(k+1), && j=k+2, k+4, k+6, \cdots, \\ &0, && j<k \text { or } j+k \text { odd}. \end{aligned}\right.$$

数值例子

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

% ChebSM1.m
% Chebyshev spectral-Galerkin method for the model equation
% -u_xx+u=f in (-1,1)
% 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 japoly(); jags(); japolym();
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]=jags(N+1,-1/2,-1/2);   % Chebyshev-Gauss points and weights
    Cm=japolym(N,-1/2,-1/2,xv)./japolym(N,-1/2,-1/2,1);  % matrix of Chebyshev 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=zeros(N-1);    % stiff matrix
    for k=1:N-1
        for j=1:N-1
            if k==j
                S(k,j)=2*pi*k*(k+1);
            elseif (k<j && mod(j-k,2)==0)
                S(k,j)=4*pi*k;
            else
                S(k,j)=0;
            end
        end
    end
    
    M=diag([3/2*pi, pi*ones(1,N-2)])+diag(-pi/2*ones(1,N-3),2)...
        +diag(-pi/2*ones(1,N-3),-2);    % mass matrix
    A=S+M;
    
    % Solving the linear system
    Pm=Cm(1:end-2,:)-Cm(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=norm(abs(un-u),inf);   % maximum pointwise error
    L2_Err=[L2_Err;L2_err];
    Max_Err=[Max_Err;Max_err];
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')
% title('Convergence of Chebyshev-Galerkin method','fontsize',12)
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(-14:2:0)
xlim([30 80])
ylim([-14 0])

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

参考书籍

• Spectral Methods: Algorithms, Analysis and Applications