matlab仿真实验-静电场

题目：用有限差分法求下图区域中各个节点的电位。试用MATLAB编写文件，计算矩形区域内电位分布

（要求：不含边界的节点数目100*100，边长设为10cm；差分方程组的计算采用超松弛迭代方法；选择适当的松弛因子，提高迭代效率。）

为了完成这个题目，我们可以按照以下仿真步骤来使用MATLAB进行计算和分析：

仿真步骤

1. 定义问题和参数：

   • 创建一个矩形区域，其尺寸为 10cm x 10cm。
   • 不含边界的节点数目为 100x100，这意味着整个网格（包含边界）的尺寸将是 102x102。

2. 初始化电位矩阵：

   • 初始化一个 102x102 的矩阵来表示电位分布。初始时，所有节点的电位值设为零。
   • 设置边界条件。例如，可以将顶部边界设为100V，其他边界设为0V。

3. 设置松弛因子：

   选择一个合适的松弛因子（alpha），它用于加速迭代过程。通常这个值需要通过试验来确定，但一般在 1.0 到 2.0 之间。

4. 应用有限差分法：

   使用有限差分法对拉普拉斯方程进行离散化。在这个情况下，每个内部节点的电位值由其邻近四个节点的电位值决定。

5. 超松弛迭代方法：

   使用超松弛迭代方法（SOR）来求解电位分布。在每次迭代中，更新每个内部节点的电位值。持续迭代直到达到预设的收敛条件，例如，当任何节点的电位变化小于某个阈值时（比如 1e-5）。

6. 迭代过程：

   在迭代过程中，跟踪并记录所需的迭代次数。这有助于分析不同松弛因子对迭代效率的影响。

7. 结果分析与可视化：

   • 绘制电位分布的三维图和等高线图，以直观地展示电位在整个矩形区域内的分布情况。
   • 分析不同松弛因子下的迭代效率，以确定最优的松弛因子。

仿真代码：

% 清空环境变量
clc; clear;

% 存储不同alpha的迭代次数
results = ones(150, 1);

% 步长0.01寻找松弛因子最优解
for alpha = 1.01:0.01:1.99
    % 设置长宽节点数
    hx = 101; hy = 101;
    
    % 用0初始化电位矩阵，并设置边界条件
    v1 = zeros(hy, hx);
    v1(hy, :) = zeros(1, hx);  % 底部边界设为0V
    v1(1, :) = ones(1, hx) * 100;  % 顶部边界设为100V
    v1(:, 1) = zeros(hy, 1);  % 左侧边界设为0V
    v1(:, hx) = zeros(hy, 1);  % 右侧边界设为0V

    % 初始化迭代过程中的电位矩阵
    v2 = v1;
    
    % 迭代控制变量
    maxt = 1;
    iteration = 0; % 本次迭代次数
    
    % 设置精度跳出迭代
    while(maxt > 1e-5)
        iteration = iteration + 1;
        maxt = 0;
        
        % 对内部节点进行迭代
        for i = 2:hy-1
            for j = 2:hx-1
                temp = v2(i, j);
                % 应用超松弛迭代公式
                v2(i, j) = v1(i, j) + ((v1(i, j+1) + v1(i+1, j) + v2(i-1, j) + v2(i, j-1) - 4 * v1(i, j)) * alpha / 4);
                t = abs(v2(i, j) - temp);
                if(t > maxt)
                    maxt = t;
                end
            end
        end
        
        % 更新电位矩阵
        v1 = v2;
    end
    
    % 存储结果
    results(int16(alpha * 100) - 100) = iteration;
end

% 绘制电位分布图形
clf;
subplot(1, 2, 1), mesh(v2);
axis([0, 101, 0, 101, 0, 100]);
title('电位分布三维图');

subplot(1, 2, 2), contour(v2, 15);
title('电位分布等高线图');
hold on;
axis([-1, 102, -1, 110]);
plot([1, 1, hx, hx, 1], [1, hy+1, hy+1, 1, 1], 'r');
text(hx/2, 0.3, '0V', 'fontsize', 11);
text(hx/2 - 0.5, hy + 0.5, '100V', 'fontsize', 11);
text(-0.5, hy/2, '0V', 'fontsize', 11);
text(hx + 0.3, hy/2, '0V', 'fontsize', 11);
hold off;

 运行结果