当前位置：首页 > news >正文

基于有限差分法的二维边值问题数值分析

news 2025/11/5 7:44:07

有限差分法是求解偏微分方程边值问题最常用的数值方法之一。

理论基础

问题描述

考虑二维泊松方程的一般形式：

-∇²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:nxU(1,i) = g(X(i), Y(1));    % 下边界U(ny,i) = g(X(i), Y(ny));  % 上边界endfor j = 1:nyU(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-1for i = 2:nx-1% 当前节点的全局索引k = (j-2)*(nx-2) + (i-1);% 主对角线A(k,k) = r;% x方向邻居if i > 2A(k,k-1) = -rx;  % 左邻居endif i < nx-1A(k,k+1) = -rx;  % 右邻居end% y方向邻居if j > 2A(k,k-(nx-2)) = -ry;  % 下邻居endif j < ny-1A(k,k+(nx-2)) = -ry;  % 上邻居end% 右端项（考虑边界贡献）F(k) = f(X(i), Y(j));% 边界贡献if i == 2F(k) = F(k) + rx * g(X(1), Y(j));endif i == nx-1F(k) = F(k) + rx * g(X(nx), Y(j));endif j == 2F(k) = F(k) + ry * g(X(i), Y(1));endif j == ny-1F(k) = F(k) + ry * g(X(i), Y(ny));endendend% 求解线性方程组u_inner = A \ F;% 将内部解填充到U矩阵for j = 2:ny-1for i = 2:nx-1k = (j-2)*(nx-2) + (i-1);U(j,i) = u_inner(k);endend
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:nxU(1,i) = g(X(i), Y(1));U(ny,i) = g(X(i), Y(ny));endfor j = 1:nyU(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_itermax_err = 0;for j = 2:ny-1for 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));endend% 检查收敛if max_err < tolfprintf('SOR收敛于第%d次迭代，误差: %e\n', iter, max_err);break;endif iter == max_iterfprintf('达到最大迭代次数%d，最终误差: %e\n', iter, max_err);endend
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:nxU(1,i) = g(x(i), y(1));U(ny,i) = g(x(i), y(ny));endfor j = 1:nyU(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-1for i = 2:nx-1k = (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 > 2A(k,k-1) = cx_left;endif i < nx-1A(k,k+1) = cx_right;endif j > 2A(k,k-(nx-2)) = cy_bottom;endif j < ny-1A(k,k+(nx-2)) = cy_top;endF(k) = f(x(i), y(j));endend% 求解u_inner = A \ F;% 填充解for j = 2:ny-1for i = 2:nx-1k = (j-2)*(nx-2) + (i-1);U(j,i) = u_inner(k);endend
end