计算科学 - 如何在 MATLAB 中形成切比雪夫微分矩阵？ - 吾爱随笔录

如何在 MATLAB 中形成切比雪夫微分矩阵？

计算科学 matlab 数值分析谱法

2021-12-18 05:51:44

我有一些代码可以做到这一点，但我不喜欢使用我不理解的东西。这是代码

N1=N+1;
cheb=cos(pi*(0:N)/N)';
unif=linspace(-1,1,N1)';
if N<3
    x=cheb;
else
    x=cheb+sin(pi*unif)./(4*N);
end
P=zeros(N1,N1);
%N1xN1 zero matrix
xold=2;
%eps= epsilon!
while max(abs(x-xold))>eps

xold=x;   
P(:,1)=1;
P(:,2)=x;
%set first collumn entries to 1 (P(:,1) = 1);
%set second collumn entries to x (P(:,2) = x);
for k=2:N
    P(:,k+1)=( (2*k-1)*x.*P(:,k)-(k-1)*P(:,k-1) )/k;
    %%Bonnets formula ;)
end
x=xold-( x.*P(:,N1)-P(:,N) )./( N1*P(:,N1));
end

%---chebyshev differentiation matrix----------------------------------------    ---------
x=-x;
X=repmat(x,1,N1);%set every collumn of X to x
Xdiff=X-X'+eye(N1);
L=repmat(P(:,N1),1,N1);
L(1:(N1+1):N1*N1)=1;
D=(L./(Xdiff.*L'));
D(1:(N1+1):N1*N1)=0;
D(1)=-(N1*N)/4;
D(N1*N1)=(N1*N)/4;

从理论上讲，我真的很想了解这个矩阵是如何形成的，然后更重要的是，如何在 MATLAB 中实现。

直到评论

% --- chebyshev differenti....

我明白。

我似乎找不到任何与这个矩阵的形成有关的文件。

注意：我可以接受的语法，它是我遇到问题的矩阵的实际形成。我已经在 MSE 上发布了这个，它被搁置了。我查看了帮助中心，这个网站似乎符合描述。

1个回答

我假设您知道 Chebyshev 搭配方法是如何工作的（但如果不知道，请告诉我，我会再解释一下）；一个很好的介绍是 Nick Trefethen在 Matlab 中的光谱方法以及他的近似理论和近似实践（特别是第 21 章）。（但请注意，此代码使用勒让德搭配，而不是切比雪夫搭配！）

使用插值多项式的拉格朗日形式，您可以推导出微分矩阵的显式形式 $D$ ：

D_{i j} = ℓ_{j}^{'} (x_{i}) = {\begin{cases} \frac{λ_{j}}{λ_{i} (x_{i} - x_{j})} & for i \neq j, \\ \frac{x_{j}}{1 - x_{j}^{2}} & for i = j, \end{cases}

$D_{ij} = \ell'_j(x_i) = \begin{cases} \frac{\lambda_j}{\lambda_i(x_i-x_j)} & \text{for }i\neq j,\\ \frac{x_j}{1-x_j^2} & \text{for }i= j, \end{cases}$ 在哪里

ℓ_{j}

$\ell_j$ 是朗朗日节点多项式（即

ℓ_{j} (x_{k}) = δ_{j k}

$\ell_j(x_k) = \delta_{jk}$ ）和

λ_{j}

$\lambda_j$ 是 Langrage 多项式中的（倒数）分母

ℓ_{j}

$\ell_j$ . 这正是最终代码块计算的内容（以特别紧凑的形式）：

Xdiff=X-X'+eye(N1)计算所有可能的差异 $x_i$ 和 $x_j$ 并为 $i=j$ （即，在对角线上）
L=repmat(P(:,N1),1,N1)计算分母 $\lambda_j$ 从切比雪夫范德蒙德矩阵
L(1:(N1+1):N1*N1)=1用一个替换对角线条目（因为没有 $\lambda_j$ 为了 $i=j$ )
D=(L./(Xdiff.*L'))现在计算所有条目 $D_ij$ 同时（对于 $i\neq j$ , 特例 $i=j$ 由上述修改照顾）
其余行修改D以考虑边界条件（代码仅计算与内部点对应的矩阵）。

（供参考，原始代码来自这里。）

其它你可能感兴趣的问题

上一篇AMG 和 Gauss Seidel 的 Python 环境作为求解器而不是预处理器下一篇3D 多边形处理库