关于使用dct求解零诺依曼边界条件PDE的一点说明

DCT的第二型定义

如我的这篇文章中所述，对于纽曼边界条件的偏微分方程，谱方法中，我们一般使用离散余弦变换，即dct来完成。对于fft和dst而言，似乎定义比较唯一。但是对于，对于dct，它的定义形式就有很多种，维基百科中就列出了它的八种形式。

一般常用的形式有两种，即DCT-I：

二维形式形如下面这样：

求和符号的两个撇撇表示求和符号中当m取到0或者J时，前面乘个系数0.5，对于n也是这样。

和DCT-II：

在Matlab中使用的便是DCT-II的定义，截图如下：

一般来说，第二种形式在图像处理等方面更为常见一些，然而，在用快速变换求解PDE时，有些情况下，可能用此定义并不满足你的要求。

第二定义是原数组关于N-1/2点展开，做离散傅里叶变换得到的，第一定义是原数组关于N点展开得到的。那么，在使用谱方法时他们所要求解的坐标节点是不一样的。什么意思呢？比如对于纽曼边界条件问题，定义域等分成N分，节点标号从0到N，共N+1个节点，包括边界。那么，第一定义求解的便是u0,u1,...,uNu_0,u_1,...,u_Nu0,u1,...,uN共N+1个值，而第二型定义求解的是网格中点处的值，即u1/2,u3/2,...,uN−1/2u_{1/2},u_{3/2},...,u_{N-1/2}u1/2,u3/2,...,uN−1/2共N个值。

很多人不明白这一点，瞎用matlab内置的dct来求解网格节点上的值，这是不对的。当你需要第二型定义来求解PDE在剖分节点上的值（包括边界）时，就需要自己写一个dct来求解。

简单用求和的定义写出的dct，idct，dct2，idct2如下，写得很挫，因为不用这个，所以随手写的：

function [ F ] = mydct( f )
%一般的列向量，我们喜欢用列向量
flag = 0;
if size(f,1)==1f = f';flag = 1;
end
N = length(f)-1;
js = 0:N;
ks = 0:N;
[js_mg,ks_mg] = meshgrid(js,ks);%前后保持一致
cos_mat = cos(pi.*js_mg.*ks_mg./N);
f(1) = f(1)*0.5;
f(end) = f(end)*0.5;F = cos_mat*f;if flag==1F = F';
endend

function [ F ] = myidct( f )
%一般的列向量，我们喜欢用列向量
flag = 0;
if size(f,1)==1f = f';flag = 1;
end
N = length(f)-1;
js = 0:N;
ks = 0:N;
[js_mg,ks_mg] = meshgrid(js,ks);%前后保持一致
cos_mat = cos(pi.*js_mg.*ks_mg./N);
f(1) = f(1)*0.5;
f(end) = f(end)*0.5;
F = cos_mat*f;F = F*2/N;if flag==1F = F';
endend

function B = mydct2( A )
[M,N] = size(A);
ms = 1:M;
ns = 1:N;
js = 1:M;
ls = 1:N;
[ms_grid,ns_grid] = meshgrid(ns,ms);
%[ms_grid,ns_grid] = meshgrid(ms,ns);
A(1,:) = A(1,:)*0.5;
A(end,:) = A(end,:)*0.5;
A(:,1) = A(:,1)*0.5;
A(:,end) = A(:,end)*0.5;for j = jsfor l = lsI = A.*cos(pi*(j-1).*(ms_grid-1)/(M-1)).*cos(pi*(l-1).*(ns_grid-1)/(N-1));B(l,j) = sum(sum(I));end
end
end

function B = myidct2( A )
[M,N] = size(A);
ms = 1:M;
ns = 1:N;
js = 1:M;
ls = 1:N;
[ms_grid,ns_grid] = meshgrid(ns,ms);
A(1,:) = A(1,:)*0.5;
A(end,:) = A(end,:)*0.5;
A(:,1) = A(:,1)*0.5;
A(:,end) = A(:,end)*0.5;for j = jsfor l = lsI = A.*cos(pi*(j-1).*(ms_grid-1)/(M-1)).*cos(pi*(l-1).*(ns_grid-1)/(N-1));B(l,j) = sum(sum(I))*4/((M-1)*(N-1));end
end
end

这么写是跑起来很慢的，我们在用matlab的时候，一个基本的原则是，有内置的函数，尽量用内置的。所以，我们根据dct的定义的来由，根据fft来计算dct，如下：

function y=dctbyfft(x)
% Discrete Cosine Transform DCT-I
[M,N]=size(x);
y=[x;flipud(x(2:M-1,:))];
yy=fft(y);
y=real(yy(1:M,:)/2);

function y=idctbyfft(x)
% Inverse Discrete Cosine Transform IDCT-I
M=size(x,1);
y=2/(M-1)*dctbyfft(x);

function [ Y ] = dct2byfft( X )
Y = dctbyfft(dctbyfft(X)')';
end

function [ Y ] = idct2byfft( X )
Y = idctbyfft(idctbyfft(X)')';
end

我们可以随意地去测试它们，保证两者的结果是一致，说明没有错误：

% a = ([1:20000]');
% myidct(a);
% idctbkl(a);
% mydct(a);
% dctbkl(a);
%A = [4 2;2 4];
A = magic(6);
B1 = myidct2(A)
B2 = idct2byfft(A)

泊松问题下的fft方法分母为零问题

在差分-傅里叶谱方法实际应用当中，如果分母出现了0，怎么办？我看到过别人直接找一个ϵ\epsilonϵ去替代0的，也有直接令那个无穷大的数为0的，我并不知道这是什么机理，不过，既然人家这么搞了，没问题，那我们不妨暂时先接受它。

差分-傅里叶变换法求解泊松方程

举个泊松方程的例子，从代码中，你们可以读出这个方法使用的一些细节，比如说下标从哪到哪等等。

function u=laplacefft(u0,bndcond)
%LAPLACEFFT Solve Laplace's Equation using Fourier's Method.
%   U = LAPLACEFFT(U0) solves Laplace's equation on a rectangle with Dirichlet
%   boundary conditions on all four boundaries ("given values"). The size of
%   the domain and the boundary conditions are given by U0. Note: Only the
%   boundary values of U0 are actually used.
%
%   U = LAPLACEFFT(U0,'DIRICHLET') is equivalent to the above.
%
%   U = LAPLACEFFT(U0,'NEUMANN') uses Neumann conditions ("insulation") on the
%   top and bottom boundaries. Note: Only U0(:,1) and U0(:,end) are used in
%   this case.
%
%   The solver might run faster if the dimensions of U0 are 2^p+1 for
%   integer p (as in the example below).
%   References:
%      G. Strang, "Introduction to Applied Mathematics", Wellesley-
%      Cambridge Press, 1986. (Section 5.5)
%
%      W. Press, S. Teukolsky, W. Vetterling, B. Flannery, "Numerical
%      Recipes in C", 2nd Edition, Cambridge University Press, 1992.
%      (Section 19.4).dirichlet=1;
if nargin>=2switch upper(bndcond)case 'DIRICHLET'case 'NEUMANN'dirichlet=0;otherwiseerror('The second input argument must be ''DIRICHLET or ''NEUMANN''.');end
end[M,N]=size(u0);
if dirichlet % Dirichlet on all four boundariesf=zeros(M-2,N-2);%先求非边界的f，非零边界条件的处理，就是f减去边界条件，就成了边界条件的处理f(:,1)=f(:,1)-u0(2:M-1,1); %f的在内部的第一列减去内部的第一列f(:,N-2)=f(:,N-2)-u0(2:M-1,N); %f在内部的倒数第一列等于自身减去u0的右边界f(1,:)=f(1,:)-u0(1,2:N-1); %f的第一行等于自身减去u0的第一行f(M-2,:)=f(M-2,:)-u0(M,2:N-1);%f的倒数第一行等于自身减去下边界K=2*(cos(pi*(1:M-2)'*ones(1,N-2)/(M-1))+cos(pi*ones(M-2,1)*(1:N-2)/(N-1))-2);fhat=dst(dst(f)')';uhat=fhat./K;u1=idst(idst(uhat')');u=u0;u(2:M-1,2:N-1)=u1;
else % Dirichlet on left and right, Neumann on top and bottomf=zeros(M,N-2);f(:,1)=f(:,1)-u0(:,1);f(:,N-2)=f(:,N-2)-u0(:,N);K=2*(cos(pi*(0:M-1)'*ones(1,N-2)/(M-1))+cos(pi*ones(M,1)*(1:N-2)/(N-1))-2);fhat=dst(dct(f)')';uhat=fhat./K;u1=idst(idct(uhat)')';u=u0;u(:,2:N-1)=u1;
end

% clc
% clear
% close all;
h=1/8;%定义步长为（0.5^5）
[xx,yy]=meshgrid(0:h:2,0:h:1);%横纵坐标的一个meshgrid
U0=xx.^2-yy.^2;     %精确解及边值
U=laplacefft(U0);  %用傅里叶方法求解
contourf(xx,yy,U),view(2),axis equal,colorbar %画图，画等高线图
error=max(abs(U(:)-U0(:))) %求无穷误
[M,N] = size(U);
error=sum(sum(abs(U(:)-U0(:))))/(M*N)

还有一些别的例子，我也贴上来：

function u = poisson1Dneumann(F,x0,xEnd)
% POISSON1DNUEMANN solves the 1D poisson equation d2U/dX2 = F, with neumann
% boundary conditions dUdX = 0 at X = 0 and X = L.
% u = poisson1Dneumann(F,x0,xEnd)
%
% u: Vector of solution
% F: Vector of right-hand-side
% x0: Begining coordinate of domain.
% xEnd: End coordinate of domain.
%% Check for compatibility
xInt = linspace(x0,xEnd,length(F));
fInt = trapz(xInt,F);
if (fInt > 0.0001) || (fInt < -0.0001)disp('Does not satisfy compatibility conditions');
end% Solution
N = length(F);
dx = (xEnd - x0) / (N - 1);b = dct(F);
m = (0:length(b)-1)';
a = dx^2 * b ./ (2 * (cos(m * pi / (N - 1)) - 1))';
a(1) = 0;
u = idct(a);

function U = poisson2Dneumann(F,L)
% POISSON2DNEUMANN solves the the 2D poisson equation d2UdX2 + d2UdY2 = F,
% with the zero neumann boundary condition on all the side walls.
% U = POISSON2DNEUMANN(F,L)
% U: Solution of d2UdX2 + d2UdY2 = F
% F: Right hand side matrix of size N*N
% L: Length of the domain in either X or Y direction. It should be noted
% that Lx = Ly = LN = length(F);
dx = L / (N - 1);
p = 0:N-1;
q = 0:N-1;
[p,q] = meshgrid(p,q);
B = dct2(F);
A = dx^2 * B ./ (2 * cos(pi * p / N) + 2 * cos(pi * q / N) - 4); A(1,1) = 0;
U = idct2(A);

更多的细节，可以参考：

Press W H, Teukolsky S A, Vetterling W T, et al. Numerical recipes 3rd edition: The art of scientific computing[M]. Cambridge university press, 2007.