matlab函数主要是计算与坐标差相关的矩阵 `xx` 和 `yy` 及其衍生矩阵

function [C,W,x,y]=GDQ(N,M)
% function [C,W,y]=GDQ(N,n1,m1,n2,m2)

Deta=1e-4; % 0; % 
Deta2=0; % 1e-4; % 
Deta3=0; % 1e-4; % 
% pi=3.1415926;
% if (n1==3 && m1==3)||(m1==3 && n2==3)||(n2==3 && m2==3)||(m2==3 && n1==3)
%     Alfa=0.5;     % FF--含自由角点
% else
%     Alfa=1;       % 不含自由角点
% end
%% ---------------------------------------------------------------
% x=zeros(N);
% y=zeros(N);

for i=1:N
    x(i)=( 1-cos(pi*(i-1)/(N-1)) )/2;
end

% for i=2: N-1
% %     xx=( 1-cos(pi*(i-1)/(N-1)) )/2;               % 思考此处 i、xx后面+Deta与最后+Deta区别
% %     x(i)=(1-Alfa)*( 3*xx^2 - 2*xx^3 ) + Alfa*xx;  % 参考 C Shu书籍：Differential Quadrature and Its Application in Engineering
% %     x(i)=(1-Alfa)*( (1-cos(pi*xx))/2 ) + Alfa*xx; % 自己构造
%     x(i)=x(i) +Deta;
% end

%% ---------------------------------------------------------------
y=x;
for i=1:N
    for j=1:N
        xx(i,j)=[x(i)-x(j)]; 
        yy(i,j)=[y(i)-y(j)];
    end
end

Xx=xx';
Yy=yy';      
Xx(Xx==0)=[];
Yy(Yy==0)=[];
for i=1:N
    for j=1:N-1
        xxx(i,j)=Xx((N-1)*(i-1)+j);
        yyy(i,j)=Yy((N-1)*(i-1)+j);
    end
end

for i=1:N
    Mx(i,:)=cumprod(xxx(i,:));   %cumprod 累积
    My(i,:)=cumprod(yyy(i,:));
end
Mx=Mx(:,N-1);
My=My(:,N-1);
C(:,:,1)=zeros(N,N);
W(:,:,1)=zeros(N,N);
for i=1:N
    for j=1:N
        if j~=i
            C(i,j,1)=Mx(i)/xx(i,j)/Mx(j);        % 一阶导数权重系数，见公式（6）
            W(i,j,1)=My(i)/yy(i,j)/My(j);
        end
    end
    C(i,i,1)=-sum(C(i,:,1)) +Deta2; %
    W(i,i,1)=-sum(W(i,:,1)) +Deta2; %
end
% C
% W
for m=2:N-1
    C(:,:,m)=zeros(N,N);
    W(:,:,m)=zeros(N,N);
    for i=1:N
        for j=1:N
            if j~=i
                C(i,j,m)=m*(C(i,j,1)*C(i,i,m-1)-C(i,j,m-1)/xx(i,j));   %二阶及其以上的导数权重系数，见公式（7）
                W(i,j,m)=m*(W(i,j,1)*W(i,i,m-1)-W(i,j,m-1)/yy(i,j));
            end
        end
        C(i,i,m)=-sum(C(i,:,m)) +Deta3; %   %所有导数，见公式（8）
        W(i,i,m)=-sum(W(i,:,m)) +Deta3; %
    end
end
% C
% W

以下是对该函数 GDQ 的功能解释：

一、函数输入参数：

• N：表示要处理的数据点的数量。
• M：具体功能在函数中未体现，可能是后续功能扩展的预留参数，目前代码中未使用。

二、函数内部主要操作：

1. 变量初始化：

   • 首先，将 Deta、Deta2 和 Deta3 初始化为一些小的数值（目前部分为 0，部分为 $10^{-4}$）。
   • 初始化 x 和 y 为长度为 N 的零向量。

2. x 坐标的生成：

   • 通过循环从 i = 1 到 N，使用公式 x(i)=( 1-cos(pi*(i-1)/(N-1)) )/2 计算 x 坐标的值。
   • 代码中存在一些注释掉的部分，这些部分可能是不同的 x 坐标生成方法的尝试，如使用多项式或其他函数的组合，但最终未被使用。

3. y 坐标的生成：

   • 简单地将 y 赋值为 x，使得 y 坐标与 x 坐标相同。

4. 差值矩阵的计算：

   • 计算 xx 和 yy 矩阵，它们存储了 x 和 y 坐标之间的差值，即 xx(i,j)=[x(i)-x(j)] 和 yy(i,j)=[y(i)-y(j)]。
   • 对 xx 和 yy 进行转置得到 Xx 和 Yy，并将其中等于 0 的元素移除。
   • 将 Xx 和 Yy 重新排列存储到 xxx 和 yyy 矩阵中，其中 xxx(i,j)=Xx((N-1)*(i-1)+j) 和 yyy(i,j)=Yy((N-1)*(i-1)+j)。

5. 累积积计算：

   • 对 xxx 和 yyy 矩阵按行进行累积积计算得到 Mx 和 My 矩阵，Mx(i,:)=cumprod(xxx(i,:)) 和 My(i,:)=cumprod(yyy(i,:))，最终只保留每行的最后一个累积积结果，即 Mx=Mx(:,N-1) 和 My=My(:,N-1)。

6. 权重系数计算：

   • 对于一阶导数的权重系数：

     • 初始化 C 和 W 矩阵为 N x N x 1 的零矩阵。
     • 对于 i 和 j 不相等的情况，根据公式 C(i,j,1)=Mx(i)/xx(i,j)/Mx(j) 和 W(i,j,1)=My(i)/yy(i,j)/My(j) 计算 C 和 W 的元素。
     • 对于 i = j 的元素，使用公式 C(i,i,1)=-sum(C(i,:,1)) +Deta2 和 W(i,i,1)=-sum(W(i,:,1)) +Deta2 计算，其中 sum(C(i,:,1)) 是对第 i 行的元素求和，并添加 Deta2 作为调整项。

   • 对于二阶及以上导数的权重系数（m 从 2 到 N-1）：

     • 扩展 C 和 W 矩阵为 N x N x m 的零矩阵。
     • 对于 i 和 j 不相等的情况，根据公式 C(i,j,m)=m*(C(i,j,1)*C(i,i,m-1)-C(i,j,m-1)/xx(i,j)) 和 W(i,j,m)=m*(W(i,j,1)*W(i,i,m-1)-W(i,j,m-1)/yy(i,j)) 计算元素。
     • 对于 i = j 的元素，使用公式 C(i,i,m)=-sum(C(i,:,m)) +Deta3 和 W(i,i,m)=-sum(W(i,:,m)) +Deta3 计算，其中 sum(C(i,:,m)) 是对第 i 行的元素求和，并添加 Deta3 作为调整项。

三、函数输出：

• C：存储不同阶导数的 x 方向的权重系数矩阵，是一个 N x N x M 的三维矩阵。
• W：存储不同阶导数的 y 方向的权重系数矩阵，是一个 N x N x M 的三维矩阵。
• x：存储生成的 x 坐标的向量，长度为 N。
• y：存储生成的 y 坐标的向量，长度为 N。

四、总结：
该函数主要是计算与坐标差相关的矩阵 xx 和 yy 及其衍生矩阵，然后通过累积积和一系列公式计算不同阶导数的权重系数矩阵 C 和 W，可能用于求解某些微分方程的数值近似，例如使用微分求积法（Differential Quadrature Method）对偏微分方程进行离散化求解。在计算过程中，引入了一些小的调整项 Deta、Deta2 和 Deta3，可能是为了避免除数为零或改善计算的稳定性和精度等。对于 y 坐标的处理目前比较简单，只是复制了 x 坐标，后续可能可以扩展为不同的计算方式。