Test při výpočtu klasického Steinerova bodu

Pomocí tohoto souboru byly generovány grafy a data k příslušným testům numerických integračních metod.

(Při vykonávání spouštět sekvenčně - inicializace funkcí, pro test 1 -> test1 ; init pro test 2 -> test 2 ; init pro test 3 -> test 3)

Test 1

Pomocne funkce (operna funkce aj.)

supportUpper[x_,y_,R_]:= R;
supportLower[x_,y_,R_] := R * Abs[x];
support[x_,y_,R_] := If[y≥0,supportUpper[x,y,R],supportLower[x,y,R]];
supportTimesY[x_,y_,R_]:=y * support[x,y,R];
hAe[{x_,y_}] := supportTimesY[x,y,1];
points[num_] := Table[{Cos[t],Sin[t]},{t,0,2*Pi,2*Pi/num}][[;;-2]];

Table[If[Mod[i,2]==0,2,4],{i,0,3}]

Metody tablkou - Simpson, Lichobeznik

trapezoid[num_] := (  (2*Pi)/num)*          Total[                                                                                                      Map[hAe,points[num]] ]/(Pi);
simpson[num_] :=      (((2*Pi)/num)/3)* Total[Table[If[Mod[i,2]==0,2,4],{i,0,num-1}] * Map[hAe,points[num]] ]/(Pi);

Rondbergova metoda

Test

compSteinerRondberg[minPoints_,maxPoints_] :=
Module[{height,table, layerPoints,i,j},
height = Log[2,maxPoints] - Log[2,minPoints] +1;
table = Array[0 &,{height,height}];

(* vypln prvni sloupce *)
For[i=1,i≤ height,i++,
table[[i,1]] = trapezoid[(2^((Log[2,minPoints]+i-1)))] ;
];

(* vypln zbytek tabulky *)
For[j=2,j≤ height, j++,
For[i=j,i≤ height, i++,
table[[i,j]] = table[[i,j-1]] + (((table[[i,j-1]] - table[[i-1,j-1]]))/((4^((j-1)))-1))
]
];
Return[table];
];

Module[{resRom,resTrap,resSimp,i,realRes,tableRom},
resTrap= {};resSimp={};resRom={};
realRes = 1/Pi;
For[i=32,i≤512,i=2*i,
AppendTo[resSimp,simpson[i]];
];
tableRom = compSteinerRondberg[32,512];
resRom = Diagonal[tableRom];

resPlot = ListPlot[{resSimp[[2;;]],resRom[[2;;]]},PlotRange->All,PlotMarkers->{Automatic,Small}];
errPlot =ListPlot[{Abs[resSimp[[1;;]]-realRes],Abs[resRom[[1;;]]-realRes]},PlotMarkers->{Automatic,Small}];
Return[N[{resSimp,resRom},20]];
];

tcks = {Table[{i,2^(i+4)},{i,1,8}],Automatic}

err = Show[resPlot,Ticks->tcks,AxesLabel->{N,"log(d)"}]

err= Show[errPlot,Ticks->tcks,AxesLabel->{N,"log(d)"}]

Table[N[simpson[i]-1/Pi],{i,{32,64,128,256,512}}]

TableForm[N[compSteinerRondberg[32,512]]]

Test 2

Funkce

PRECISION = 50;

shift[{x_,y_}] := If[y≥0,{x,y},{x,0}]

pointsA[num_] := Join[
(*Map[shift,Table[{Cos[t],Sin[t]},{t,0,2*Pi,2*Pi/num}][[;;-2]]][[;;-(1+(num-2)/2)]],*)
Map[shift,Table[{Cos[t],Sin[t]},{t,0,2*Pi,2*Pi/num}][[;;-2]]][[;;-(1+(num-2)/2)]],
Table[{-1+i,0},{i,0,2,2/(1+(num-2)/2)}][[2;;-2]]
];

pointsShift[num_] := Table[{Cos[t+(Pi/(Sqrt[2]))],Sin[t+(Pi/(Sqrt[2]))]},{t,0,2*Pi,2*Pi/num}][[;;-2]];

findSupport[{ex_,ey_},pts_] := Max[Map[({ex,ey} . #) &, pts]]

findSupportPoint[{ex_,ey_},pts_] := 
Module[{maxIndex,supports},
supports= Map[({ex,ey} . #) &, pts] ;
maxIndex = Position[supports,Max[supports],1][[1,1]];
Return[pts[[maxIndex]]];
];

hAe2[{ex_,ey_},pts_] := ey* findSupport[{ex,ey},pts];
points[num_] := Table[{Cos[t],Sin[t]},{t,0,2*Pi,2*Pi/num}][[;;-2]];

Ntrapezoid2shift[num_] := ((2*Pi)/num)*Total[ParallelMap[Quiet[hAe2[#,N[pointsA[num],PRECISION]]] &,N[pointsShift[num],PRECISION]]]/(Pi);

Nsimpson2shift[num_] := (((2*Pi)/num)/3)* Total[Table[If[Mod[i,2]==0,2,4],{i,0,num-1}] * ParallelMap[(hAe2[ #,N[pointsA[num],PRECISION]]) &,N[pointsShift[num],PRECISION]]] / (Pi);

Nromberg2shift[minPoints_,maxPoints_] :=
Module[{height,table, layerPoints,i,j},
height = Log[2,maxPoints] - Log[2,minPoints] +1;
table = Array[0 &,{height,height}];

(* vypln prvni sloupce *)
For[i=1,i≤ height,i++,table[[i,1]] = Ntrapezoid2shift[(2^((Log[2,minPoints]+i-1)))] ;];

(* vypln zbytek tabulky *)
For[j=2,j≤ height, j++,
For[i=j,i≤ height, i++,
table[[i,j]] = table[[i,j-1]] + (((table[[i,j-1]] - table[[i-1,j-1]]))/((4^((j-1)))-1))
]];
Return[table];
];

Obrázky

aproxA = ListPlot[pointsA[32],AspectRatio->Automatic,PlotRange->{Automatic,{-1.1,1.1}},PlotMarkers->{Automatic,Tiny}]

aproxDir = ListPlot[pointsShift[32],AspectRatio->Automatic, PlotRange->{Automatic,{-1.1,1.1}},PlotMarkers->{Automatic,Tiny}]

ListPlot[{pointsA[32],pointsShift[32]},AspectRatio->Automatic,PlotMarkers->{Automatic,Tiny}]

TEST

Module[{resRom,resTrap,resSimp,i,realRes,tableRom},
resSimp={};resRom={};
realRes = 1/Pi;
For[i=32,i≤512,i=2*i,AppendTo[resSimp,Nsimpson2shift[i]];Print["Bodů Simpson: ",i];];
tableRom = Nromberg2shift[32,512];
resRom = Diagonal[tableRom];

resPlot = ListPlot[{resSimp[[1;;]],resRom[[1;;]]},PlotRange->All,PlotMarkers->{Automatic,Small}];
errPlot =ListPlot[{Abs[resSimp[[1;;]]-realRes],Abs[resRom[[1;;]]-realRes]},PlotMarkers->{Automatic,Small}];

Return[N[{resSimp,resRom}]];
];

tcks = {Table[{i,2^(i+4)},{i,1,8}],Automatic}

Show[errPlot,AxesLabel->{N,"d"},Ticks->tcks,PlotRange->All]

Show[resPlot,AxesLabel->{N,"d"},Ticks->tcks,PlotRange->All]

Table[N[Nsimpson2shift[i]-1/Pi],{i,{32,64,128,256,512,1024}}]

Diagonal[N[Nromberg2shift[32,4096]-1/Pi]]

TableForm[N[Nromberg2shift[32,1024]]]

Test 3

Definice

INTERVAL = {-(10^(-4)),(10^(-4))};
PRECISION = 50;

points[num_] := Table[{Cos[t],Sin[t]},{t,0,2*Pi,2*Pi/num}][[;;-2]];

pointsAfuzzy[num_] := Join[
(*Map[shift,Table[{Cos[t],Sin[t]},{t,0,2*Pi,2*Pi/num}][[;;-2]]][[;;-(1+(num-2)/2)]],*)
Map[shift,Table[(1+RandomReal[INTERVAL])*{Cos[t],Sin[t]},{t,0,2*Pi,2*Pi/num}][[;;-2]]][[;;-(1+(num-2)/2)]],
Table[{-1+i,RandomReal[INTERVAL]},{i,0,2,2/(1+(num-2)/2)}][[2;;-2]]
];

shift[{x_,y_}] := If[y≥0,{x,y},{x,0}]

pointsShift[num_] := Table[{Cos[t+(Pi/(Sqrt[2]))],Sin[t+(Pi/(Sqrt[2]))]},{t,0,2*Pi,2*Pi/num}][[;;-2]];

findSupport[{ex_,ey_},pts_] := Max[Map[({ex,ey} . #) &, pts]]

hAe2[{ex_,ey_},pts_] := ey* findSupport[{ex,ey},pts];

ListPlot[pointsAfuzzy[32],AspectRatio->Automatic]

N[pointsShift[4]]

pointsAfuzzy[4]

Nsimpson3shift[4,{{0.9999538387080701`,0.`},{0.`,1.0000644618970964`},{-0.9999315776899205`,0.`},{0,4.9942263720838655`*^-6}}]

Obrázky

ListPlot[{pointsA[32],pointsAfuzzy[32]},AspectRatio->Automatic,PlotMarkers->{Automatic,Tiny}]

a = Graphics[{Circle[{0,0},1,{0,Pi}], Line[{{-1,0},{1,0}}]},PlotRange->{{-1.2,1.2},{-0.2,1.2}}]

b = Graphics[{Blue,PointSize[Large],Point[pointsA[16]], Gray,Dashed,Circle[{0,0},1,{0,Pi}], Line[{{-1,0},{1,0} }]},PlotRange->{{-1.2,1.2},{-0.2,1.2}}]

im = Graphics[{
PointSize[Large],
Orange,
Point[pointsAfuzzy[16]],
Gray,
Arrowheads[Small],
Arrow[Table[{{-1+i,0},{-1+i,0.15}},{i,0,2,2/(1+(16-2)/2)}][[2;;-2]]],
Arrow[Table[{{-1+i,0},{-1+i,-0.15}},{i,0,2,2/(1+(16-2)/2)}][[2;;-2]]],
Arrow[Table[{{Cos[t],Sin[t]},1.15{Cos[t],Sin[t]}},{t,0,Pi,(2Pi)/16}]],
Arrow[Table[{{Cos[t],Sin[t]},0.85{Cos[t],Sin[t]}},{t,0,Pi,(2Pi)/16}]],
Dashed,
Circle[{0,0},1,{0,Pi}],
Line[{{-1,0},{1,0}}]
},PlotRange->{{-1.2,1.2},{-0.2,1.2}}]

Export["test3-shift-a.pdf",a,"AllowRasterization"->True,ImageSize->250,ImageResolution->600];
Export["test3-shift-b.pdf",b,"AllowRasterization"->True,ImageSize->250,ImageResolution->600];
Export["test3-shift-c.pdf",im,"AllowRasterization"->True,ImageSize->250,ImageResolution->600];

Funkce

Nsimpson3shift[num_,pts_] := (((2*Pi)/num)/3)* Total[Table[If[Mod[i,2]==0,2,4],{i,0,num-1}] * ParallelMap[(hAe2[ #,N[pts,PRECISION]]) &,N[pointsShift[num],PRECISION]]] / (Pi);

Ntrapezoid3shift[num_,pts_] := ((2*Pi)/num)*Total[ParallelMap[Quiet[hAe2[#,N[pts[[;;;;Length[pts]/num]],PRECISION]]] &,N[pointsShift[num],PRECISION]]]/(Pi);

Nromberg3shift[minPoints_,maxPoints_,pts_] :=
Module[{height,table, layerPoints,i,j},
height = Log[2,maxPoints] - Log[2,minPoints] +1;
table = Array[0 &,{height,height}];

(* vypln prvni sloupce *)
For[i=1,i≤ height,i++,table[[i,1]] = Ntrapezoid3shift[(2^((Log[2,minPoints]+i-1))),pts] ;];

(* vypln zbytek tabulky *)
For[j=2,j≤ height, j++,
For[i=j,i≤ height, i++,
table[[i,j]] = table[[i,j-1]] + (((table[[i,j-1]] - table[[i-1,j-1]]))/((4^((j-1)))-1))
]];
Return[table];
];

Test

Module[{resRom,resTrap,resSimp,i,realRes,tableRom, pts,j,k},
resSimp=Array[0 &,4];resRom=Array[0 &,4];
realRes = 1/Pi;
k=1;
For[j=1,j≤20,j++,
pts = pointsAfuzzy[256];
For[i=32,i≤256,i=2*i,resSimp[[k]] += Nsimpson3shift[i,pts[[;;;;Length[pts]/i]] ]; k++;];
resRom += Diagonal[Nromberg3shift[32,256,pts]];
k=1;
];

resRom /= 20; resSimp /= 20;

resPlot = ListPlot[{resSimp[[1;;]],resRom[[1;;]]},PlotRange->All,PlotMarkers->{Automatic,Small}];
errPlot =ListPlot[{Abs[resSimp[[1;;]]-realRes],Abs[resRom[[1;;]]-realRes]},PlotRange->All,PlotMarkers->{Automatic,Small}];

Print["S:",resSimp-1/Pi];
Print["R:",resRom-1/Pi];
Return[{resSimp,resRom}];
];

tcks = {Table[{i,2^(i+4)},{i,1,8}],Automatic}

Show[errPlot,PlotRange->All,Ticks->tcks,AxesLabel->{N,"d"}]

Show[resPlot,PlotRange->All,Ticks->tcks,AxesLabel->{N,Subscript["s","y"]},GridLines->{None,{1/Pi}}]

Label

gg= Graphics[{Text[Style["Simpsonova metoda",FontSize->15],{0.55,0}],Inset[Graphics[{Blue,Disk[{0,0},1] },ImageSize->{15,15}],{-0.6,0}],
Text[Style["Rombergova metoda",FontSize->15],{2.5,0}],Inset[ Graphics[{Orange,Rectangle[{1,0},{2,1}]},ImageSize->{15,15}],{1.4,0} ]},PlotRange->All];

Show[gg,PlotRangePadding->0]

gg=

Export["test-label.pdf",gg,"AllowRasterization"->True,ImageResolution->600];