SLX=ZK(ii,1); SLY=ZK(ii,2); rr=sqrt((SLX-JSX)^2+(SLY-JSY)^2); Xt=JSX-SLX; %计算点到钻孔的x距离 for j=1:nj if (j==1) z=1; elseif(j==nj) z=H-1; else z=(j-1)*dz; end for k=1:nj if(k==1) a=0; b=dz/2; elseif(k==nj) a=H-dz/2; b=H; else a=(2*k-3)*dz*0.5; b=(2*k-1)*dz*0.5; end rydis=(a+b)/2; jsdis=z; [v,Rap,Iap,Rlamd,Ilamd] = untitled55(rydis,jsdis,rr); end end end end end向量化代码

可以将代码中的循环替换为向量化运算，如下所示：

SLX = ZK(:,1);
SLY = ZK(:,2);
rr = sqrt((SLX-JSX).^2 + (SLY-JSY).^2);
Xt = JSX - SLX;

z = (1:nj).' * dz; % 生成 z 向量
z(1) = 1; % 第一个元素赋值为 1
z(nj) = H - 1; % 最后一个元素赋值为 H-1
a = (2.*(1:nj)-3) .* dz./2; % 生成 a 向量
b = (2.*(1:nj)-1) .* dz./2; % 生成 b 向量
a(1) = 0; % 第一个元素赋值为 0
b(nj) = H; % 最后一个元素赋值为 H

[k, j] = meshgrid(1:nj, 1:nj); % 生成 k 和 j 矩阵
rydis = (a+b)/2;
jsdis = z(j);
[v, Rap, Iap, Rlamd, Ilamd] = untitled55(rydis, jsdis, rr);

Tj = zeros(LL, nj, nj, length(ZK), 2);
t = (1:LL) * dt; % 生成 t 向量

for ii = 1:length(ZK)
    aa = integral(@(x)0.25*exp(v(:,:,ii).*Xt(ii).*0.5./Rap(:,:,ii)).*exp(-v(:,:,ii).*sqrt(Rap(:,:,ii).^2 + (z - x).*(z - x)).*0.5./Rap(:,:,ii)).*erfc((sqrt(Rap(:,:,ii).^2 + (z - x).*(z - x)) - v(:,:,ii).*t.')*0.5./sqrt(Rap(:,:,ii).*t.'))...
        ./sqrt(Rap(:,:,ii).^2 + (z - x).*(z - x))./(2.*3.1415926.*Rlamd(:,:,ii)), a, b);

    ab = integral(@(x)0.25*exp(v(:,:,ii).*Xt(ii).*0.5./Rap(:,:,ii)).*exp(v(:,:,ii).*sqrt(Rap(:,:,ii).^2 + (z - x).*(z - x)).*0.5./Rap(:,:,ii)).*erfc((sqrt(Rap(:,:,ii).^2 + (z - x).*(z - x)) + v(:,:,ii).*t.')*0.5./sqrt(Rap(:,:,ii).*t.'))...
        ./sqrt(Rap(:,:,ii).^2 + (z - x).*(z - x))./(2.*3.1415926.*Rlamd(:,:,ii)), a, b);

    ac = integral(@(x)0.25*exp(v(:,:,ii).*Xt(ii).*0.5./Iap(:,:,ii)).*exp(-v(:,:,ii).*sqrt(Iap(:,:,ii).^2 + (z + x).*(z + x)).*0.5./Iap(:,:,ii)).*erfc((sqrt(Iap(:,:,ii).^2 + (z + x).*(z + x)) - v(:,:,ii).*t.')*0.5./sqrt(Iap(:,:,ii).*t.'))...
        ./sqrt(Iap(:,:,ii).^2 + (z + x).*(z + x))./(2.*3.1415926.*Ilamd(:,:,ii)), a, b);

    ad = integral(@(x)0.25*exp(v(:,:,ii).*Xt(ii).*0.5./Iap(:,:,ii)).*exp(v(:,:,ii).*sqrt(Iap(:,:,ii).^2 + (z + x).*(z + x)).*0.5./Iap(:,:,ii)).*erfc((sqrt(Iap(:,:,ii).^2 + (z + x).*(z + x)) + v(:,:,ii).*t.')*0.5./sqrt(Iap(:,:,ii).*t.'))...
        ./sqrt(Iap(:,:,ii).^2 + (z + x).*(z + x))./(2.*3.1415926.*Ilamd(:,:,ii)), a, b);

    aa(isnan(aa)) = 0;
    ab(isnan(ab)) = 0;
    ac(isnan(ac)) = 0;
    ad(isnan(ad)) = 0;

    Tj(:, :, :, ii, 1) = aa + ab - ac - ad;
    Tj(:, :, :, ii, 2) = repmat(rydis.', LL, 1);
end