[v,Rap,Iap,Rlamd,Ilamd] = untitled55(rydis,jsdis); for i=1:LL t=i*dt; aa=integral(@(x)0.25*exp(v*Xt*0.5 / Rap)*exp(-v * sqrt(rr*rr + (z - x).*(z - x))*0.5 / Rap).*erfc((sqrt(rr*rr + (z - x).*(z - x)) - v * t)*0.5 / sqrt(Rap*t))./sqrt(rr*rr + (z - x).*(z - x))/(2 * 3.1415926*Rlamd),a,b); ab=integral(@(x)0.25*exp(v*Xt*0.5 / Rap)*exp(v*sqrt(rr*rr + (z - x).*(z - x))*0.5 / Rap).*erfc((sqrt(rr*rr + (z - x).*(z - x)) + v * t)*0.5 / sqrt(Rap*t))./sqrt(rr*rr + (z - x).*(z - x))/(2 * 3.1415926*Rlamd),a,b); ac=integral(@(x)0.25*exp(v*Xt*0.5 / Iap)*exp(-v * sqrt(rr*rr + (z + x).*(z + x))*0.5 / Iap).*erfc((sqrt(rr*rr + (z + x).*(z + x)) - v * t)*0.5 / sqrt(Iap*t))./sqrt(rr*rr + (z+ x).*(z + x))/(2 * 3.1415926*Ilamd),a,b); ad=integral(@(x)0.25*exp(v*Xt*0.5 / Iap)*exp(v*sqrt(rr*rr + (z + x).*(z + x))*0.5 / Iap).*erfc((sqrt(rr*rr + (z + x).*(z + x)) + v * t)*0.5 / sqrt(Iap*t))./sqrt(rr*rr + (z + x).*(z + x))/(2 * 3.1415926*Ilamd),a,b); aa(isnan(aa)) = 0;ab(isnan(ab)) = 0;ac(isnan(ac)) = 0; ad(isnan(ad)) = 0; Tj(i,j,k,ii,a1,b1)=(aa+ab-ac-ad); end向量化代码

可以使用 MATLAB 的向量化操作来改写上述代码，如下所示：

[v,Rap,Iap,Rlamd,Ilamd] = untitled55(rydis,jsdis);

x = linspace(a,b,LL);
t = (1:LL)*dt;

rr2 = rr^2;
z2 = z.^2;

aa = zeros(1, LL);
ab = zeros(1, LL);
ac = zeros(1, LL);
ad = zeros(1, LL);
for i = 1:LL
    aa(i) = trapz(x, 0.25*exp(v*Xt*0.5 / Rap).*exp(-v * sqrt(rr2 + (z - x).^2)*0.5 / Rap).*erfc((sqrt(rr2 + (z - x).^2) - v * t(i))*0.5 / sqrt(Rap*t(i)))./sqrt(rr2 + (z - x).^2))/(2 * 3.1415926*Rlamd);
    ab(i) = trapz(x, 0.25*exp(v*Xt*0.5 / Rap).*exp(v*sqrt(rr2 + (z - x).^2)*0.5 / Rap).*erfc((sqrt(rr2 + (z - x).^2) + v * t(i))*0.5 / sqrt(Rap*t(i)))./sqrt(rr2 + (z - x).^2))/(2 * 3.1415926*Rlamd);
    ac(i) = trapz(x, 0.25*exp(v*Xt*0.5 / Iap).*exp(-v * sqrt(rr2 + (z + x).^2)*0.5 / Iap).*erfc((sqrt(rr2 + (z + x).^2) - v * t(i))*0.5 / sqrt(Iap*t(i)))./sqrt(rr2 + (z+ x).^2))/(2 * 3.1415926*Ilamd);
    ad(i) = trapz(x, 0.25*exp(v*Xt*0.5 / Iap).*exp(v*sqrt(rr2 + (z + x).^2)*0.5 / Iap).*erfc((sqrt(rr2 + (z + x).^2) + v * t(i))*0.5 / sqrt(Iap*t(i)))./sqrt(rr2 + (z + x).^2))/(2 * 3.1415926*Ilamd);
end
Tj(:,j,k,ii,a1,b1) = aa + ab - ac - ad;

在这个向量化版本中，我们使用 `linspace` 函数来代替 `for` 循环生成 `x` 向量，使用 `trapz` 函数来代替 `integral` 函数求解积分。这样可以显著提高代码的运行速度。