Get["~/doc/research/first_order_singularities/mma/schofield.wl"]

SetAttributes[untilSuccess,HoldAll];
untilSuccess[e_]:=Quiet[Check[e,untilSuccess[e]]]

randomEval[eq_]:=eq/.Around[x_,y_]:>RandomVariate[NormalDistribution[x,y]]

Matching all

Matching all (Subscript[θ, c]) fixed

Match ends

nAs=2;

s0=#1->Re[#2]&@@@untilSuccess[newSol[randomEval[eqAround[2,h[0]][1,-1,2]],{θc->1+RandomReal[],θ0->0.2},{AH,A[1],A[2]},.1]]

Timing[dGdξLowList[nAs,h[0]][6]/(Range[0,6]!)//.rules[h[0]]/.s0]

ListLogPlot[{
Abs[Gls],
Chop[Abs[dGdξLowList[nAs,h[0]][6]]//.rules[h[0]]/.s0]
}]

ListLogPlot[{
Abs[Ghs],
Chop[Abs@dGdξList[nAs,h[0]][6,0]//.rules[h[0]]/.s0]
}]

res[2,h[0],∞][2]

m0=annealFit[Most@res[2,h[0],d][2]//.rules[h[0]],d,s0]

m01=FindRoot[Evaluate[D[Total[Most@((res[2,h[0],(10^(-15))][2])^2)]//.rules[h[0]],{m0[[2,All,1]]}]],Evaluate[m0[[2]]/.Rule->List],WorkingPrecision->20,MaxIterations->5000]

fitCov[res_,s_]:=With[{J=D[res,{s[[All,1]]}]},Symmetrize[Inverse[SetPrecision[Transpose[J].J/.s,60]]]]

cov0=fitCov[Most@res[2,h[0],(10^(-15))][2]//.rules[h[0]],m01]

Eigenvalues[cov0]

m1=annealFit[Most@res[2,h[1],δ][3]//.rules[h[1]],δ,Join[m0[[2]],{gC[1]->0}]]

m11=FindRoot[Evaluate[D[Total[Most@((res[2,h[1],(10^(-15))][3])^2)]//.rules[h[1]],{m1[[2,All,1]]}]],Evaluate[m1[[2]]/.Rule->List],WorkingPrecision->20,MaxIterations->5000]

cov1=fitCov[Most@res[2,h[1],(10^(-15))][3]//.rules[h[1]],m11];

Eigenvalues[cov1]

Total[Most@((res[2,h[1],(10^(-15))][3])^2)]//.rules[h[1]]/.m11

m2=annealFit[Most@res[3,h[2],δ][4]//.rules[h[2]],δ,Join[m1[[2]],{gC[2]->0,A[3]->0}]]

m21=FindRoot[Evaluate[D[Total[Most@((res[2,h[3],(10^(-13))][4])^2)]//.rules[h[3]],{m2[[2,All,1]]}]],Evaluate[m2[[2]]/.Rule->List],WorkingPrecision->15,MaxIterations->5000]

cov2=fitCov[Most@res[2,h[3],(10^(-13))][4]//.rules[h[3]],m21];

Eigenvalues[cov2]

Total[Most@((res[2,h[3],(10^(-15))][4])^2)]//.rules[h[3]]/.m21

m3=annealFit[Join[Most@res[3,h[3],δ][5]]//.rules[h[3]],δ,Join[m2[[2]],{gC[3]->0}]]

m31=FindRoot[Evaluate[D[Total[Most@((res[2,h[4],(10^(-15))][5])^2)]//.rules[h[4]],{m3[[2,All,1]]}]],Evaluate[m3[[2]]/.Rule->List],WorkingPrecision->20,MaxIterations->5000]

cov3=fitCov[Most@res[2,h[4],(10^(-15))][5]//.rules[h[4]],m31];

Eigenvalues[cov3]

Total[Most@((res[2,h[3],(10^(-15))][5])^2)]//.rules[h[3]]/.m31

m4=annealFit[Join[Most@res[4,h[4],δ][6]//.rules[h[4]],Table[gC[i](((i!))^2),{i,1,4}]],δ,Join[m3[[2]],{gC[4]->0,A[4]->0}]]

m41=FindRoot[Evaluate[D[Total[Most@((res[2,h[6],(10^(-15))][])^2)]//.rules[h[6]],{m4[[2,All,1]]}]],Evaluate[m4[[2]]/.Rule->List],WorkingPrecision->20,MaxIterations->5000]

cov4=fitCov[Most@res[2,h[5],(10^(-13))][6]//.rules[h[5]],m41];

Total[Most@((res[2,h[5],(10^(-15))][6])^2)]//.rules[h[5]]/.m41

m5=annealFit[Most@res[4,h[5],δ][7]//.rules[h[5]],δ,Join[m4[[2]],{}]]

res[2,h[6],(10^(-15))][7]//Length

m5[[2]]//Length

m51=FindRoot[Evaluate[D[Total[((res[2,h[7],(10^(-15))][7])^2)]//.rules[h[7]],{m5[[2,All,1]]}]],Evaluate[m5[[2]]/.Rule->List],WorkingPrecision->20,MaxIterations->5000]

Total[Most@((res[2,h[7],(10^(-15))][7])^2)]//.rules[h[7]]/.m51

m6=annealFit[Join[Most@res[5,h[6],δ][8]//.rules[h[6]],Table[gC[i](i!),{i,1,6}]],δ,Join[m5[[2]],{gC[6]->0,A[5]->0}]]

m6={2272.734589158684`,{θc->1.3053061442708302`,θ0->0.19663560107595623`,AH->-2.5368261599243653`,A[1]->1.7108160252330362`,A[2]->0.5926167832063681`,gC[1]->-0.2812358144995081`,gC[2]->-0.04198949750026452`,A[3]->-0.17345114501721076`,gC[3]->0.0065596668614836505`,gC[4]->0.004317004580958006`,gC[5]->-0.0014157584142083521`,A[4]->0.006658268409302458`,gC[6]->0.00012139778301924035`,A[5]->0.0009751352780293256`}}

m61=FindRoot[Evaluate[D[Total[((res[2,h[9],(10^(-15))][8])^2)]//.rules[h[9]],{m6[[2,All,1]]}]],Evaluate[m6[[2]]/.Rule->List],WorkingPrecision->20,MaxIterations->5000]

Total[((res[2,h[9],(10^(-15))][8])^2)]//.rules[h[9]]/.m61

m7=annealFit[Most@res[5,h[7],δ][9]//.rules[h[7]],δ,Join[m6[[2]],{gC[7]->0}]]

m8=annealFit[res[6,h[8],δ][10]//.rules[h[8]],δ,Join[m7[[2]],{gC[8]->0,A[6]->0}]]

ntest=3;n1=2;
ourΦ1=Map[With[{df=D[#,{m01[[All,1]]}]},Around[#,df.cov0.df]]&,dΦdηList[n1,h[0]][ntest,1]]//.rules[h[0]]/.m01;
ourΦ2=Map[With[{df=D[#,{m11[[All,1]]}]},Around[#,df.cov1.df]]&,dΦdηList[n1,h[1]][ntest,1]]//.rules[h[1]]/.m11;
ourΦ3=Map[With[{df=D[#,{m21[[All,1]]}]},Around[#,df.cov2.df]]&,dΦdηList[n1,h[3]][ntest,1]]//.rules[h[3]]/.m21;
ourΦ4=Map[With[{df=D[#,{m31[[All,1]]}]},Around[#,df.cov3.df]]&,dΦdηList[n1,h[3]][ntest,1]]//.rules[h[3]]/.m3[[2]];
ourΦ5=dΦdηList[n1,h[5]][ntest,1]//.rules[h[5]]/.m41;
ourΦ6=dFdηList[n1,h[5]][ntest][1]/Range[0,ntest]!//.rules[h[5]]/.m5[[2]];
ourΦ7=dFdηList[n1,h[7]][ntest][1]/Range[0,ntest]!//.rules[h[7]]/.m6[[2]];
test=((Φs[[;;ntest+1]]-#)/Φs[[;;ntest+1]])&/@{ourΦ1,ourΦ2,ourΦ3,ourΦ4,ourΦ5,ourΦ6,ourΦ7};

ntest=2;n1=2;
ourΦ1=dΦdηList[n1,h[0]][ntest,1]//.rules[h[0]]/.m01;
ourΦ2=dΦdηList[n1,h[1]][ntest,1]//.rules[h[1]]/.m11;
ourΦ3=dΦdηList[3,h[2]][ntest,1]//.rules[h[2]]/.m2[[2]];
ourΦ4=dΦdηList[3,h[2]][ntest,1]//.rules[h[2]]/.m3[[2]];
ourΦ5=dΦdηList[4,h[4]][ntest,1]//.rules[h[4]]/.m4[[2]];
ourΦ6=dΦdηList[4,h[5]][ntest,1]//.rules[h[5]]/.m5[[2]];
ourΦ7=dΦdηList[5,h[6]][ntest,1]//.rules[h[6]]/.m6[[2]];
ourΦ8=dΦdηList[5,h[7]][ntest,1]//.rules[h[7]]/.m7[[2]];
test=((Φs[[;;ntest+1]]-#)/Φs[[;;ntest+1]])&/@{ourΦ1,ourΦ2,ourΦ3,ourΦ4,ourΦ5,ourΦ6,ourΦ7,ourΦ8};

((Φs[[;;2+1]]-#)/Φs[[;;2+1]])&@(dΦdηList[5,h[6]][2,1]//.rules[h[6]]/.m6[[2]])

((Φs[[;;2+1]]-#)/Φs[[;;2+1]])&@(dΦdηList[5,h[7]][2,1]//.rules[h[7]]/.m7[[2]])

Show[ListLogLogPlot[Abs@Re@Transpose[test],Joined->True],LogLogPlot[0.0005(x^(-2)),{x,0.5,10},PlotStyle->Black]]