这一段代码中的循环嵌套比较复杂，但也可以通过向量化的方式进行改写。具体实现如下： 复制 [a1, b1] = meshgrid(1:jx, 1:jy); JSX = a1*gj-r0; JSY = b1*gj; SLX = ZK(:,1); SLY = ZK(:,2); SLX = repmat(SLX, 1, jx*jy); SLY = repmat(SLY, 1, jx*jy); r = sqrt((SLX-JSX(:)').^2 + (SLY-JSY(:)').^2); rr = reshape(r, [], 1); Xt = repmat(JSX(:)'-SLX, nj*LL, 1); z = (0:nj)*dz; z = repmat(z', 1, jx*jy); z = reshape(z, [], 1); a = (2*(1:nj)-3)*dz*0.5; b = (2*(1:nj)-1)*dz*0.5; a = repmat(a', 1, jx*jy); b = repmat(b', 1, jx*jy); a = reshape(a, [], 1); b = reshape(b, [], 1); Tj = zeros(LL, jx*jy); for ii = 1:zks for j = 1:nj if (j == 1) zj = 1; elseif (j == nj) zj = H-1; else zj = (j-1)*dz; end zj = repmat(zj, 1, jx*jy); Z1 = rr.^2 + (z-zj).^2; Z2 = rr.^2 + (z+zj).^2; aa = trapz(x, 0.25*exp(v*Xt*0.5 / Rap).*exp(-v * sqrt(Z1)*0.5 / Rap).*erfc((sqrt(Z1) - v * t).*0.5 ./ sqrt(Rap*t))./sqrt(Z1))/(2 * 3.1415926*Rlamd); ab = trapz(x, 0.25*exp(v*Xt*0.5 / Rap).*exp(v*sqrt(Z1)*0.5 / Rap).*erfc((sqrt(Z1) + v * t).*0.5 ./ sqrt(Rap*t))./sqrt(Z1))/(2 * 3.1415926*Rlamd); ac = trapz(x, 0.25*exp(v*Xt*0.5 / Iap).*exp(-v * sqrt(Z2)*0.5 / Iap).*erfc((sqrt(Z2) - v * t).*0.5 ./ sqrt(Iap*t))./sqrt(Z2))/(2 * 3.1415926*Ilamd); ad = trapz(x, 0.25*exp(v*Xt*0.5 / Iap).*exp(v*sqrt(Z2)*0.5 / Iap).*erfc((sqrt(Z2) + v * t).*0.5 ./ sqrt(Iap*t))./sqrt(Z2))/(2 * 3.1415926*Ilamd); Tj(:, :) = Tj(:, :) + (aa + ab - ac - ad); end end Tj = reshape(Tj, LL, jx, jy);向量化这段代码

可以将循环嵌套的部分替换为矩阵运算，从而实现向量化。具体实现如下：

[a1, b1] = meshgrid(1:jx, 1:jy);
JSX = a1*gj-r0;
JSY = b1*gj;
SLX = ZK(:,1);
SLY = ZK(:,2);
SLX = repmat(SLX, 1, jx*jy);
SLY = repmat(SLY, 1, jx*jy);
r = sqrt((SLX-JSX(:)').^2 + (SLY-JSY(:)').^2);
rr = reshape(r, [], 1);
Xt = repmat(JSX(:)'-SLX, nj*LL, 1);
z = (0:nj)*dz;
z = repmat(z', 1, jx*jy);
z = reshape(z, [], 1);
a = (2*(1:nj)-3)*dz*0.5;
b = (2*(1:nj)-1)*dz*0.5;
a = repmat(a', 1, jx*jy);
b = repmat(b', 1, jx*jy);
a = reshape(a, [], 1);
b = reshape(b, [], 1);

ZJ = repmat((1:nj)', 1, jx*jy);
ZJ(ZJ == 1) = 0;
ZJ(ZJ == nj) = H-1;
zj = ZJ*dz;
zj = repmat(zj, LL, 1);

Z1 = rr.^2 + (z-zj(:)').^2;
Z2 = rr.^2 + (z+zj(:)').^2;
exp1 = exp(v*Xt*0.5 / Rap);
exp2 = exp(v*Xt*0.5 / Iap);
sqrt1 = sqrt(Z1)*0.5 / Rap;
sqrt2 = sqrt(Z2)*0.5 / Iap;
erfc1 = erfc((sqrt(Z1) - v * t).*0.5 ./ sqrt(Rap*t));
erfc2 = erfc((sqrt(Z1) + v * t).*0.5 ./ sqrt(Rap*t));
erfc3 = erfc((sqrt(Z2) - v * t).*0.5 ./ sqrt(Iap*t));
erfc4 = erfc((sqrt(Z2) + v * t).*0.5 ./ sqrt(Iap*t));
aa = trapz(x, 0.25*exp1.*exp(-v*sqrt1).*erfc1./sqrt(Z1))/(2 * 3.1415926*Rlamd);
ab = trapz(x, 0.25*exp1.*exp(v*sqrt1).*erfc2./sqrt(Z1))/(2 * 3.1415926*Rlamd);
ac = trapz(x, 0.25*exp2.*exp(-v*sqrt2).*erfc3./sqrt(Z2))/(2 * 3.1415926*Ilamd);
ad = trapz(x, 0.25*exp2.*exp(v*sqrt2).*erfc4./sqrt(Z2))/(2 * 3.1415926*Ilamd);
Tj = reshape(aa + ab - ac - ad, LL, jx*jy);
Tj = sum(Tj, 2);
Tj = reshape(Tj, LL, jx, jy);

这样，循环嵌套的部分就被替换为矩阵运算，可以大大提高代码的效率。