In[]:=
CompoundExpression[
]
​​deploy
Tue 9 May 2023 00:16:08
Using exp(-At) to approximate (I-A) math.SE post
In[]:=
d=10000;​​maxSteps=1000;​​p=1.;​​h=Table[
-p
i
,{i,1,d}];​​h=h/Max[h];​​exact[t_]:=Total[
t
(1-h)
];​​approx[t_]:=Total[Exp[-th]];​​error[truth_,est_]:=(est-truth);​​relError[truth_,est_]:=
(est-truth)
truth
;​​Print["Tr(A)=",Total[h]]​​Print["Tr(​
2
A
​)=",Total[h*h]]​​ListPlot[h,PlotLabel->"A eigenvalues",ScalingFunctions->{"Log","Log"},Filling->Axis];​​LogPlot[error[exact[s],approx[s]],{s,1,maxSteps},PlotLabel->"Error of approximating (I-A
t
)\),
with
-tA
e
",AxesLabel->{"t","error"}]​​
Tr(A)=9.78761
Tr(​
2
A
​)=1.64483
Out[]=
In[]:=
xvals=Table[t,{t,1,maxSteps}];​​yvals=Table[error[exact[t],approx[t]],{t,xvals}];​​observedPlot=ListPlot[{xvals,yvals},PlotLegends->{"
f
A
​(t)"}];​​formula=
1
2
Exp[-x/d];​​predictedPlot=Plot[formula,{x,0,d},Filling->Axis,PlotLegends->{formula}];​​observedPlot=ListPlot[{xvals[[;;;;100]],yvals[[;;;;100]]},PlotStyle->{Red}];​​Show[observedPlot,predictedPlot,PlotLabel->"Error of approximating (I-A
t
)\),
with
-tA
e
",AxesLabel->{"t","error"}]​​
Out[]=
-
x
10000

2
pdf=
-1/p
(y)
py
;​​int=Inactive[Integrate][pdfExp[-ys],{y,
-p
(d+1)
,1}];​​formulaError=Assuming[{p>=1,d>1,s>1},Activate[int]]
Out[]=
-
1
p
s
Gamma-
1
p
,s-Gamma-
1
p
,
-p
(1+d)
s
p
if
p
(1+d)
≠1
In[]:=
errorPlot=Block[{p=1,d=10000},Plot[formulaError,{s,1,maxSteps}]]
Out[]=

Power-Law fits

In[]:=
powerLawFit[xvals_,yvals_]:=Module[{logData,linearFit,a,b},logData=Transpose[{Log[xvals],Log[yvals]}];​​linearFit=LinearModelFit[logData,x,x];​​{a,b}=linearFit["BestFitParameters"];​​fit=Exp[a+b*Log[x]];​​fitModel=Function@@{{x},fit};​​fitModel​​];​​exponentialFit[xvals_,yvals_]:=Module[{logData,linearFit,a,b},​​logData=Transpose[{xvals,Log[yvals]}];​​linearFit=LinearModelFit[logData,x,x];​​{a,b}=linearFit["BestFitParameters"];​​fit=Exp[a+b*x];​​fitModel=Function@@{{x},fit};​​fitModel​​];​​​​f1=powerLawFitxvals;;
d
,yvals;;
d
;​​f2=exponentialFit[xvals[[maxSteps/2;;maxSteps-maxSteps/4]],yvals[[maxSteps/2;;maxSteps-maxSteps/4]]];​​fittedPlot=Plot[{None,f1[t],f2[t]},{t,0,maxSteps},PlotLegends->{None,f1[t],f2[t]}];​​Show[observedPlot,fittedPlot,PlotRange->{{0,maxSteps},{0,1}}]
Out[]=
0.548933
0.0243381
t
-0.000100435t-0.692587

In[]:=
f2[x]
Out[]=
-0.692587-0.000100435x
