Andy123t
0%

Legendre 谱方法及其算例

发表于 更新于 分类于

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

两点边值问题

考虑两点边值问题：

Legendre-Galerkin 方法

弱形式：

令 $\phi_{k}(x)=L_{k}(x)+a_{k} L_{k+1}(x)+b_{k} L_{k+2}(x)$ 满足边界条件，[Shen2011, P146] 可得 $a_{k}=0$, $b_{k}=-1$，所以

设 $\{\phi_{k}\}_{k=0}^{N-2}$ 是 $X_{N}$ 的一组基函数，我们将逼近解展开为

谱格式：

将 $u_{N}(x)$ 的展开式代入谱格式方程，取 $v_N=\phi_{k}$。令

可得以下方程组

刚度矩阵 $S=(s_{jk})$ 是一个对角矩阵 [P146-4.22]，它的非零元素为

质量矩阵 $M=(m_{jk})$ 是一个对称的五对角矩阵 [P146-4.23]，它的非零元素为

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

数值例子

对于前面提到的两点边值问题，取 $\alpha=1$，如果取精确解为 $\sin (k \pi x)$，则右端项为 $f=k^2 \pi^2 \sin(k\pi x)+\sin(k\pi x)$, $k=1$.

首先选取 ${\phi}_{k}(x)$ 作为基函数.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
% 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

LegenSM.png

进一步，取 $\alpha=1$，如果取精确解为 $\sin (k \pi x)$，则右端项为 $f=k^2 \pi^2 \sin(k\pi x)+\sin(k\pi x)$, $k=10$.

选取 $\tilde{\phi}_{k}(x)$ 作为基函数.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
% 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

LegenSM1.png

参考书籍

  • Spectral Methods: Algorithms, Analysis and Applications.