MATLAB数值积分求立方体电势与可视化

空间中有一均匀带电的立方体，电荷密度 ρ = 1，边长为 a = 1。 使用数值积分求空间中各点的电场和电势

主程序（根据输入坐标值求该点电势，场强）：

% 开始计时
tic;    

% 接受用户输入的坐标
x_input = input('请输入 x 坐标：');
y_input = input('请输入 y 坐标：');
z_input = input('请输入 z 坐标：');

% 调用函数计算电场和电势
[E_field, V_potential] = caculate_E(x_input, y_input, z_input);


 % 停止计时并获取经过的时间
 time_elapsed = toc;

% 输出结果
disp(['Electric Field at (', num2str(x_input), ', ', num2str(y_input), ', ', num2str(z_input), '): ', num2str(E_field)]);
disp(['Electric Potential at (', num2str(x_input), ', ', num2str(y_input), ', ', num2str(z_input), '): ', num2str(V_potential)]);
disp(['Time Elapsed: ', num2str(time_elapsed), ' seconds']);

计算电势函数：

function [E_x, E_y, E_z] = caculate_E(x_input, y_input, z_input)
    % 该函数计算给定坐标处的电场分量
    
    % 定义常数
    epsilon_0 = 8.854e-12; % 真空中的介电常数，单位：法拉/米
    a = 1; % 立方体的边长，假设为1
    rho = 1; % 电荷密度，假设为1
    
    % 定义电场的积分函数
    % x 方向的电场积分函数
    integrand_E_x = @(x0, y0, z0) rho .* (x_input - x0) ./ ((x_input - x0).^2 + (y_input - y0).^2 + (z_input - z0).^2).^(3/2);
    % y 方向的电场积分函数
    integrand_E_y = @(x0, y0, z0) rho .* (y_input - y0) ./ ((x_input - x0).^2 + (y_input - y0).^2 + (z_input - z0).^2).^(3/2);
    % z 方向的电场积分函数
    integrand_E_z = @(x0, y0, z0) rho .* (z_input - z0) ./ ((x_input - x0).^2 + (y_input - y0).^2 + (z_input - z0).^2).^(3/2);
    
    % 计算电场分量
    % 使用 integral3 函数对三维区域进行积分，分别计算 x、y、z 方向的电场分量
    % integral3 参数依次为积分函数、x 范围、y 范围、z 范围、以及积分的精度选项
    E_x = integral3(integrand_E_x, -a/2, a/2, -a/2, a/2, -a/2, a/2, 'AbsTol', 1e-9, 'RelTol', 1e-9);
    E_y = integral3(integrand_E_y, -a/2, a/2, -a/2, a/2, -a/2, a/2, 'AbsTol', 1e-9, 'RelTol', 1e-9);
    E_z = integral3(integrand_E_z, -a/2, a/2, -a/2, a/2, -a/2, a/2, 'AbsTol', 1e-9, 'RelTol', 1e-9);
    
    % 应用库仑定律的常数
    k = 1 / (4 * pi * epsilon_0);
    
    % 计算电场
    % 应用比例常数，计算最终的电场分量
    E_x = k * E_x; % x 方向的电场分量
    E_y = k * E_y; % y 方向的电场分量
    E_z = k * E_z; % z 方向的电场分量
end

 可视化

% 定义空间范围和点数量
x_range = linspace(1, 2, 10); % x 范围从 1 到 2，10 个点
y_range = linspace(1, 2, 10); % y 范围从 1 到 2，10 个点
z_range = linspace(1, 2, 10); % z 范围从 1 到 2，10 个点

% 创建均匀网格
[x_points, y_points, z_points] = meshgrid(x_range, y_range, z_range);

% 将网格展平成一维数组
x_points = x_points(:);
y_points = y_points(:);
z_points = z_points(:);

% 预分配存储电场和电势值的数组
E_field_values = zeros(length(x_points), 3);
V_potential_values = zeros(length(x_points), 1);

% 计算电场和电势值
for i = 1:length(x_points)
    [E_field_values(i, :), V_potential_values(i)] = caculate(x_points(i), y_points(i), z_points(i));
end

% 创建三维散点图 - 电势值
figure;
scatter3(x_points, y_points, z_points, 50, V_potential_values, 'filled');
colorbar;
xlabel('X');
ylabel('Y');
zlabel('Z');
title('Potential Distribution (3D Scatter)');
colormap('jet'); % 使用 jet 颜色映射增强对比度

% 创建存储电场值的数组
E_x_values = zeros(length(x_points), 1);
E_y_values = zeros(length(x_points), 1);
E_z_values = zeros(length(x_points), 1);

% 计算电场值
for i = 1:length(x_points)
    [E_x_values(i), E_y_values(i), E_z_values(i)] = caculate_E(x_points(i), y_points(i), z_points(i));
end

% 创建三维散点图 - 电场值
figure;
quiver3(x_points, y_points, z_points, E_x_values, E_y_values, E_z_values);
xlabel('X');
ylabel('Y');
zlabel('Z');
title('Electric Field Distribution (3D Quiver)');

结果(第一卦限)：