求助:在程序中途出现“精度不够”的情况,附代码文本
mathematica吧
全部回复
仅看楼主
level 1
ClearAll["`*"]
eq0[\[Mu]_] :=
1 - \[Omega] *t^(1/2)*
Sum[Sum[(-1)^(l + 1)*
l^(-1/2)*(Exp[-l*\[Omega]*\[Beta]*(1/2 -
g*\[Sigma] - \[Mu]/\[Omega])])/(1 -
Exp[-l*\[Omega]*\[Beta]]), {l, 1, 400}], {\[Sigma], {-(1/2),
1/2}}];
g = 0.1;
\[Omega] = 0.01;
\[Beta] = 1/t;
x = \[Omega]/t;
K = {};
SetDirectory["E:\\yangjiayingsan\\"];
outfile = OpenAppend["cv-b g0.1 w0.01", PageWidth -> Infinity];
{\[Mu]i} = {-1};
Do[
{t = 10 - 0.01*nn;
{\[Mu]i} = {\[Mu]} = {\[Mu]} /.
FindRoot[
{
eq0[\[Mu]]
},
{{\[Mu], \[Mu]i}}
];
cv =
x*t^(3/2)*
Sum[((3/4)*
Sum[(-1)^(l + 1)*l^(-3/2)*
Exp[-(l/t)*(1/2*\[Omega] - g*\[Sigma]*\[Omega] - \[Mu])]/(
1 - Exp[-l*x]), {l, 1, 400}] +
x*Sum[(-1)^(l + 1)*l^(-1/2)*
Exp[-(l/
t)*(1/2*\[Omega] -
g*\[Sigma]*\[Omega] - \[Mu] + \[Omega])]/(1 -
Exp[-l*x])^2, {l, 1, 400}] +
x^2*Sum[(-1)^(l + 1)*l^(1/2)*
Exp[-(l/
t)*(1/2*\[Omega] -
g*\[Sigma]*\[Omega] - \[Mu] + \[Omega])]/(1 -
Exp[-l*x])^2, {l, 1, 400}] +
2*x^2*Sum[(-1)^(l + 1)*l^(1/2)*
Exp[-(l/
t)*(1/2*\[Omega] - g*\[Sigma]*\[Omega] - \[Mu] +
2*\[Omega])]/(1 - Exp[-l*x])^3, {l, 1, 400}] -
x^2*(Sum[(-1)^(l + 1)*l^(1/2)*
Exp[-(l/
t)*(1/2*\[Omega] -
g*\[Sigma]*\[Omega] - \[Mu] + \[Omega])]/(1 -
Exp[-l*x])^2, {l, 1, 400}])^2/
Sum[(-1)^(l + 1)*l^(1/2)*
Exp[-(l/t)*(1/2*\[Omega] - g*\[Sigma]*\[Omega] - \[Mu])]/(
1 - Exp[-l*x]), {l, 1,
400}] - (Sum[(-1)^(l + 1)*l^(-1/2)*
Exp[-(l/
t)*(1/2*\[Omega] - g*\[Sigma]*\[Omega] - \[Mu])]/(
1 - Exp[-l*x]), {l, 1, 400}])^2/(4*
Sum[(-1)^(l + 1)*l^(1/2)*
Exp[-(l/
t)*(1/2*\[Omega] - g*\[Sigma]*\[Omega] - \[Mu])]/(
1 - Exp[-l*x]), {l, 1, 400}]) -
x*(Sum[(-1)^(l + 1)*l^(-1/2)*
Exp[-(l/t)*(1/2*\[Omega] - g*\[Sigma]*\[Omega] - \[Mu])]/(
1 - Exp[-l*x]), {l, 1,
400}])*(Sum[(-1)^(l + 1)*l^(1/2)*
Exp[-(l/
t)*(1/2*\[Omega] -
g*\[Sigma]*\[Omega] - \[Mu] + \[Omega])]/(1 -
Exp[-l*x])^2, {l, 1, 400}])/
Sum[(-1)^(l + 1)*l^(1/2)*
Exp[-(l/t)*(1/2*\[Omega] - g*\[Sigma]*\[Omega] - \[Mu])]/(
1 - Exp[-l*x]), {l, 1, 400}]), {\[Sigma], {-(1/2), 1/2}}];
K = AppendTo[K, {\[Beta], cv}];
Write[outfile, {\[Omega], g, \[Beta], t, \[Mu], \[Beta]*\[Mu], cv}];
Print[{\[Omega], g, \[Beta], t, \[Mu], \[Beta]*\[Mu], cv}]}, {nn,
0, 999}];
ListPlot[K, Joined -> True, AxesLabel -> {"\[Beta]", "cv"},
PlotRange -> All]
Close[outfile];
2020年12月22日 02点12分 1
level 1
其中t是关于u的函数,g,w,β为变量,β=1/t
附:出错位置,在温度t等于0.75时,理论上最后一项cv项应该保持下降趋势。
2020年12月22日 02点12分 2
level 1
不会
2020年12月23日 07点12分 3
1