求助 一个令人费解的错误，它是如何产生的？

halfIntegerQ[j_] := IntegerQ[2 j] && j >= 1/2;
jplus[jmax_?halfIntegerQ] :=
SparseArray[
Band[{1, 2}] ->
Table[N@Sqrt[jmax (jmax + 1) - m (m + 1)], {m,
jmax - 1, -jmax, -1}], {2 jmax + 1, 2 jmax + 1}];
jminus[jmax_?halfIntegerQ] := Transpose@jplus[jmax];
jx[jmax_?halfIntegerQ] := (jplus[jmax] + jminus[jmax])/2;
jy[jmax_?halfIntegerQ] := (jplus[jmax] - jminus[jmax])/(2 I);
jz[jmax_?halfIntegerQ] :=
SparseArray[
Band[{1, 1}] -> N@Range[jmax, -jmax, -1], {2 jmax + 1, 2 jmax + 1}];
jid[jmax_?halfIntegerQ] := IdentityMatrix[2 jmax + 1, SparseArray];
hamiltonian[h_, spinNumber_?IntegerQ] := -h*jx[spinNumber/2] -
MatrixPower[jx[spinNumber/2], 2]/spinNumber;
initialState[spinNumber_?IntegerQ] :=
Block[{\[CapitalPi], e, v, psi01, psi02}, \[CapitalPi] =
MatrixExp[
I*Pi*(spinNumber/2*jid[spinNumber/2] + jz[spinNumber/2])];
{e, v} = Eigensystem[jx[spinNumber/2], 2];
{psi01, psi02} = {Normalize[v[[1]] + v[[2]]],
Normalize[v[[1]] - v[[2]]]};
If[Total@Abs[\[CapitalPi] . psi01 - psi01] < 10.^(-5), psi01,
psi02]];
ketToRho[psi_] := Outer[Times, psi, Conjugate@psi];
problem[rho_?MatrixQ] :=
Block[{jmax, opj, \[CapitalDelta], points, grid, data, ifun, w0,
w1}, jmax = (First@Dimensions@rho - 1)/2;
opj[\[Phi]_, \[Theta]_] :=
Sin[\[Theta]]*Cos[\[Phi]]*jx[jmax] +
Sin[\[Theta]]*Sin[\[Phi]]*jy[jmax] + Cos[\[Theta]]*jz[jmax];
\[CapitalDelta][index_] := \[CapitalDelta][index] =
Block[{m = index - 1 - jmax}, (1./(2.*jmax + 1.))*
Sum[(2 l + 1)*ClebschGordan[{jmax, m}, {l, 0}, {jmax, m}], {l,
0, 2*jmax, 1}]];
points[\[Phi]_, \[Theta]_] :=
Block[{e, v, kernel}, {e, v} =
Transpose@
SortBy[Transpose@Eigensystem[opj[\[Phi], \[Theta]]], First];
kernel =
Total@MapIndexed[\[CapitalDelta][First[#2]]*
Outer[Times, #1, Conjugate@#1] &, v];
Re@Tr[rho . kernel]];
grid =
Flatten[CoordinateBoundsArray[{{0., 2.*Pi}, {0, 1.*Pi}}, {Into[
401], Into[201]}], 1];
data = MapApply[{##, Quiet@points[##]} &, grid];
ifun = Interpolation[data];
w0 = (2.*jmax + 1.)/(4*Pi)*
Quiet@NIntegrate[
ifun[\[Phi], \[Theta]]*Sin[\[Theta]], {\[Phi], 0,
2 Pi}, {\[Theta], 0, Pi}];
w1 = 0.5*(2.*jmax + 1.)/(4*Pi)*(Quiet@
NIntegrate[
Abs[ifun[\[Phi], \[Theta]]]*Sin[\[Theta]], {\[Phi], 0,
2 Pi}, {\[Theta], 0, Pi}]) - 0.5;
{w0, w1}];
test[n_?IntegerQ] := problem@ketToRho@initialState[n];

……你的Parallelize在哪呢？

halfIntegerQ[j_] := IntegerQ[2 j] && j >= 1/2;
jplus[jmax_?halfIntegerQ] :=
SparseArray[
Band[{1, 2}] ->
Table[N@Sqrt[jmax (jmax + 1) - m (m + 1)], {m,
jmax - 1, -jmax, -1}], {2 jmax + 1, 2 jmax + 1}];
jminus[jmax_?halfIntegerQ] := Transpose@jplus[jmax];
jx[jmax_?halfIntegerQ] := (jplus[jmax] + jminus[jmax])/2;
jy[jmax_?halfIntegerQ] := (jplus[jmax] - jminus[jmax])/(2 I);
jz[jmax_?halfIntegerQ] :=
SparseArray[
Band[{1, 1}] -> N@Range[jmax, -jmax, -1], {2 jmax + 1,
2 jmax + 1}];
jid[jmax_?halfIntegerQ] := IdentityMatrix[2 jmax + 1, SparseArray];
hamiltonian[h_, spinNumber_?IntegerQ] := -h*jx[spinNumber/2] -
MatrixPower[jx[spinNumber/2], 2]/spinNumber;
initialState[spinNumber_?IntegerQ] :=
Block[{\[CapitalPi], e, v, psi01, psi02}, \[CapitalPi] =
MatrixExp[
I*Pi*(spinNumber/2*jid[spinNumber/2] + jz[spinNumber/2])];
{e, v} = Eigensystem[jx[spinNumber/2], 2];
{psi01, psi02} = {Normalize[v[[1]] + v[[2]]],
Normalize[v[[1]] - v[[2]]]};
If[Total@Abs[\[CapitalPi] . psi01 - psi01] < 10.^(-5), psi01,
psi02]];
ketToRho[psi_] := Outer[Times, psi, Conjugate@psi];
problem[rho_?MatrixQ] :=
Block[{jmax, opj, \[CapitalDelta], points, grid, data, ifun, w0,
w1}, jmax = (First@Dimensions@rho - 1)/2;
opj[\[Phi]_, \[Theta]_] :=
Sin[\[Theta]]*Cos[\[Phi]]*jx[jmax] +
Sin[\[Theta]]*Sin[\[Phi]]*jy[jmax] + Cos[\[Theta]]*jz[jmax];
\[CapitalDelta][index_] := \[CapitalDelta][index] =
Block[{m = index - 1 - jmax}, (1./(2.*jmax + 1.))*
Sum[(2 l + 1)*ClebschGordan[{jmax, m}, {l, 0}, {jmax, m}], {l,
0, 2*jmax, 1}]];
points[\[Phi]_, \[Theta]_] :=
Block[{e, v, kernel}, {e, v} =
Transpose@
SortBy[Transpose@Eigensystem[opj[\[Phi], \[Theta]]], First];
kernel =
Total@MapIndexed[\[CapitalDelta][First[#2]]*
Outer[Times, #1, Conjugate@#1] &, v];
Re@Tr[rho . kernel]];
grid =
Flatten[CoordinateBoundsArray[{{0., 2.*Pi}, {0, 1.*Pi}}, {Into[
401], Into[201]}], 1];
data = Parallelize@MapApply[{##, Quiet@points[##]} &, grid];
ifun = Interpolation[data];
w0 = (2.*jmax + 1.)/(4*Pi)*
Quiet@NIntegrate[
ifun[\[Phi], \[Theta]]*Sin[\[Theta]], {\[Phi], 0,
2 Pi}, {\[Theta], 0, Pi}];
w1 = 0.5*(2.*jmax + 1.)/(4*Pi)*(Quiet@
NIntegrate[
Abs[ifun[\[Phi], \[Theta]]]*Sin[\[Theta]], {\[Phi], 0,
2 Pi}, {\[Theta], 0, Pi}]) - 0.5;
{w0, w1}];
test[n_?IntegerQ] := problem@ketToRho@initialState[n];
我的问题: 第一个积分w0正确的结果是1, 以这个结果为准绳的话, 上面的代码只有在首次运算时才能返回正确的结果, 我想不明白错误是如何产生的

“第一个积分w1在数学上的结果是1”，你是想说“第一个积分w0在数学上的结果是1”吧？
那你大可以把w1的定义给去掉啊，这不就少了一行了吗？

我努力缩减了代码
problem[n_?IntegerQ] := Block[{rho, \[CapitalDelta], points, data},
rho = (1/n)*DiagonalMatrix@Range[1, n];
\[CapitalDelta][index_] := \[CapitalDelta][index] =
First@Dimensions@rho + index;
points[\[Phi]_, \[Theta]_] := \[CapitalDelta][1]*rho[[1, 1]];
data = Parallelize@MapApply[points[##] &, {{1, 1}}];
data];
{problem[1], problem[2], problem[3]}
这个代码输入整数n, 预期结果应该是 1 + (1/n) 但我们发现结果不是这样的, 只有第一项的结果符合预期
我们删除代码中的Parallelize, 得到problem2, 会给出正确结果
problem2[n_?IntegerQ] := Block[{rho, \[CapitalDelta], points, data},
rho = (1/n)*DiagonalMatrix@Range[1, n];
\[CapitalDelta][index_] := \[CapitalDelta][index] =
First@Dimensions@rho + index;
points[\[Phi]_, \[Theta]_] := \[CapitalDelta][1]*rho[[1, 1]];
data = MapApply[points[##] &, {{1, 1}}];
data];
{problem2[1], problem2[2], problem2[3]}
另外我们将代码中 Δ[index_] := Δ[index] = … 改为Δ[index_] := , 得到problem3, 也会给出正确结果
problem3[n_?IntegerQ] := Block[{rho, \[CapitalDelta], points, data},
rho = (1/n)*DiagonalMatrix@Range[1, n];
\[CapitalDelta][index_] := First@Dimensions@rho + index;
points[\[Phi]_, \[Theta]_] := \[CapitalDelta][1]*rho[[1, 1]];
data = MapApply[points[##] &, {{1, 1}}];
data];
{problem3[1], problem3[2], problem3[3]}
实在难以想清楚为什么Parallelize和Δ[index_] := Δ[index] =为什么有冲突, 还烦请点拨一下

\[CapitalDelta]需要share，代码里加一行
SetSharedFunction[\[CapitalDelta]];
就没事了。参SE帖子
《Assignments in ParallelDo when using indexed variables》（编号80325）
顺便，联用并行和记忆式的函数定义似乎可能会导致效率问题，参SE帖子
《Parallelize evaluation of function with memoization》（编号1259）