In[]:=
CompoundExpression[
]
​​deploy
Thu 9 Mar 2023 16:37:36

Laplace transform of sum

Util


Main

https://mathematica.stackexchange.com/questions/282011/efficient-powers-of-dpr1-matrices
Compute
tr(
s
A
)=
∑
i
s
(1-h(i))
​​Steps 1. Write
k
D
ii
as
k
(1-g(i))
≈exp(-kg(i))
​2. Replace sum over $i$ with integral and do a change of variables $y=g(i)$3. Use Heavy-side function to expand integration limits​
In[]:=
d=5;​​h[i_]=
-1
i
;​​hvals=Table[h[i],{i,1,d}];​​A=DiagonalMatrix[Table[1-h[i],{i,1,d}]];​​Print["A=",MatrixForm[A]];​​​​sum[s_]:=Tr[MatrixPower[A,s]]/d;​​intOriginal=Inactive[Integrate][Exp[-sh[i]]/d,{i,1,d}];​​​​(*changeintointegraloverhdensity*)​​intLaplace=IntegrateChangeVariables[intOriginal,hi,hi==h[i]];​​(*expandlimitsto0,∞*)​​intLaplace=expandIntegrationLimits[intLaplace];​​​​(*minustoworkaroundbughttps://mathematica.stackexchange.com/questions/270358/integratechangevariables-producing-incorrect-result*)​​intLaplace=-LaplaceTransform[First@intLaplace/Exp[-his],hi,s];​​​​Print["Exact value diag: ",nf@sum[d]]​​Print["Integral: ",nf@Activate[intOriginal/.{Integrate->NIntegrate,s->d}]];​​Print["Laplace integral: ",nf[intLaplace/.s->d]];​​Print["Laplace formula: ",intLaplace];​​​​matrixPowers=NestList[A.#&,A,d];​​actualPlot=DiscretePlot[sum[s],{s,1,d},AxesOrigin->{0,0},PlotLegends->{"Tr
s
A
"}];​​laplacePlot=Plot[intLaplace,{s,1,d}];​​Show[actualPlot,laplacePlot]
A=
0
0
0
0
0
0
1
2
0
0
0
0
0
2
3
0
0
0
0
0
3
4
0
0
0
0
0
4
5
Exact value diag: 0.146
Integral: 0.148
Laplace integral: 0.148
Laplace formula:
1
5
-
-s

+5
-s/5

-sGamma0,
s
5
+sGamma[0,s]
Out[]=
Tr
s
A

DPR1 version (not done)

(*normalizehtoensureDPR1matrixisconvergent*)​​hvals=hvals/Total[hvals];​​dpr[s_]:=Tr[MatrixPower[DiagonalMatrix[1-hvals]+{hvals}.{hvals},s]]/d;​​s0=5;​​Print["Exact value dpr1: ",dpr[s0]//N];​​Print["Approx value dpr1: ",dprApprox[s0]];

(original version of Laplace sum, manual change of variables)
