基于有限差分法的二维边值问题数值分析
有限差分法是求解偏微分方程边值问题最常用的数值方法之一。
理论基础
问题描述
考虑二维泊松方程的一般形式：
-∇²u = f(x,y), 在区域Ω内
u = g(x,y), 在边界∂Ω上
其中∇²是拉普拉斯算子，f是源项，g是边界条件。
有限差分法基本原理
有限差分法的核心思想是用差商代替微商，将连续的偏微分方程转化为离散的代数方程组。
拉普拉斯算子的离散化：
∇²u ≈ [u(x+h,y) + u(x-h,y) + u(x,y+h) + u(x,y-h) - 4u(x,y)] / h²
数值方法实现
五点差分格式
对于均匀网格，最常用的是五点差分格式：
function [U, X, Y] = finiteDifference2D(f, g, a, b, c, d, nx, ny)
% 二维边值问题的有限差分法求解
% 输入参数：
% f: 源项函数句柄 f(x,y)
% g: 边界条件函数句柄 g(x,y)
% a,b: x方向区间 [a,b]
% c,d: y方向区间 [c,d]
% nx,ny: x和y方向的网格点数
% 网格参数
hx = (b-a)/(nx-1);
hy = (d-c)/(ny-1);
% 网格坐标
X = linspace(a, b, nx);
Y = linspace(c, d, ny);
[Xm, Ym] = meshgrid(X, Y);
% 初始化解矩阵
U = zeros(ny, nx);
% 应用边界条件
for i = 1:nx
U(1,i) = g(X(i), Y(1)); % 下边界
U(ny,i) = g(X(i), Y(ny)); % 上边界
end
for j = 1:ny
U(j,1) = g(X(1), Y(j)); % 左边界
U(j,nx) = g(X(nx), Y(j)); % 右边界
end
% 构造系数矩阵和右端项（使用向量化方法）
N = (nx-2)*(ny-2); % 内部节点数
A = sparse(N, N);
F = zeros(N, 1);
% 网格参数
rx = 1/hx^2;
ry = 1/hy^2;
r = 2*(rx + ry);
% 填充系数矩阵
for j = 2:ny-1
for i = 2:nx-1
% 当前节点的全局索引
k = (j-2)*(nx-2) + (i-1);
% 主对角线
A(k,k) = r;
% x方向邻居
if i > 2
A(k,k-1) = -rx; % 左邻居
end
if i < nx-1
A(k,k+1) = -rx; % 右邻居
end
% y方向邻居
if j > 2
A(k,k-(nx-2)) = -ry; % 下邻居
end
if j < ny-1
A(k,k+(nx-2)) = -ry; % 上邻居
end
% 右端项（考虑边界贡献）
F(k) = f(X(i), Y(j));
% 边界贡献
if i == 2
F(k) = F(k) + rx * g(X(1), Y(j));
end
if i == nx-1
F(k) = F(k) + rx * g(X(nx), Y(j));
end
if j == 2
F(k) = F(k) + ry * g(X(i), Y(1));
end
if j == ny-1
F(k) = F(k) + ry * g(X(i), Y(ny));
end
end
end
% 求解线性方程组
u_inner = A \ F;
% 将内部解填充到U矩阵
for j = 2:ny-1
for i = 2:nx-1
k = (j-2)*(nx-2) + (i-1);
U(j,i) = u_inner(k);
end
end
end
逐次超松弛迭代法(SOR)
对于大规模问题，迭代法更有效：
function U = SOR_method(f, g, a, b, c, d, nx, ny, omega, max_iter, tol)
% SOR方法求解二维泊松方程
% omega: 松弛因子 (1 < omega < 2)
% 网格参数
hx = (b-a)/(nx-1);
hy = (d-c)/(ny-1);
X = linspace(a, b, nx);
Y = linspace(c, d, ny);
% 初始化
U = zeros(ny, nx);
% 边界条件
for i = 1:nx
U(1,i) = g(X(i), Y(1));
U(ny,i) = g(X(i), Y(ny));
end
for j = 1:ny
U(j,1) = g(X(1), Y(j));
U(j,nx) = g(X(nx), Y(j));
end
% 迭代参数
rx = 1/hx^2;
ry = 1/hy^2;
r = 2*(rx + ry);
% SOR迭代
for iter = 1:max_iter
max_err = 0;
for j = 2:ny-1
for i = 2:nx-1
% 五点差分格式
residual = (rx*(U(j,i-1) + U(j,i+1)) + ...
ry*(U(j-1,i) + U(j+1,i)) - ...
f(X(i), Y(j))) / r - U(j,i);
% SOR更新
U(j,i) = U(j,i) + omega * residual;
% 更新最大误差
max_err = max(max_err, abs(residual));
end
end
% 检查收敛
if max_err < tol
fprintf('SOR收敛于第%d次迭代，误差: %e\n', iter, max_err);
break;
end
if iter == max_iter
fprintf('达到最大迭代次数%d，最终误差: %e\n', iter, max_err);
end
end
end
完整仿真示例
示例1：简单泊松方程
% 示例1：-∇²u = 2π²sin(πx)sin(πy)，边界为0
% 真解：u(x,y) = sin(πx)sin(πy)
% 定义问题参数
a = 0; b = 1; c = 0; d = 1;
nx = 51; ny = 51;
% 定义函数
f = @(x,y) 2*pi^2 * sin(pi*x) .* sin(pi*y);
g = @(x,y) zeros(size(x)); % 齐次边界条件
exact_solution = @(x,y) sin(pi*x) .* sin(pi*y);
% 有限差分求解
[U, X, Y] = finiteDifference2D(f, g, a, b, c, d, nx, ny);
% 计算精确解
U_exact = exact_solution(X, Y);
% 计算误差
error = abs(U - U_exact);
max_error = max(error(:));
rms_error = sqrt(mean(error(:).^2));
fprintf('最大误差: %.6e\n', max_error);
fprintf('RMS误差: %.6e\n', rms_error);
% 可视化结果
figure('Position', [100, 100, 1200, 400]);
subplot(1,3,1);
surf(X, Y, U);
title('数值解');
xlabel('x'); ylabel('y'); zlabel('u(x,y)');
subplot(1,3,2);
surf(X, Y, U_exact);
title('精确解');
xlabel('x'); ylabel('y'); zlabel('u(x,y)');
subplot(1,3,3);
surf(X, Y, error);
title('绝对误差');
xlabel('x'); ylabel('y'); zlabel('误差');
示例2：非齐次边界条件
% 示例2：-∇²u = 0，非齐次边界条件
% 真解：u(x,y) = x² - y²
% 定义问题参数
a = -1; b = 1; c = -1; d = 1;
nx = 41; ny = 41;
% 定义函数
f = @(x,y) zeros(size(x)); % 拉普拉斯方程
g = @(x,y) x.^2 - y.^2; % 边界条件
exact_solution = @(x,y) x.^2 - y.^2;
% 使用SOR方法求解
omega = 1.8; % 最优松弛因子
max_iter = 1000;
tol = 1e-6;
U_sor = SOR_method(f, g, a, b, c, d, nx, ny, omega, max_iter, tol);
% 计算精确解和误差
[X, Y] = meshgrid(linspace(a,b,nx), linspace(c,d,ny));
U_exact = exact_solution(X, Y);
error_sor = abs(U_sor - U_exact);
fprintf('SOR方法最大误差: %.6e\n', max(error_sor(:)));
% 收敛性分析
figure('Position', [100, 100, 1000, 400]);
subplot(1,2,1);
contourf(X, Y, U_sor, 20);
title('数值解等值线');
xlabel('x'); ylabel('y');
colorbar;
axis equal;
subplot(1,2,2);
contourf(X, Y, error_sor, 20);
title('误差分布');
xlabel('x'); ylabel('y');
colorbar;
axis equal;
示例3：网格收敛性分析
% 网格收敛性分析
% 验证有限差分法的收敛阶
a = 0; b = 1; c = 0; d = 1;
f = @(x,y) 2*pi^2 * sin(pi*x) .* sin(pi*y);
g = @(x,y) zeros(size(x));
exact_solution = @(x,y) sin(pi*x) .* sin(pi*y);
% 不同网格尺寸
grid_sizes = [11, 21, 41, 81];
errors = zeros(size(grid_sizes));
h_values = zeros(size(grid_sizes));
for i = 1:length(grid_sizes)
nx = grid_sizes(i);
ny = grid_sizes(i);
h = 1/(nx-1);
h_values(i) = h;
[U, X, Y] = finiteDifference2D(f, g, a, b, c, d, nx, ny);
U_exact = exact_solution(X, Y);
errors(i) = max(abs(U(:) - U_exact(:)));
end
% 计算收敛阶
orders = zeros(length(grid_sizes)-1, 1);
for i = 2:length(grid_sizes)
orders(i-1) = log(errors(i-1)/errors(i)) / log(h_values(i-1)/h_values(i));
end
% 绘制收敛图
figure;
loglog(h_values, errors, 'o-', 'LineWidth', 2, 'MarkerSize', 8);
hold on;
loglog(h_values, h_values.^2, '--', 'LineWidth', 1.5);
xlabel('网格尺寸 h');
ylabel('最大误差');
title('有限差分法收敛性分析');
legend('计算误差', 'O(h^2)参考线', 'Location', 'northwest');
grid on;
fprintf('收敛阶估计:\n');
for i = 1:length(orders)
fprintf('从h=%.3f到h=%.3f: %.4f\n', h_values(i), h_values(i+1), orders(i));
end
高级应用
非均匀网格处理
function U = nonuniformFD(f, g, x, y)
% 非均匀网格的有限差分法
% x, y: 非均匀网格坐标向量
nx = length(x);
ny = length(y);
U = zeros(ny, nx);
% 应用边界条件
for i = 1:nx
U(1,i) = g(x(i), y(1));
U(ny,i) = g(x(i), y(ny));
end
for j = 1:ny
U(j,1) = g(x(1), y(j));
U(j,nx) = g(x(nx), y(j));
end
% 构造系数矩阵（考虑非均匀网格间距）
N = (nx-2)*(ny-2);
A = sparse(N, N);
F = zeros(N, 1);
for j = 2:ny-1
for i = 2:nx-1
k = (j-2)*(nx-2) + (i-1);
% x方向间距
hx_left = x(i) - x(i-1);
hx_right = x(i+1) - x(i);
% y方向间距
hy_bottom = y(j) - y(j-1);
hy_top = y(j+1) - y(j);
% 非均匀网格系数
cx_left = 2/(hx_left*(hx_left + hx_right));
cx_right = 2/(hx_right*(hx_left + hx_right));
cx_center = -2/(hx_left*hx_right);
cy_bottom = 2/(hy_bottom*(hy_bottom + hy_top));
cy_top = 2/(hy_top*(hy_bottom + hy_top));
cy_center = -2/(hy_bottom*hy_top);
% 填充矩阵
A(k,k) = cx_center + cy_center;
if i > 2
A(k,k-1) = cx_left;
end
if i < nx-1
A(k,k+1) = cx_right;
end
if j > 2
A(k,k-(nx-2)) = cy_bottom;
end
if j < ny-1
A(k,k+(nx-2)) = cy_top;
end
F(k) = f(x(i), y(j));
end
end
% 求解
u_inner = A \ F;
% 填充解
for j = 2:ny-1
for i = 2:nx-1
k = (j-2)*(nx-2) + (i-1);
U(j,i) = u_inner(k);
end
end
end
参考代码 基于有限差分法的二维边值问题的数值分析 www.youwenfan.com/contentcnl/79244.html
结果分析与验证
精度验证
- 五点差分格式具有二阶精度 O(h²)
- 误差主要来源于截断误差
- 收敛性可通过网格细化验证
稳定性考虑
- 直接法（如反斜杠运算）对小规模问题稳定
- 迭代法（如SOR）需要选择合适的松弛因子
- 条件数随网格细化而增大
计算效率
- 直接法：O(N³) 时间复杂度，适合中小规模问题
- 迭代法：O(N²) 每次迭代，适合大规模问题
- 稀疏矩阵存储可显著减少内存需求
浙公网安备 33010602011771号