基于有限差分法的Klein-Gordon方程Matlab求解错误排查

问题背景 #

尝试用有限差分法求解Klein-Gordon方程，并用解析解验证结果，但当前数值解与解析解图形不一致，需修正代码。

错误分析 #

原代码存在以下核心问题：

• 差分格式推导错误：未正确离散Klein-Gordon方程$W_{tt} = a^2 W_{xx} + b W$，迭代公式混淆了时间与空间的计算顺序，索引逻辑混乱。
• 初始条件设置错误：错误地将解析解的边界空间点所有时间值作为初始条件，而非初始时刻的空间分布及初始时间导数对应的第二个时间层。
• 循环迭代逻辑错误：时间步循环中引用了未计算的未来时间层（j+1），导致计算逻辑完全偏离正确的时间推进方向。

修正后代码 #

close all
clear 
clc

%% 常量参数初始化
A      = 2; 
B      = 3; 
lambda = 2; 
mu     = 3; 
a      = 4; 
b      = mu^2 - (lambda^2)/a^2;  % 修正表达式顺序，与方程定义一致

%% 时空坐标系构建
number_of_discrete_time_steps = 300; 
t = linspace(0, 2, number_of_discrete_time_steps); 
dt = t(2) - t(1); 

number_of_discrete_space_steps = 100; 
x = transpose(linspace(0, 1, number_of_discrete_space_steps)); 
dx = x(2) - x(1); 

%% 解析解计算与绘制
Wa = cos(lambda * x) * (A * cos(mu * t) + B * sin(mu * t));

figure('Name', '解析解'); 
surface(t, x, Wa, 'edgecolor', 'none'); 
colormap(jet(256)); 
colorbar; 
xlabel('t'); 
ylabel('x'); 
title('Wa(x, t) - 解析解');

%% 数值解计算
Wn = zeros(number_of_discrete_space_steps, number_of_discrete_time_steps);

% 设置初始条件：t=0时的空间分布
Wn(:, 1) = cos(lambda * x) * A;  
% 设置t=dt时的空间分布（泰勒展开，利用初始时间导数W_t(x,0)=cos(lambda*x)*B*mu）
Wn(:, 2) = Wn(:, 1) + dt * cos(lambda * x) * B * mu + 0.5 * dt^2 * (a^2 * (-lambda^2 * cos(lambda * x)*A) + b * Wn(:,1));

% 时间推进循环（从第2个时间步到倒数第1个，计算下一个时间步）
for j = 2 : number_of_discrete_time_steps - 1  
    % 空间内部点计算（边界点这里假设为解析解边界值，也可按需设置边界条件）
    for i = 2 : number_of_discrete_space_steps - 1 
        % 从Klein-Gordon方程离散得到的显式格式：W(i,j+1) = 2W(i,j) - W(i,j-1) + (dt^2/a^2)*(a^2*(W(i+1,j)-2W(i,j)+W(i-1,j))/dx^2 + b*W(i,j))
        % 简化后：
        Wn(i, j+1) = 2*Wn(i,j) - Wn(i,j-1) + (dt^2)*( (Wn(i+1,j)-2*Wn(i,j)+Wn(i-1,j))/dx^2 + (b/a^2)*Wn(i,j) );
    end
    % 边界点直接用解析解（保证边界条件一致）
    Wn(1, j+1) = Wa(1, j+1);
    Wn(end, j+1) = Wa(end, j+1);
end

%% 数值解绘制
figure('Name', '数值解'); 
surface(t, x, Wn, 'edgecolor', 'none'); 
colormap(jet(256)); 
colorbar; 
xlabel('t'); 
ylabel('x'); 
title('Wn(x, t) - 数值解');

%% 误差分析（可选）
figure('Name', '数值解与解析解误差');
surface(t, x, abs(Wn - Wa), 'edgecolor', 'none');
colormap(jet(256));
colorbar;
xlabel('t');
ylabel('x');
title('|Wn - Wa| - 误差分布');

关键修正说明 #

1. 初始条件修正：
   • 初始时刻$t=0$：$W(x,0) = A\cos(\lambda x)$，对应Wn(:,1)。
   • 第二个时间层$t=dt$：利用泰勒展开结合初始时间导数$W_t(x,0)=B\mu\cos(\lambda x)$，以及方程本身$W_{tt}=a^2W_{xx}+bW$推导得到，确保初始时间推进的正确性。
2. 差分格式修正：
   从Klein-Gordon方程的二阶时间和空间离散出发，推导得到显式推进格式：
   $$W_{i,j+1} = 2W_{i,j} - W_{i,j-1} + \Delta t^2\left( \frac{W_{i+1,j}-2W_{i,j}+W_{i-1,j}}{\Delta x^2} + \frac{b}{a^2}W_{i,j} \right)$$
   确保时间步按顺序推进，每次用已计算的$j-1$和$j$时间层计算$j+1$层。
3. 边界条件处理：直接使用解析解的边界值，保证边界条件与解析解一致，避免边界误差累积。

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