MATLAB求解1D Poisson方程数值解与解析解不符的问题排查

核心错误原因 #

1. 绘图x轴未对齐
   原代码中，plot(u(X))使用[0,1]区间的均匀内部网格点X作为x轴，但plot(linsolve(A, b))未指定x坐标，MATLAB默认用1到m的索引作为x轴，导致两条曲线的x轴范围完全错位，视觉上呈现巨大差异——这并非数值解本身错误，而是绘图坐标系不匹配。

2. tridiag函数正确性存疑
   MATLAB无内置tridiag函数，若自定义函数的参数顺序错误（比如主对角线、次对角线顺序颠倒），会生成错误的Poisson矩阵，直接导致数值解偏离解析解。

修正后的代码 #

m = 32;
h = 1/(m+1);
X = h:h:(1-h);

% 定义右端项f(x)
f = @(X) 25*pi^2*(sin(5*pi*X) + 9*sin(15*pi*X));
b = h^2 * f(X)';  % 生成线性系统右端向量

% 解析解
u = @(X) sin(5*pi*X) + sin(15*pi*X);

% 直接生成正确的三对角Poisson矩阵（无需依赖自定义tridiag）
A = diag(2*ones(m,1)) + diag(-1*ones(m-1,1),1) + diag(-1*ones(m-1,1),-1);
% 大m场景建议用稀疏矩阵提升效率：A = sparse(diag(2*ones(m,1)) + diag(-1*ones(m-1,1),1) + diag(-1*ones(m-1,1),-1));

hold on
plot(X, u(X), 'b-', LineWidth=1.5)  % 指定x轴为网格点X
plot(X, linsolve(A, b), 'r--', LineWidth=1.5)  % 对齐x轴绘制数值解
legend('解析解', '数值解')
xlabel('x')
ylabel('u(x)')
grid on
hold off

验证说明 #

修正x轴对齐后，数值解与解析解会几乎完全重合（误差在中心差分法的截断误差范围内）。若仍有偏差，需检查自定义tridiag函数的实现逻辑，确保生成的矩阵主对角线为2，上下次对角线为-1，对应离散化后的方程：
$$-u_{i-1} + 2u_i - u_{i+1} = h^2 f(x_i)$$

内容的提问来源于stack exchange，提问作者Cymek3