Spectral Method 2011版学习日记 1.Introduction

1.Introduction
1.5.1 Finite-Difference Versus Spectral-Collocation P16

初学小白，书16页数值实现，有限差分法和谱配点法的对比。
生疏了，很不精简，有什么奇怪的地方，就忽略它吧！

function p16
clear;clc;
u0=inline('log(2+sin(x))','x');
du0=inline('(cos(x))./(2+sin(x))','x');%点除点乘啊啊啊，不然代入就是输出矩阵啊啊
% algebraic convergence for FEM and FDM
figure('Position',[0 0 1000 500] ), subplot(1,2,1); 
for Nvec =10:4:100;
errfd =fdm(u0,du0,Nvec);%有限差分
errsp= spm(u0,du0,Nvec);
loglog(Nvec, errfd, 'r+', 'markersize', 8), hold on%双对数坐标变换
loglog(Nvec, errsp, 'bo', 'markersize', 8), hold on
end
title('Convergence of 4th-order FD &Fourier-collection methed');
grid on, xlabel N, ylabel error
for Nvec =10:0.1:100;
    loglog(Nvec, (2.*pi./Nvec).^4, 'b-','markersize', 14),hold on%对数坐标系中h^{4}
end
text(100, 1e-5, 'h^{4}', 'fontsize', 18);

% spectral convergence for spectral method
subplot(1,2,2);%1行2列第二图
for Nvec =2:4:50;%图二横坐标 
errsp= spm(u0,du0,Nvec);%Fourier 谱方法
[Nvec errsp];
semilogy(Nvec, errsp, 'r+', 'markersize', 8)
hold on;%y坐标系为对数坐标系
end
grid on, xlabel N, ylabel error
title('Convergence of Fourier-collection method');


end


function [u,du,x]= EvalSolEx(uinit,duin,N)
    h=2*pi./N;  x=(0:h:2*pi-h)';%xi
    u=feval(uinit, x);%已知数据代入函数句柄，uinit即u0
    du=feval(duin,x);
end

function  errsp = spm(uinit,du0,N)%Fourier spectral
%uinit原函数,Nvec取点数
      [u,du,x] = EvalSolEx(uinit,du0,N);
      h=2*pi./N;
      H=zeros(N,N);
      for k=0:1:N-1
          for j=0:1:N-1
              if k~=j
                  H(j+1,k+1)=(1/2).*((-1).^(k+j))*(cot((j-k).*pi./N));
              else
                  H(j+1,k+1)=0;
              end
              
          end
      end
      w=H*u;
      errsp= norm(du-w, 'inf');%返回矩阵的无穷范数，即最大一行的和 
end

% algebraic convergence for f1
function err= fdm(uinit,duin,N)%差分
       errfd =[];
       A=diag(repmat([0], 1, N))+diag(repmat([8],1, N-1), 1)+diag(repmat([-1], 1, N-2), 2)+diag(repmat([-8], 1, N-1), -1)+diag(repmat([1], 1, N-2), -2);
       A(1,N-1)=1;A(2,N)=1;
       A(1,N)=-8;
       A(N-1,1)=-1;A(N,2)=-1;
       A(N,1)=8;
       [u,du,x] = EvalSolEx(uinit,duin, N);
       h=2.*pi./N;
       w=(1./(12.*h))*A*u;
       err = norm(du-w, 'inf');
       errfd = [errfd err];
end