我在拟合阻尼摆振幅数据时，FindMaximum和InterpolatingFunction等报错。FindMaximum应该是精度有问题，但提高FindMaximum工作精度或者修改系统默认精度后，问题仍然未完全解决；InterpolatingFunction是超出了范围，但我不清楚哪里超出了范围。这两种错出现在某些特定的循环次数。在循环进行20次左右后，振幅数据突然变的很大。
以下是我的代码：
Clear["Global`*"]
g = 9.8;(* 重力加速度 *)
L = 1;(* 摆长 *)
\[CapitalOmega] = Sqrt[g/L];(* 角频率 *)
time = 80;(* 计算时间 *)
\[Eta] = 0.3;(* 阻力系数 *)
{\[Theta]0, \[Omega]0} = {\[Pi]/3, 0};(* 初始角度和初始角速度 *)
s = NDSolve[{\[Theta]''[
t] == -\[CapitalOmega]^2 Sin[\[Theta][t]] - \[Eta] \[Theta]'[
t], \[Theta][0] == \[Theta]0, \[Theta]'[
0] == \[Omega]0}, \[Theta], {t, 0, time}];
\[Theta] = \[Theta] /. s[[1]];
t0 = 0;
\[Delta]t = 2.1; peak = {};
Do[
\[Theta]m =
FindMaximum[\[Theta][t], {t, \[Delta]t + t0},
WorkingPrecision -> 20];
(* Print["\[Theta]m=",\[Theta]m]; *)
AppendTo[peak, {t /. \[Theta]m[[2]], \[Theta]m[[1]]}];
(* Print["peak=",peak]; *)
t1 = t /. \[Theta]m[[2]];
\[Delta]t = t1 - t0;
t0 = t1,
{j, 37}]
g1 = ListPlot[peak, PlotRange -> All];
T = Table[peak[[j, 1]] - peak[[j - 1, 1]], {j, 2, 37}];(* period *)
Partition[T, 4] // TableForm
peak = FindFit[peak, {a E^(-b t), a > 0 && b < 0}, {a, b}, t]
g2 = Plot[a E^(-b t) /. peak, {t, 0, time}, PlotRange -> All];
Show[{g1, g2}, AxesLabel -> {"t/s", "\[Theta]m/rad"}] 
