TightQuad.nb

\[CapitalDelta]=85/100
m[k_]:=\[CapitalDelta]^(-2k)/(2k+1)
weights[n_]:=Table[Symbol["w"<>ToString[i]],{i,n}];
nodes[n_]:=Table[Symbol["xi"<>ToString[i]],{i,n}];
eqns[n_,p_]:=With[{w=weights[n],\[Xi]=nodes[n]},Table[Sum[w[[i]]\[Xi][[i]]^k,{i,1,n}]-m[k],{k,0,p}]];
ineqs[n_]:=With[{\[CapitalXi]=nodes[n]},Prepend[\[CapitalXi],0]-Append[\[CapitalXi],1]];
obj[n_]:=Total@Abs@weights[n];
vars[n_]:=Flatten[{nodes[n],weights[n]}];
(*Produce GAMS/BARON input*)
baron[n_,p_]:=Module[{f,eq,less,k},eq=eqns[n,p];
less=ineqs[n];
f=OpenWrite["/tmp/quad_n"<>ToString[n]<>"_p"<>ToString[p]<>".bar"];
WriteString[f,"$Offdigit\n"];
WriteString[f,"Positive Variables ",StringJoin@Riffle[ToString/@nodes[n],","],";\n"];
WriteString[f,"Variables ",StringJoin@Riffle[ToString/@weights[n],","],";\n"];
WriteString[f,"Variable norm;\n\n"];
Do[WriteString[f,ToString[v]<>".up = 1; ",ToString[v]<>".lo = 0;\n"],{v,nodes[n]}];
Do[WriteString[f,ToString[v]<>".up = 100; ",ToString[v]<>".lo = -100;\n"],{v,weights[n]}];
WriteString[f,"norm.up = 100; norm.lo = 0;\n"];
WriteString[f,"\n\nEquations ",StringJoin@Riffle[("e"<>ToString[#])&/@Range[1+Length[eq]+Length[less]],","],";\n\n"];
WriteString[f,"e1.. ",FortranForm[Total[Abs[weights[n]]]]," - norm =e= 0;\n"];
k=2;
Do[WriteString[f,"e",k,".. ",FortranForm[e]," =e= 0;\n"];k++,{e,eq}];
Do[WriteString[f,"e",k,".. ",FortranForm[e]," =l= 0;\n"];k++,{e,less}];
WriteString[f,"\nModel m / all /;\n\n"];WriteString[f,"Option dnlp = baron;\nOption decimals = 8;\nOption optcr = 0;\n\n"];
WriteString[f,"Solve m minimizing norm using dnlp;\n"];
Do[WriteString[f,"display "<>ToString@v<>".l;\n"],{v,vars[n]}];
Close[f];]
getsol[n_,p_]:=Module[{f,lines},f=OpenRead["/tmp/quad_n"<>ToString[n]<>"_p"<>ToString[p]<>".lst",BinaryFormat->True];
lines=ReadString[f];
lines=Select[StringSplit[lines,"\n"],StringMatchQ[#,"----*VARIABLE*"]&];
lines=Cases[StringSplit[#," "],Except[""]]&/@lines;
Close[f];
Thread[vars[n]->Internal`StringToDouble/@lines[[;;,-1]]]]
refinesol[n_,p_]:=Module[{approx},approx=getsol[n,p];
FindMinimum[Flatten@{obj[n],Thread[ineqs[n]<=0],Thread[eqns[n,p]==0]},Table[{v,SetPrecision[v/.approx,50]},{v,vars[n]}],WorkingPrecision->50,AccuracyGoal->20,PrecisionGoal->20]]
pretty[n_,p_]:=Module[{s,v,\[Xi],w,x,ord},v=vars[n];
s=v/.Last@Chop@N[refinesol[n,p],20];
\[Xi]=s[[;;n]];
w=s[[n+1;;2n]];
x=\[CapitalDelta] Sqrt[\[Xi]];
ord=2n;
If[\[Xi][[1]]==0,ord--;
Do[Print["\\xi_{",i,",",ord+1-i,"}&=",\[Xi][[n+1-i]],"\\\\"],{i,n-1}];
Print["\\xi_{",n,"}&=",\[Xi][[1]],"\\\\"];
Do[Print["x_{",i,",",ord+1-i,"}&=\\pm ",x[[n+1-i]],"\\\\"],{i,n-1}];
Print["x_{",n,"}&=",x[[1]],"\\\\"];
Do[Print["w_{",i,",",ord+1-i,"}&=",w[[n+1-i]],"\\\\"],{i,n-1}];
Print["w_{",n,"}&=",w[[1]],"\\\\"];,Do[Print["\\xi_{",i,",",ord+1-i,"}&=",\[Xi][[n+1-i]],"\\\\"],{i,n}];
Do[Print["x_{",i,",",ord+1-i,"}&=\\pm ",x[[n+1-i]],"\\\\"],{i,n}];
Do[Print["w_{",i,",",ord+1-i,"}&=",w[[n+1-i]],"\\\\"],{i,n}];];]
baron[5,5]
refinesol[5,5]
sol55=Last@refinesol[5,5]
Quad55[f_]:=Module[{x,w,\[Xi]},\[Xi]=Table[0,{5}];
w=Table[0,{5}];
\[Xi]=vars[5][[;;5]]/.sol55;
w=vars[5][[6;;]]/.sol55;
x=\[CapitalDelta] Sqrt[\[Xi]];
(w.f/@x)+(w.f/@-x)]
Quad55[#^12&]-Integrate[t^12,{t,-1,1}]
Quad55[#^10&]-Integrate[t^10,{t,-1,1}]