二维有限元方法模拟静电场透镜在matlab中用泊松方程的代码

以下是一个简单的二维静电场透镜模拟的 Matlab 代码，使用有限元方法和泊松方程来解决问题：

% 定义透镜的参数
r1 = 2; % 内径
r2 = 4; % 外径
d = 1; % 厚度
n = 2; % 折射率

% 定义模型的参数和网格
L = 10; % x, y的长度
nx = 50; % x方向网格数量
ny = 50; % y方向网格数量
x = linspace(-L/2, L/2, nx);
y = linspace(-L/2, L/2, ny);
[X,Y] = meshgrid(x,y);
h = x(2) - x(1); % 网格大小

% 定义边界条件
bc = zeros(nx,ny);
bc(1,:) = 1; % 左侧
bc(end,:) = 1; % 右侧
bc(:,1) = 1; % 底部
bc(:,end) = 1; % 顶部

% 初始化解向量
u = zeros(nx,ny);

% 循环求解泊松方程
for k = 1:1000
    % 计算更新后的解
    u_new = zeros(nx,ny);
    for i = 2:nx-1
        for j = 2:ny-1
            if (i >= nx/2-r2/h && i <= nx/2-r1/h && sqrt((j-ny/2)^2+(i-nx/2)^2) >= r1/h && sqrt((j-ny/2)^2+(i-nx/2)^2) <= r2/h)
                % 透镜区域内的节点
                u_new(i,j) = ((u(i+1,j)+u(i-1,j))/h^2 + (u(i,j+1)+u(i,j-1))/h^2 - n^2/h^2) / (2/h^2);
            else
                % 非透镜区域内的节点
                u_new(i,j) = ((u(i+1,j)+u(i-1,j))/h^2 + (u(i,j+1)+u(i,j-1))/h^2 - bc(i,j)*n^2/h^2) / (2/h^2);
            end
        end
    end
    
    % 判断收敛性
    if max(max(abs(u_new-u))) < 1e-4
        break;
    end
    
    % 更新解向量
    u = u_new;
end

% 绘制结果
figure;
surf(X,Y,u');

这个代码中，我们首先定义了透镜的几何参数，然后定义了模型的网格和边界条件。然后我们使用一个循环来求解泊松方程，并在每次循环中更新解向量。最后，我们绘制了透镜的电势分布。