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[nAs,h[0]][1,-1,2]],{θc->1+RandomReal[],θ0->0.2},{AH,A[1],A[2]},.1]]

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

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

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

m0=levenburgMarquardt[res[2,h[0],(10^(-6))][2]//.rules[h[0]],s0]

m01=levenburgMarquardt[res[2,h[0],(10^(-15))][2]//.rules[h[0]],m0[[2]]]

m1=levenburgMarquardt[res[2,h[1],(10^(-7))][3]//.rules[h[1]],Append[m0[[2]],gC[1]->0]]

m11=levenburgMarquardt[res[2,h[1],(10^(-15))][3]//.rules[h[1]],m1[[2]]]

m2=levenburgMarquardt[res[2,h[3],(10^(-6))][4]//.rules[h[3]],Join[m1[[2]],{gC[2]->0,gC[3]->0}]]

m21=levenburgMarquardt[res[2,h[3],(10^(-15))][4]//.rules[h[3]],Join[m2[[2]]]]

m3=levenburgMarquardt[res[2,h[4],(10^(-5))][5]//.rules[h[4]],Join[m2[[2]],{gC[4]->0}]]

m31=levenburgMarquardt[res[2,h[4],(10^(-7))][5]//.rules[h[4]],Join[m3[[2]]]]

m32=levenburgMarquardt[res[2,h[4],(10^(-9))][5]//.rules[h[4]],Join[m31[[2]]]]

m33=levenburgMarquardt[res[2,h[4],(10^(-15))][5]//.rules[h[4]],Join[m32[[2]]]]

m4=levenburgMarquardt[res[2,h[6],(10^(-5))][6]//.rules[h[6]],Join[m31[[2]],{gC[5]->0,gC[6]->0}]]

m41=levenburgMarquardt[res[2,h[6],(10^(-9))][6]//.rules[h[6]],Join[m4[[2]],{}]]

m42=levenburgMarquardt[res[2,h[6],(10^(-15))][6]//.rules[h[6]],Join[m41[[2]],{}]]

m5=levenburgMarquardt[res[2,h[7],(10^(-5))][7]//.rules[h[7]],Join[m4[[2]],{gC[7]->0}]]

m51=levenburgMarquardt[res[2,h[7],(10^(-9))][7]//.rules[h[7]],Join[m5[[2]],{}]]

m52=levenburgMarquardt[res[2,h[7],(10^(-15))][7]//.rules[h[7]],Join[m51[[2]],{}],1,2,(10^(-16))]

m5={1.9017425705133395`,{θc->1.4205660602855785`,θ0->0.19819854132109802`,AH->-2.4450078675776648`,A[1]->1.5723900032479108`,A[2]->0.2774470302292178`,gC[1]->-0.17288287974478225`,gC[2]->0.13146623874166138`,gC[3]->-0.10479974533776427`,gC[4]->0.05727932658366424`,gC[5]->-0.019428278795753182`,gC[6]->0.003684613297584252`,gC[7]->-0.0002991412844663506`}};

m6=levenburgMarquardt[res[2,h[9],(10^(-5))][8]//.rules[h[9]],Join[m5[[2]],{gC[8]->0,gC[9]->0}]]

m61=levenburgMarquardt[res[2,h[9],(10^(-9))][8]//.rules[h[9]],Join[m6[[2]],{}]]

ntest=3;n1=2;
ourΦ1=dFdηList[n1,h[0]][ntest][1]/Range[0,ntest]!//.rules[h[0]]/.m01[[2]];
ourΦ2=dFdηList[n1,h[1]][ntest][1]/Range[0,ntest]!//.rules[h[1]]/.m11[[2]];
ourΦ3=dFdηList[n1,h[3]][ntest][1]/Range[0,ntest]!//.rules[h[3]]/.m21[[2]];
ourΦ4=dFdηList[n1,h[4]][ntest][1]/Range[0,ntest]!//.rules[h[4]]/.m33[[2]];
ourΦ5=dFdηList[n1,h[6]][ntest][1]/Range[0,ntest]!//.rules[h[6]]/.m42[[2]];
ourΦ6=dFdηList[n1,h[7]][ntest][1]/Range[0,ntest]!//.rules[h[7]]/.m52[[2]];
test=(Φs[[;;ntest+1]]-#)/Φs[[;;ntest+1]]&/@{ourΦ1,ourΦ2,ourΦ3,ourΦ4,ourΦ5,ourΦ6};

ListLogPlot[Abs@Re@Transpose[test],Joined->True]

min0=FindMinimum[Evaluate[chiSquared[2,h[0],(10^(-6))][2]//.rules[h[0]]],s0/.Rule->List,Method->"LevenbergMarquardt"]

min1=FindMinimum[Evaluate[chiSquared[2,h[0],(10^(-15))][2]//.rules[h[0]]],min0[[2]]/.Rule->List]

sol[2]=min1[[2]];hf[2]=h[0];

min3=FindMinimum[Evaluate[chiSquared[2,h[0],(10^(-6))][3]//.rules[h[0]]],Join[{#1,Re[#2]}&@@@min0[[2]],{}],MaxIterations->1000]

min4=FindMinimum[Evaluate[chiSquared[2,h[1],(10^(-6))][3]//.rules[h[1]]],Join[{#1,Re[#2]}&@@@min3[[2]],{{gC[1],0}}],MaxIterations->(10^5)]

min5=FindMinimum[Evaluate[chiSquared[2,h[1],(10^(-15))][3]//.rules[h[1]]],Join[{#1,Re[#2]}&@@@min4[[2]],{}],MaxIterations->(10^5)]

sol[3]=min5[[2]];hf[3]=h[1];

min6=FindMinimum[chiSquared[2,h[1],(10^(-6))][4]//.rules[h[1]],Join[{#1,Re[#2]}&@@@min4[[2]],{}],MaxIterations->(10^5)]

min7=FindMinimum[chiSquared[2,h[2],(10^(-5))][4]//.rules[h[2]],Join[{#1,Re[#2]}&@@@min6[[2]],{{gC[2],0}}],MaxIterations->(10^5)]

min8=FindMinimum[chiSquared[2,h[3],(10^(-5))][4]//.rules[h[3]],Join[{#1,Re[#2]}&@@@min7[[2]],{{gC[3],0}}],MaxIterations->(10^5)]

min9=FindMinimum[chiSquared[2,h[3],(10^(-15))][4]//.rules[h[3]],Join[{#1,Re[#2]}&@@@min8[[2]],{}],MaxIterations->(10^5)]

sol[4]=min9[[2]];hf[4]=h[3];

min10=FindMinimum[chiSquared[2,h[3],(10^(-5))][5]//.rules[h[3]],Join[{#1,Re[#2]}&@@@min8[[2]],{}],MaxIterations->(10^5)]

min11=FindMinimum[chiSquared[2,h[4],(10^(-5))][5]//.rules[h[4]],Join[{#1,Re[#2]}&@@@min10[[2]],{{gC[4],0}}],MaxIterations->(10^5)]

min12=FindMinimum[chiSquared[2,h[5],(10^(-5))][5]//.rules[h[5]],Join[{#1,Re[#2]}&@@@min11[[2]],{{gC[5],0}}],MaxIterations->(10^5)]

min13=FindMinimum[chiSquared[2,h[5],(10^(-15))][5]//.rules[h[5]],Join[{#1,Re[#2]}&@@@min12[[2]],{}],MaxIterations->(10^5)]

sol[5]=min13[[2]];hf[5]=h[5];

min14=FindMinimum[chiSquared[2,h[5],(10^(-5))][6]//.rules[h[5]],Join[{#1,Re[#2]}&@@@min12[[2]],{}],MaxIterations->(10^5)]

min16=FindMinimum[chiSquared[2,h[6],(10^(-15))][6]//.rules[h[6]],Join[{#1,Re[#2]}&@@@min15[[2]],{}],MaxIterations->(10^5)]

min15=FindMinimum[chiSquared[2,h[6],(10^(-5))][6]//.rules[h[6]],Join[{#1,Re[#2]}&@@@min14[[2]],{{gC[6],0}}],MaxIterations->(10^5)]

min16=FindMinimum[chiSquared[2,h[7],(10^(-5))][6]//.rules[h[7]],Join[{#1,Re[#2]}&@@@min42[[2]],{{gC[7],0}}],MaxIterations->(10^5)]

sol[6]=min15[[2]];hf[6]=h[6];

h5=h[5];

ListLogPlot[{
Abs[Gls[[;;8]]-dFdξLowList[2,h5][7]/(Range[0,7]!)]//.rules[h5]/.min42[[2]]
}]

h6=h[6];

min51=FindMinimum[chiSquared[2,h[5],(10^(-5))][7]//.rules[h[5]],Join[{#1,Re[#2]}&@@@min41[[2]],{}],MaxIterations->(10^5)]

min52=FindMinimum[chiSquared[2,h[6],(10^(-5))][7]//.rules[h[6]],Join[{#1,Re[#2]}&@@@min51[[2]],{{gC[6],0}}],MaxIterations->(10^5)]

min53=FindMinimum[chiSquared[2,h[7],(10^(-5))][7]//.rules[h[7]],Join[{#1,Re[#2]}&@@@min52[[2]],{{gC[7],0}}],MaxIterations->(10^5)]

min54=FindMinimum[chiSquared[2,h[6],(10^(-4))][7]//.rules[h[6]],Join[{#1,Re[#2]}&@@@min51[[2]],{{gC[6],0}}],MaxIterations->(10^5)]

h7=h[7];

min53=FindMinimum[chiSquared[2,h7,(10^(-13))][7]//.rules[h7],Join[{#1,Re[#2]}&@@@min52[[2]],{{gC[7],0}}],MaxIterations->(10^5)]

min52=FindMinimum[chiSquared[2,h6,(10^(-11))][7]//.rules[h6],Join[{#1,Re[#2]}&@@@min51[[2]],{}],MaxIterations->(10^5)]

h6=h[6];

min61=FindMinimum[chiSquared[2,h6,(10^(-11))][6]//.rules[h6],Join[{#1,Re[#2]}&@@@min51[[2]],{{gC[6],0}}],MaxIterations->(10^5)]

h7=h[7];

min71=FindMinimum[chiSquared[2,h7,(10^(-11))][6]//.rules[h7],Join[{#1,Re[#2]}&@@@min61[[2]],{{gC[7],0}}],MaxIterations->(10^5)]

h1'[0]

sol1s=Table[untilSuccess[newSol[randomEval[eq1],sol1,{},0,(10^(-3))]],{200}];

n2=2;
F2=RF[n2];
h2=h[2];
eq2=eqAround[n2,h2][3,-1,2];

sol2=untilSuccess[newSol[randomEval[eq2],sol1,{gC[2]},(10^(-1)),(10^(-2))]]

sol2s=Table[untilSuccess[newSol[randomEval[eq2],sol2,{},0,(10^(-3))]],{100}];

mt=8;ListLogPlot[{
Abs[Gls],
Chop[Abs[dFdξLowList[n2,h2][6]/(Range[0,6]!)]//.rules[h2]/.sol2]
}]

n3=2;
F3=RF[n3];
h3=h[4];
eq3=eqAround[n3,h3][4,-1,4];

sol3=untilSuccess[newSol[randomEval[eq3],sol2,{gC[3],gC[4]},(10^(-2))]]

sol3s=Table[untilSuccess[newSol[randomEval[eq3],sol3,{},0,(10^(-4))]],{100}];

mt=8;ListLogPlot[{
Abs[Gls],
Chop[Abs[dFdξLowList[n3,h3][6]/(Range[0,6]!)]//.rules[h3]/.sol3]
}]

ListLogPlot[{
Abs[Ghs],
Chop[Abs[dFdξHighList[n3,h3][8]/(Range[0,8]!)]//.rules[h3]/.sol3]
}]

h4=h[5];
eq4=eqAround[n3,h4][5,-1,4];

sol4=untilSuccess[newSol[randomEval[eq4],sol3,{gC[5]},(10^(-4))]]

sol4s=Table[untilSuccess[newSol[randomEval[eq4],sol4,{},0,(10^(-5))]],{100}];

h5=h[7];
eq5=eqAround[2,h5][6,-1,6];

sol5=untilSuccess[newSol[randomEval[eq5],sol4,{gC[6],gC[7]},(10^(-4))]]

sol5s=Table[untilSuccess[newSol[randomEval[eq5],sol5,{},0,(10^(-5))]],{100}];

h6=h[8];
eq6=eqAround[n3,h6][7,-1,6];

sol6=untilSuccess[newSol[randomEval[eq6],sol5,{gC[8]},(10^(-4))]]

sol6s=Table[untilSuccess[newSol[randomEval[eq6],sol6,{},0,(10^(-4))]],{100}];