一开始以为是太复杂，因为我知道代码整的比较粗糙计算比较笨。但是少一个变量的大约6×10×100个类似的式子算得还是比较快的，而且换到工作站用96核上还是跑不动。
另外分这么多变量事实上是另一步定积分的数值计算我给手动拆区间取点取均值了，我目前只是单纯想要相关计算结果，看了下精度和之前耗时还可以接受所以暂时先将就用着在，准备后面再慢慢看怎么做合理。

下面是代码
(*定义 Psi 函数*)\[Psi][theta_, gamma_, delta_] :=
PDF[TruncatedDistribution[{0, Pi/2},
NormalDistribution[gamma, delta]], theta]
(*定义 nz0*)
nz0[gamma_, delta_, y_, lf_] :=
NIntegrate[\[Psi][theta, gamma, delta]*Cos[theta], {theta, 0,
ArcSin[2*y/lf]}]
(*定义 \[CapitalPhi]z0*)
\[CapitalPhi]z0[gamma_, delta_, y_, lf_] :=
NIntegrate[\[Psi][theta, gamma, delta], {theta, 0, ArcSin[2*y/lf]}]
(*定义 A0*)
A0[lf_, theta_] := (Pi*lf^2/2)*Sin[theta]
(*定义 A1*)
A1[lf_, y_, theta_] := lf^2*Sin[theta]*ArcSin[2*y/(lf*Sin[theta])]
(*定义 nz1*)
nz1[gamma_, delta_, y_, z_, lf_] :=
NIntegrate[(A1[lf, y, theta]/A0[lf, theta])*\[Psi][theta, gamma,
delta]*Cos[theta], {theta, ArcSin[2*y/lf], ArcSin[2*z/lf]}]
(*定义 \[CapitalPhi]y1*)
\[CapitalPhi]y1[gamma_, delta_, y_, lf_] :=
NIntegrate[(A1[lf, y, theta]/A0[lf, theta])*\[Psi][theta, gamma,
delta], {theta, ArcSin[2*y/lf], Pi/2}]
(*定义 A2*)
A2[lf_, z_, y_, theta_] :=
Abs[lf^2*Sin[theta]*(ArcSin[2*y/lf] - ArcCos[2*z/lf])]
(*定义 Az*)
Az[lf_, z_, theta_] := lf^2*Sin[theta]*ArcSin[2*z/(lf*Sin[theta])]
(*定义 nz2*)
nz2[gamma_, delta_, y_, z_, lf_] :=
NIntegrate[(A2[lf, y, z, theta]/A0[lf, theta])*\[Psi][theta, gamma,
delta]*Cos[theta], {theta, ArcSin[2*z/lf], Pi/2}]
(*定义 \[CapitalPhi]z2*)
\[CapitalPhi]z2[gamma_, delta_, y_, z_, lf_] :=
NIntegrate[(A2[lf, y, z, theta]/A0[lf, theta])*\[Psi][theta, gamma,
delta], {theta, ArcSin[2*z/lf], Pi/2}]
(*定义 \[CapitalPhi]z3*)
\[CapitalPhi]z3[gamma_, delta_, y_, lf_] :=
NIntegrate[(A0[lf, theta] - A1[lf, y, theta])/A0[lf, theta]*\[Psi][
theta, gamma, delta], {theta, ArcSin[2*y/lf], Pi/2}]
(*定义 \[CapitalPhi]z4*)
\[CapitalPhi]z4[gamma_, delta_, z_, lf_] :=
NIntegrate[(A0[lf, theta] - Az[lf, z, theta])/A0[lf, theta]*\[Psi][
theta, gamma, delta], {theta, ArcSin[2*z/lf], Pi/2}]
(*定义 nz3*)
nz3[gamma_, delta_, y_,
lf_] := \[CapitalPhi]z3[gamma, delta, y,
lf]*(2/Pi)*(Sqrt[lf^2 - 4*y^2]/lf)
(*定义 nz4*)
nz4[gamma_, delta_, y_,
lf_] := \[CapitalPhi]z4[gamma, delta, y,
lf]*(2/Pi)*(Sqrt[lf^2 - 4*y^2]/lf)
(*定义 nyz*)
nyz[gamma_, delta_, y_, z_, lf_] :=
nz1[gamma, delta, y, z, lf] + nz2[gamma, delta, y, z, lf] +
nz3[gamma, delta, y, lf] + nz4[gamma, delta, z, lf]
(*手动输入参数*)
delta = Pi/10;
lf = 13;
(*生成 gamma 的取值列表*)
gammaValues = Table[Pi*i/20, {i, 1, 10}];
(*生成 y 和 z 的取值列表,确保 y 始终大于 z*)
yValues = Table[lf*i/100, {i, 1, 50}]; (*减少取值范围*)
zValues = Table[lf*j/100, {j, 1, 49}]; (*减少取值范围*)
(*生成所有 (y,z) 组合,使得 y>z*)
combinations =
Select[Flatten[Table[{y, z}, {y, yValues}, {z, zValues}],
1],
#[[1]] > #
[[2]] &];
(*提取符合条件的 y 和 z 值*)
yValuesFiltered = combinations[[All, 1]];
zValuesFiltered = combinations[[All, 2]];
(*启动并行计算内核*)
LaunchKernels[96];
(*并行计算结果并存储在列表中*)
results =
ParallelTable[
Module[{gamma = gammaValue, nz1Vals, nz2Vals, nz3Vals,
nz4Vals, \[CapitalPhi]y1Vals, \[CapitalPhi]z2Vals, \
\[CapitalPhi]z3Vals, \[CapitalPhi]z4Vals, nyzVals},
nz1Vals =
Mean[Table[
nz1[gamma, delta, y, z, lf], {y, yValuesFiltered}, {z,
zValuesFiltered}]];
nz2Vals =
Mean[Table[
nz2[gamma, delta, y, z, lf], {y, yValuesFiltered}, {z,
zValuesFiltered}]];
nz3Vals =
Mean[Table[nz3[gamma, delta, y, lf], {y, yValuesFiltered}]];
nz4Vals =
Mean[Table[nz4[gamma, delta, z, lf], {z, zValuesFiltered}]];
\[CapitalPhi]y1Vals =
Mean[Table[\[CapitalPhi]y1[gamma, delta, y, lf], {y,
yValuesFiltered}]];
\[CapitalPhi]z2Vals =
Mean[Table[\[CapitalPhi]z2[gamma, delta, y, z, lf], {y,
yValuesFiltered}, {z, zValuesFiltered}]];
\[CapitalPhi]z3Vals =
Mean[Table[\[CapitalPhi]z3[gamma, delta, y, lf], {y,
yValuesFiltered}]];
\[CapitalPhi]z4Vals =
Mean[Table[\[CapitalPhi]z4[gamma, delta, z, lf], {z,
zValuesFiltered}]];
nyzVals =
Mean[Table[
nyz[gamma, delta, y, z, lf], {y, yValuesFiltered}, {z,
zValuesFiltered}]];
{gammaValue, nz1Vals, nz2Vals, nz3Vals,
nz4Vals, \[CapitalPhi]y1Vals, \[CapitalPhi]z2Vals, \
\[CapitalPhi]z3Vals, \[CapitalPhi]z4Vals, nyzVals}], {gammaValue,
gammaValues}];
(*关闭并行计算内核*)
CloseKernels[];
(*输出结果*)
results

这是我在pc上截图下来卡主的位置，在工作站上是96核，20分钟还是同一个状态

10个γ值是对应10组结果

Mean[Table[nyz[gamma, delta, y, z, lf], {y, yValuesFiltered}, {z, zValuesFiltered}]
Table的结果会是一个二维列表, Mean作用在高维列表上不会得出所有元素的平均值, 我想应该为
Mean[Flatten@Table[nyz[gamma, delta, y, z, lf], {y, yValuesFiltered}, {z, zValuesFiltered}]
这个代码在我这里可以通过并行运算得出结果, 我只是删除了LaunchKernels[96];和CloseKernels[];并且把zValuesFiltered和yValuesFiltered中元素改得相对少