Legendre 配置法及其算例

谱配置方法也称配点法，选取节点为 Legendre 高斯节点称为 Legendre 配置方法。这里介绍两点边值问题的 Legendre 配置方法，然后 MATLAB 编程算出数值解以及误差。

两点边值问题

考虑两点边值问题：

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

Legendre Collocation 方法

在区间 上选取配置点 为 Legendre-Gauss-Lobatto 节点，其中 ，. 以上方程的配置法 (\label{equcm})

$$\left\{\begin{aligned} &\text{Find } u \in P_N \text{ such that} \\ &-\mu u_{N}^{\prime \prime}(x_{i})+\nu u_{N}^{\prime}(x_{i})+\rho u_{N}(x_{i})=f(x_{i}), \quad 1 \leq i \leq N-1, \\ &u_N(-1)=0, \quad u_N(1)=0. \end{aligned} \right.$$

逼近解展开

$$u_{N}(x)=\sum_{j=0}^{N} u_{N}\left(x_{j}\right) h_{j}(x).$$

其中 是 Lagrange 基多项式（也称为节点基函数），即 ，，令 . 记 , 将 展开为 ， 可知

$$u_{N}\left(x_{k}\right)=\sum_{j=0}^{N} w_{j} h_{j}\left(x_{k}\right)=w_{k},$$ $$\begin{aligned} u_{N}^{\prime}(x_{k})&=\sum_{j=0}^{N} w_{j} h_{j}^{\prime}(x_{k})=\sum_{j=0}^{N} d_{kj} w_{j} \\&=\sum_{j=1}^{N-1} d_{kj} w_{j}+d_{k0} w_{0}+d_{kN} w_{N}, \end{aligned}$$ $$\begin{aligned} u_{N}^{\prime \prime}(x_{k})&=\sum_{j=0}^{N} w_{j} h_{j}^{\prime \prime}(x_{k})=\sum_{j=0}^{N}(D^{2})_{k j} w_{j} \\ &=\sum_{j=1}^{N-1}\left(D^{2}\right)_{k j} w_{j}+\left(D^{2}\right)_{k 0} w_{0}+\left(D^{2}\right)_{k N} w_{N}. \end{aligned}$$

将以上公式代入方程 (equcm) 可得

$$\left\{\begin{aligned} &\sum_{j=0}^{N}\left[-\mu \left(D^{2}\right)_{ij}+\nu d_{i j}+\rho \delta_{i j}\right] w_{j}=f(x_{i}), \quad i=1,2, \cdots, N-1, \\ &w_{0}=0, \quad w_{N}=0. \end{aligned} \right.$$

记

$$\begin{aligned} & a_{i j}=-\mu (D^{2})_{ij}+\nu d_{ij}+\rho \delta_{ij}, \quad 1 \leq i \leq N-1,\, 0 \leq j \leq N, \\ &a_{0 j}=\delta_{0 j},\quad a_{N j}=\delta_{N j}, \quad 0 \leq j \leq N, \\ &\mathbf{b}=(0, f(x_{1}), f(x_{2}), \cdots, f(x_{N-1}), 0)^{T}, \\ &\mathbf{w}=(w_{0}, w_{1}, \ldots, w_{N})^{T}, \quad A=\left(a_{i j}\right)_{0 \leq i, j \leq N}. \end{aligned}$$

线性系统可以简化为

$$A\mathbf{w}=\mathbf{b}.$$

数值例子

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

% LegenCM2.m
% Legendre-collocation method for the model equation: 
% -mu*u''(x)+nu*u'(x)+rho*u(x)=f(x),  x in (-1,1),
% boundary condition: u(-1)=u(1)=0;
% exact solution: u=sin(kw*pi*x); 
% RHS: f(x)=mu*kw^2*pi^2*sin(kw*pi*x)+nu*kw*pi*cos(kw*pi*x)+rho*sin(kw*pi*x)
% Rmk: Use routines lepoly(); legslb(); legslbdm(); 
clear all; close all;
kw=10;
nu=1;
mu=1;
rho=1;
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]=legslb(N);    % Legendre-Gauss-Lobatto points and weights
    u=sin(kw*pi*xv);      % exact solution
    f=mu*kw*kw*pi^2*sin(kw*pi*xv)+nu*kw*pi*cos(kw*pi*xv)+rho*sin(kw*pi*xv);  % RHS
    % Solve the collocation system
    D1=legslbdiff(N,xv);  % 1st order differentiation matrix
    %D1=legslbdm(N);      % 1st order differentiation matrix
    D2=D1*D1;             % 2nd order differentiation matrix
    % Compositing the coefficient matrix
    D=-mu*D2(2:N-1,2:N-1)+nu*D1(2:N-1,2:N-1)+rho*eye(N-2);
    b=f(2:N-1);           % RHS
    un=D\b;
    un=[0;un;0];
    
    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)
% title('Convergence of Legendre-collocation method','fontsize',12)

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

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

参考书籍

• Spectral Methods: Algorithms, Analysis and Applications