计算科学 - 使用 Chebyshev 多项式通过谱法求解 ODE - 吾爱随笔录

我想按照线性弹性的基本方程求解（为了一维的简单）

\frac{d}{d x} (E \frac{d u}{d x}) = 0 with u (1) = 0, u (- 1) = b

$\frac{d}{dx} \left( E \frac{du}{dx} \right) = 0 \quad \textrm{with} \quad u(1)=0, \; u(-1)=b$

按照 Nick Trefethen 的书“Matlab 中的光谱方法”（https://people.maths.ox.ac.uk/trefethen/spectral.html）中第 138 页（程序 33）中的示例，不应该太复杂以适应这与上述指定的弹性方程有关。根据边界条件，我去掉了微分矩阵的第一列第一行。但是，实现出现了问题，因为我的结果振荡得很重（期望线性行为）。如有任何建议，我将不胜感激！

clear all; clc;
N = 16; % Number of collocation points
[D,x] = cheb(N); % Set up differentiation matrix

% Insert boundary condition: u(-1) = 0, u(1) = 2
D2 = D(2:N+1,2:N+1);
D2(N,:) = zeros(1,N);
D2(N,N)=ones(1,1);

% right hand side
f = zeros(N,1);
f(N) = 2.0;
u = (D2.*(D2))\f;
u = [0;u];

% Plot solution
clf
%subplot ('position', [.1 .4 .8 .5])
plot(x,u,'.', 'markersize',16)
xx = -1:0.01:1;
uu = polyval(polyfit(x,u,N+1),xx);
line(xx,uu), grid on

和：

% CHEB compute D = differentitation matrix, x = Chebyshev grid

function [D, x] = cheb(N)
  if N == 0, D = 0; x = 1; return, end
  x = cos(pi*(0:N)/N)';
  c = [2; ones(N-1,1); 2].*(-1).^(0:N)';
  X = repmat(x,1,N+1);
  dX = X-X';
  D = (c*(1./c)')./(dX+(eye(N+1)));
  D = D - diag(sum(D'));