Optimal approximation of spherical squares by tensor product quadratic Bezier patches

2. Preliminaries

(Subscript[b, 0,0])={-a,-a,(Sqrt[1-2(a^2)])};(Subscript[b, 2,0])={a,-a,(Sqrt[1-2(a^2)])};(Subscript[b, 0,2])={-a,a,(Sqrt[1-2(a^2)])};(Subscript[b, 2,2])={a,a,(Sqrt[1-2(a^2)])};
(Subscript[b, 1,0])=γ((Subscript[b, 0,0])+(Subscript[b, 2,0]));(Subscript[b, 0,1])=γ((Subscript[b, 0,0])+(Subscript[b, 0,2])); (Subscript[b, 2,1])=γ((Subscript[b, 2,0])+(Subscript[b, 2,2]));(Subscript[b, 1,2])=γ((Subscript[b, 0,2])+(Subscript[b, 2,2])); 
(Subscript[b, 1,1])=δ{0,0,1};
(Subscript[B, 2,2])[u_]=(((((1+u)/2)))^2); (Subscript[B, 2,1])[u_]=2(((1+u)/2))(((1-u)/2)); (Subscript[B, 2,0])[u_]=(((((1-u)/2)))^2);
p[u_,v_,a_,γ_,δ_]=Sum[(Subscript[B, 2,i])[u](Subscript[B, 2,j])[v](Subscript[b, i,j]),{i,0,2},{j,0,2}]//Simplify;
f[u_,v_,a_,γ_,δ_]=p[u,v,a,γ,δ].p[u,v,a,γ,δ]-1//Simplify;
(Subscript[f, s])[u_,a_,γ_]=f[u,1,a,γ,0]//Simplify;
(Subscript[f, d])[u_,a_,γ_,δ_]=f[u,u,a,γ,δ]//Simplify;
g[u_,v_,a_,γ_,δ_]=(Sqrt[p[u,v,a,γ,δ].p[u,v,a,γ,δ]])-1//Simplify;
(Subscript[g, s])[u_,a_,γ_]=g[u,1,a,γ,0]//Simplify;
(Subscript[g, d])[u_,a_,γ_,δ_]=g[u,u,a,γ,δ]//Simplify;

Figure 1: The unit sphere  (light blue), boundary curves of a spherical square  (black) and a tensor product quadratic Bézier approximant  (orange) together with its control points (black).

geodesicArc[v1_List,v2_List,center_List,n_Integer]:=Module[{α0=VectorAngle[v1-center,v2-center],β,f},f[β_]=(EuclideanDistance[v1,center] Sec[1/2 (α0-2 β)] Sin[β])/EuclideanDistance[v1,v2];
Line[Table[center+Norm[v1-center]*Normalize[t*(v1-center)+(1-t) (v2-center)],{t,Map[f,Subdivide[α0,n]]}]]];

aa=1;γγ=0.61;
δδ=2.2;Show[Graphics3D[{Gray,Opacity[0.3],Sphere[{0,0,0},1]}],Graphics3D[{Black,Thickness[0.005],geodesicArc[(Subscript[b, 0,0])/.a->(1/(Sqrt[3])),(Subscript[b, 0,2])/.a->(1/(Sqrt[3])),{0,0,0},10]//N}],Graphics3D[{Black,Thickness[0.005],geodesicArc[(Subscript[b, 0,2])/.a->(1/(Sqrt[3])),(Subscript[b, 2,2])/.a->(1/(Sqrt[3])),{0,0,0},10]//N}],Graphics3D[{Black,Thickness[0.0075],geodesicArc[(Subscript[b, 2,2])/.a->(1/(Sqrt[3])),(Subscript[b, 2,0])/.a->(1/(Sqrt[3])),{0,0,0},10]//N}],Graphics3D[{Black,Thickness[0.005],geodesicArc[(Subscript[b, 2,0])/.a->(1/(Sqrt[3])),(Subscript[b, 0,0])/.a->(1/(Sqrt[3])),{0,0,0},10]//N}],
Graphics3D[{PointSize[0.025],Black,Point[{{(1/(Sqrt[3])),(1/(Sqrt[3])),(1/(Sqrt[3]))},{(1/(Sqrt[3])),-(1/(Sqrt[3])),(1/(Sqrt[3]))},{-(1/(Sqrt[3])),(1/(Sqrt[3])),(1/(Sqrt[3]))},{-(1/(Sqrt[3])),-(1/(Sqrt[3])),(1/(Sqrt[3]))}}],Point[{(Subscript[b, 1,0]),(Subscript[b, 0,1]),(Subscript[b, 1,2]),(Subscript[b, 2,1]),{0,0,1.2}}/.{a->(1/(Sqrt[3])),γ->γγ,δ->δδ}]}],
ParametricPlot3D[p[u,v,(1/(Sqrt[3])),γγ,δδ],{u,-1,1},{v,-1,1},Mesh->None,PlotStyle->{ColorData[97,"ColorList"][[2]]},PlotPoints->50],
Boxed->False]

Figure 2: Graphs of simplified radial errors  (blue) and  (orange) for “the optimal" approximant from [1] (left) and for the optimal approximant constructed in this paper (right). In all cases .  It is clearly seen that the one on the right has a smaller error on .

γ1=γ/.Solve[
D[(Subscript[f, s])[u,(1/(Sqrt[3])),γ],u]==0&&
(Subscript[f, s])[0,(1/(Sqrt[3])),γ]+(Subscript[f, s])[u,(1/(Sqrt[3])),γ]==0&&
0<u<1&&γ>0,{u,γ}][[1]];
δ1=δ/.Solve[
D[(Subscript[f, d])[u,(1/(Sqrt[3])),γ1,δ],u]==0&&
(Subscript[f, d])[0,(1/(Sqrt[3])),γ1,δ]+(Subscript[f, d])[u,(1/(Sqrt[3])),γ1,δ]==0&&
0<u<1&&δ>0,{u,δ}][[1]];

{γ2,δ2}={γ,δ}/.Solve[
D[(Subscript[f, d])[u,(1/(Sqrt[3])),γ,δ],u]==0&&
(Subscript[f, d])[0,(1/(Sqrt[3])),γ,δ]+(Subscript[f, d])[u,(1/(Sqrt[3])),γ,δ]==0&&
(Subscript[f, s])[0,(1/(Sqrt[3])),γ]==(Subscript[f, d])[0,(1/(Sqrt[3])),γ,δ]&&
0<u<1&&δ>0,{u,γ,δ}][[1]];

GraphicsRow[{Plot[{(Subscript[f, s])[u,(1/(Sqrt[3])),γ1],(Subscript[f, d])[u,(1/(Sqrt[3])),γ1,δ1]},{u,-1,1},PlotRange->{-0.023,0.023}],Plot[{(Subscript[f, s])[u,(1/(Sqrt[3])),γ2],(Subscript[f, d])[u,(1/(Sqrt[3])),γ2,δ2]},{u,-1,1},PlotRange->{-0.023,0.023}]}]

3. Main results

Lemma 3.1. For every  and  functions , , and  are strictly increasing.

D[(Subscript[f, s])[u,a,γ],γ]==(1-(a^2)) (1-(u^2)) (1+(u^2)+2γ (1-(u^2)))//Simplify

D[(Subscript[f, d])[u,a,γ,δ],γ]==(1/2)(Sqrt[1-2 (a^2)]) δ(1-(u^4))(((1-(u^2)))^2)+2γ(1-2(a^2))(((1-(u^4)))^2)+4 (a^2) γ (u^2)(((1-(u^2)))^2)+(1/2) (1-2 (a^2) )(1-(u^8))+(u^2) (1-(u^4))//Simplify

D[(Subscript[f, d])[u,a,γ,δ],δ]==(1/8) (((1-(u^2)))^2) (δ(((1-(u^2)))^2)+(Sqrt[1-2 (a^2)])(((1+(u^2)))^2)+4 (Sqrt[1-2 (a^2)]) γ(1-(u^4)))//Simplify

Lemma 3.2. Let   and let . Then  for all .

Solve[(Subscript[f, s])[0,a,γ]==(Subscript[f, d])[0,a,γ,δ],δ]//Simplify

{{δ->-(Sqrt[1-2 (a^2)])-4 (Sqrt[1-2 (a^2)]) γ-2 (Sqrt[-((-1+(a^2)) (((1+2 γ))^2))])},{δ->-(Sqrt[1-2 (a^2)])-4 (Sqrt[1-2 (a^2)]) γ+2 (Sqrt[-((-1+(a^2)) (((1+2 γ))^2))])}}

Assuming[γ>0&&0<a<(1/(Sqrt[2])),-(Sqrt[1-2 (a^2)])-4 (Sqrt[1-2 (a^2)]) γ+2 (Sqrt[-((-1+(a^2)) (((1+2 γ))^2))])==2 (1+2 γ)(Sqrt[1-(a^2)])-(1+4γ)(Sqrt[1-2 (a^2)])//Simplify]

γt=((a^2)/(2(1-(a^2))))t+(1/2);
δγt=2 (1+2 γ)(Sqrt[1-(a^2)])-(1+4γ)(Sqrt[1-2 (a^2)])/.γ->γt;

(Subscript[f, s])[u,a,γt]-(Subscript[f, d])[u,a,γt,δγt]==(u^2)(1-(u^2))(2(1-(Sqrt[1-2(a^2)])(Sqrt[1-(a^2)])-3((a^2)/2))(1-(u^2))(2-(u^2))+(a^2))(1-t)+(((u^2)(((1-(u^2)))^2))/(1-(a^2)))((1-(Sqrt[1-2(a^2)])(Sqrt[1-(a^2)])-3((a^2)/2))(4-2(u^2)-(a^2)(u^2))+(a^4)((u^2)/2))t(1-t)+(((u^2)(((1-(u^2)))^2))/(4(((1-(a^2)))^2)))((1-(Sqrt[1-2(a^2)])(Sqrt[1-(a^2)])-3((a^2)/2)-((a^4)/8))(8-4(a^2))(2-(a^2)-(u^2))+(1/2)(a^6)((a^2)+3(u^2)))(t^2)//Simplify

Equations (4)

(Subscript[f, s])[u,a,(1/2)]==-(a^2) (1-(u^2))//Simplify

(Subscript[f, s])[u,a,(1/(2-2(a^2)))]==(((a^4) (((1-(u^2)))^2))/(4(1-(a^2))))//Simplify

Lemma 3.5. If , then  and  for all .

(Subscript[δ, 0])[γ_]=2 (1+2 γ)(Sqrt[1-(a^2)])-(1+4γ)(Sqrt[1-2 (a^2)]);

(Subscript[f, s])[u,a,(1/(Sqrt[1-(a^2)]))-(1/2)]==-(((1-(Sqrt[1-(a^2)])))^2) (u^2) (1-(u^2))//Simplify

(Subscript[f, d])[u,a,(1/(Sqrt[1-(a^2)]))-(1/2),(Subscript[δ, 0])[(1/(Sqrt[1-(a^2)]))-(1/2)]]==-((2 (u^2) (1-(u^2)))/(1-(a^2)))(2(1-(Sqrt[1-(a^2)])-((a^2)/2))(u^2)( (u^2)(1-(a^2))+(a^2)(1-(u^2)))+(((1-(u^2)))^2)((((Sqrt[1-2 (a^2)])- (Sqrt[1-(a^2)])))^2)+(1-(Sqrt[1-2 (a^2)])-(a^2))(u^2)(1-(u^2))(1-(a^2)))//Simplify

Lemma 3.6. Let . Then

for all .

k=((((Sqrt[1-2 (a^2)]) -(Sqrt[1-(a^2)])))^2)(1-(u^2)) (2(((1-(u^2)))^2)+ (u^2)+ (u^4))+2 ((Sqrt[1-2 (a^2)]) (Sqrt[1-(a^2)])) (u^2)(1-(u^4))+(u^6)(1-(a^2))+(a^2) (u^2) (1-(u^2))(3-(u^2));
n=((((Sqrt[1-2 (a^2)])-(Sqrt[1-(a^2)])))^2)(2-5 (u^2)+5(u^4)-(3/2)(u^6))+ (Sqrt[1-2 (a^2)]) (Sqrt[1-(a^2)]) (u^2)+(a^2) ((u^2)-(u^4)+((u^6)/2));

D[(Subscript[f, d])[0,a,γ,(Subscript[δ, 0])[γ]],γ]-D[(Subscript[f, d])[u,a,γ,(Subscript[δ, 0])[γ]],γ]==2(u^2)(k(γ-(1/2))+n)//Simplify

Corollary 3.7. If  then .

(Subscript[f, s])[u,a,(1/(Sqrt[1-(a^2)]))-(1/2)]==-(((1-(Sqrt[1-(a^2)])))^2) (u^2) (1-(u^2))//Simplify

Lemma 3.8.  TFor the optimal parameters  and  we have f.

Lemma 3.9. Let . Let  be such that .
Then  for all . In particular, .

γz=((1-(a^2)-2 (Sqrt[1+z])(Sqrt[1-(a^2)]))/(2 (-1+(a^2))));
δz=((-(Sqrt[1-2 (a^2)])-4 (Sqrt[1+z])+4(Sqrt[1+z]) (Sqrt[1-2 (a^2)]) (Sqrt[1-(a^2)])+(a^2) ((Sqrt[1-2 (a^2)])+4 (Sqrt[1+z])))/(-1+(a^2)));
δγz=-(Sqrt[1-2 (a^2)])+4 (Sqrt[1+z])-4 (Sqrt[1-2 (a^2)]) γ;
(Subscript[A, 1])=2 (Sqrt[1-2 (a^2)]) (1-(a^2)) (((1-(u^2)))^2);
(Subscript[A, 2])=2 (1-(u^2)) (2(1-(u^2)) (Sqrt[1-2 (a^2)]) (Sqrt[1-(a^2)])+2  (u^2)(1-2(a^2))+(a^2));
(Subscript[A, 3])=2 ((Sqrt[1-2 (a^2)]) (Sqrt[1-(a^2)])2(((1-(u^2)))^2)(1-(Sqrt[1-(a^2)]))+(Sqrt[1-(a^2)])((u^4)(1-2(a^2))+(a^2)(u^2))+(1-(u^2))((1-2(a^2))2(u^2)+(a^2)));

D[(Subscript[f, d])[u,a,γz,δz],z]-D[(Subscript[f, d])[u,a,γ,δγz],z]==(((u^2) (1-(u^2))((Subscript[A, 1])((3/2)-γ)+(Subscript[A, 2])((Sqrt[1+z])-1)+(Subscript[A, 3])))/((1-(a^2))(Sqrt[1+z])))//Simplify

New parametrization with symmetric polynomials

g[e1_,e2_,γ_]=-1+(1/4) (1-(a^2)) (((1+2 γ))^2)+(1/4) ((((Sqrt[1-(a^2)])-(Sqrt[1-2(a^2)])))^2) (((1+2 γ))^2)((e1^2)-e1-e1 e2)+(1/2) (1+2 γ)(((((Sqrt[1-(a^2)])-(Sqrt[1-2(a^2)])))^2)2γ+1+(a^2)-2 γ)e2-(((((Sqrt[1-(a^2)])-(Sqrt[1-2(a^2)])))^2) (1+2 γ)γ+(a^2) γ (1-2 γ))e1 e2+(1/4) ((((Sqrt[1-(a^2)])(1+2γ)-4 (Sqrt[1-2 (a^2)]) γ))^2)(e2^2);

f[u,v,a,γ,(Subscript[δ, 0])[γ]]==g[(u^2) +(v^2), (u^2)(v^2),γ]//Simplify

has no local extrema on

has one local extrema on

{x,y}={e1,e2}/.Solve[D[g[e1,e2,γ],e1]==0&&D[g[e1,e2,γ],e2]==0,{e1,e2}][[1]]//Simplify;

x==((4 (1+2 γ) (γ (1-2 γ) (2-4 γ+(a^2) (3+2 γ))+(1+7γ+8(γ^2)-20(γ^3))((((Sqrt[1-2 (a^2)])- (Sqrt[1-(a^2)])))^2)))/((((1-2 γ))^2) (1+6 γ) ((a^2) (2 γ-1)-2 (1+2 γ)((((Sqrt[1-2 (a^2)]) -(Sqrt[1-(a^2)])))^2))))//Simplify

y==(((((1+2 γ))^2) ((a^2) (1-4 (γ^2))+2(5+2 γ-8 (γ^2))((((Sqrt[1-2 (a^2)])- (Sqrt[1-(a^2)])))^2)))/((((1-2 γ))^2) (1+6 γ) ((a^2) (2 γ-1)-2 (1+2 γ)((((Sqrt[1-2 (a^2)]) -(Sqrt[1-(a^2)])))^2))))//Simplify

(Subscript[γ, l])=((1-(1/5)(a^4))/(2 (1-(a^2))));
(Subscript[γ, r])=((1+(1/5)(a^4))/(2 (1-(a^2))));

(Subscript[f, s])[0,a,(Subscript[γ, l])]==(((a^4) (5+10 (a^2)+(a^4)))/(100 (1-(a^2))))//Simplify

For[i=0,i<=8,i++,(Subscript[η, i])={1033449984000,9650307744000,37616758528225,83055891992150,127733097043935,164350922247220,170358020376975,109671637130390,29988229774801}[[i+1]]];For[i=0,i<=4,i++,(Subscript[ζ, i])={2028600,13803300,34512660,37529100 ,14988708}[[i+1]]];
For[i=0,i<=5,i++,(Subscript[ξ, i])={2028600,14817600,41415455,53704805,30526277,5476151}[[i+1]]];
η=(a^4)Sum[(Subscript[η, i])(a^(2i))(((1-2(a^2)))^(8-i)),{i,0,8}];
ζ=Sum[(Subscript[ζ, i])(a^(2i))(((1-2(a^2)))^(4-i)),{i,0,4}];
ξ=Sum[(Subscript[ξ, i])(a^(2i))(((1-2(a^2)))^(5-i)),{i,0,5}];

(Subscript[f, d])[(3/4),a,(Subscript[γ, l]),(Subscript[δ, 0])[(Subscript[γ, l])]]+(Subscript[f, d])[0,a,(Subscript[γ, l]),(Subscript[δ, 0])[(Subscript[γ, l])]]==-(η/(6553600 (((1-(a^2)))^2)(ζ(Sqrt[1-2 (a^2)])(Sqrt[1-(a^2)])+ξ)))//Simplify

ζ=4 (10-5 (a^2)+(a^4)) (u^2) (((1-(u^2)))^2) (5(2- (a^2)-(u^2))+(a^4) (1-(u^2)));
ξ=(200 (u^2) (-2+(u^2)) (((-1+(u^2)))^2)-200 (a^2) (u^2) (((-1+(u^2)))^2) (-5+2 (u^2))+(a^10) (-2+6 (u^2)-18 (u^4)+22 (u^6)-9 (u^8))+(a^8) (22-64 (u^2)+150 (u^4)-112 (u^6)+15 (u^8))-5 (a^6) (22-62 (u^2)+134 (u^4)-110 (u^6)+27 (u^8))+5 (a^4) (18-156 (u^2)+330 (u^4)-220 (u^6)+37 (u^8)));

(Subscript[f, d])[0,a,(Subscript[γ, r]),(Subscript[δ, 0])[(Subscript[γ, r])]]+(Subscript[f, d])[u,a,(Subscript[γ, r]),(Subscript[δ, 0])[(Subscript[γ, r])]]==((ζ(Sqrt[1-2 (a^2)]) (Sqrt[1-(a^2)])+ξ)/(100(((1-(a^2)))^2)))//Simplify

For[i=0,i<=4,i++,
For[j=0,j<=9-i,j++,
(Subscript[k, i,j])={{644,1952,1236,888,9462,21754,22112,11596,3058,322},
{1646,4488,4940,11572,29861,37620,23474,7064,823},
{1522,3274,4170,11650,20601,16455,5847,761},{610,860,1350,4020,4685,2070,305},
{90,50,180 ,440,265,45}}[[i+1,j+1]]]]

ζ(1-(3/2)(a^2)-(a^4))+ξ==Sum[Sum[(Subscript[k, i,j])(u^(2j))(((1-(u^2)))^(9-i-j)),{j,0,9-i}](a^(12-2i))(((1-2(a^2)))^i),{i,0,4}]//Simplify

(Subscript[y, n])=(((1+2 γ))^2) ((a^2) (1-4 (γ^2))+2(5+2 γ-8 (γ^2))((((Sqrt[1-2 (a^2)])- (Sqrt[1-(a^2)])))^2));
(Subscript[y, d])=(((1-2 γ))^2) (1+6 γ) ((a^2) (2 γ-1)-2 (1+2 γ)((((Sqrt[1-2 (a^2)]) -(Sqrt[1-(a^2)])))^2));
y==((Subscript[y, n])/(Subscript[y, d]))//Simplify

ζ=4(-10+5 (a^2)+(a^4) (1-2 t));
ξ=40-80 (a^2)+7(a^6) (1-2 t)+(a^4)(21+8t );

(-((5(1-(a^2)))/((((1-2 γ))^2) (1+6 γ)))(Subscript[y, d])/.γ->(1-t)(Subscript[γ, l])+t (Subscript[γ, r]))==ξ+ζ(Sqrt[1-2 (a^2)]) (Sqrt[1-(a^2)])//Simplify

ξ==(40(((1-2(a^2)))^3)+160(a^2)(((1-2(a^2)))^2)+181(a^4)(1-2(a^2))+49(a^6)+8 (a^4) (1-2 (a^2)) t+2 (a^6) t)//Simplify

(ξ^2)-(((ζ(Sqrt[1-2 (a^2)]) (Sqrt[1-(a^2)])))^2)==(a^6) (17 (a^6) (((1-2 t))^2)+480 (1-(a^2))(1-t)+320t+12(a^4)t+25 (a^2)(1-2(a^2))t+105(a^2) (1-2(a^2))(1- t)+232(a^4) (1-t)+32t (a^4)(1-t))//Simplify

Remove[ξ];
For[i=0,i<=15,i++,
For[j=0,j<=4,j++,
(Subscript[ξ, i,j])={{40000,40000,40000,40000,40000},
{1050000,1050000,1050000,1050000,1050000},
{12788000,12792000,12792000,12792000,12792000},
{95809500,95885500,95885500,95885500,95885500},
{493565700,494185900,494182700,494182700,494182700},
{1850729200,1853443400,1853378200,1853378200,1853378200},
{5214530120,5220369600,5219770240,5219770080,5219770080},
{11232362060,11230397300,11227119620,11227117540,11227117540},
{18631647684,18577150372,18565319428,18565309060,18565309124},
{23771986296,23586788408,23557248632,23557228920,23557229816},
{23108818256,22747536848,22695612752,22695637840,22695643216},
{16781812960,16318939360,16254922720,16255137760,16255155680},
{8796756160,8400009920,8346057920,8346559680,8346595520},
{3135352448,2915641984,2886297216,2886907520,2886950528},
{677756672,606707456,597575424,597968640,597997312},
{66801152,56649216,55453184,55559680,55567872}}[[i+1,j+1]]]]
ξ=(1/(625(((1-(a^2)))^4)))(Sum[(Subscript[ξ, i,j])(a^(2i))(((1-2(a^2)))^(15-i))(t^j)(1-t),{i,0,15},{j,0,3}]+Sum[(Subscript[ξ, i,4])(a^(2i))(((1-2(a^2)))^(15-i))(t^4),{i,0,15}]);

k=(103901/256);
For[i=0,i<=15,i++,
For[j=0,j<=4,j++,
(Subscript[ζ, i,j])={{10447802,10447802,10447802,10447802,10447802},
{259674060,259674060,259674060,259674060,259674060},
{2987244840,2988268840,2988268840,2988268840,2988268840},
{21086895280,21104303280,21104303280,21104303280,21104303280},
{102082894880,102203266080,102202446880,102202446880,102202446880},
{358756614592,359130272192,359114502592,359114502592,359114502592},
{944866245760,944788135040,944650980480,944650939520,944650939520},
{1897672426240,1892199535360,1891487384320,1891486810880,1891486810880},
{2928265946624,2903827287552,2901370265088,2901366586880,2901366603264},
{3469548612096,3407439129088,3401508608512,3401493895680,3401494084096},
{3129407856640,3023970844672,3013710462976,3013668552704,3013669437440},
{2109908940800,1985424128000,1972627650560,1972538521600,1972540651520},
{1029768847360,927851028480,916493967360,916355686400,916358307840},
{343888683008,288232628224,281344598016,281199370240,281200549888},
{70420529152,52041220096,49451761664,49361059840,49360535552},
{6687981568,3889856512,3430318080,3405152256,3404627968}}[[i+1,j+1]]]]
ζ=(1/(625(((1-(a^2)))^4)))((1/256)(Sum[(Subscript[ζ, i,j])(a^(2i))(((1-2(a^2)))^(15-i))(t^j)(1-t),{i,0,15},{j,0,3}]+Sum[(Subscript[ζ, i,4])(a^(2i))(((1-2(a^2)))^(15-i))(t^4),{i,0,15}])-2k);

((Subscript[y, d])-(Subscript[y, n])/.γ->(1-t)(Subscript[γ, l])+t (Subscript[γ, r]))==ξ(Sqrt[1-(a^2)])(Sqrt[1-2(a^2)])-ζ//Simplify

4. Optimal (G^1) approximation

Cross[D[p[u,v,((Sqrt[3])/3),γ,δ],u],D[p[u,v,((Sqrt[3])/3),γ,δ],v]]/.{u->-1,v->-1}//Simplify

{(2/3) (1-2 γ) γ,(2/3) (1-2 γ) γ,-(4/3) (-1+γ) γ}

Solve[Cross[{(2/3) (1-2 γ) γ,(2/3) (1-2 γ) γ,-(4/3) (-1+γ) γ},{-((Sqrt[3])/3),-((Sqrt[3])/3),((Sqrt[3])/3)}]=={0,0,0},γ]

{{γ->0},{γ->(3/4)}}

(Subscript[γ, G])=(3/4);

(Cross[(D[p[u,v,((Sqrt[3])/3),(Subscript[γ, G]),δ],u]/.v->-1),(D[p[u,v,((Sqrt[3])/3),(Subscript[γ, G]),δ],v]/.v->-1)].Cross[p[0-1,-1,((Sqrt[3])/3),(Subscript[γ, G]),δ],p[1,-1,((Sqrt[3])/3),(Subscript[γ, G]),δ]])//Simplify

(1/18) (-1+(u^2)) (-7+2 (Sqrt[3]) δ)

(Subscript[δ, G])=((7(Sqrt[3]))/6);

g[0,0,((Sqrt[3])/3),(Subscript[γ, G]),(Subscript[δ, G])]//Simplify

N[g[0,0,((Sqrt[3])/3),(Subscript[γ, G]),(Subscript[δ, G])]]

Plot3D[g[u,v,((Sqrt[3])/3),(Subscript[γ, G]),(Subscript[δ, G])],{u,-1,1},{v,-1,1}]

5. Rectangular case

Clear[a,b,γ1,γ2,δ];
(Subscript[b, {0,0}])={-a,-b,(Sqrt[1-(a^2)-(b^2)])};(Subscript[b, {2,0}])={a,-b,(Sqrt[1-(a^2)-(b^2)])};(Subscript[b, {0,2}])={-a,b,(Sqrt[1-(a^2)-(b^2)])};(Subscript[b, {2,2}])={a,b,(Sqrt[1-(a^2)-(b^2)])};
(Subscript[b, {1,0}])=γ1((Subscript[b, {0,0}])+(Subscript[b, {2,0}]));(Subscript[b, {0,1}])=γ2((Subscript[b, {0,0}])+(Subscript[b, {0,2}])); (Subscript[b, {2,1}])=γ2((Subscript[b, {2,0}])+(Subscript[b, {2,2}]));(Subscript[b, {1,2}])=γ1((Subscript[b, {0,2}])+(Subscript[b, {2,2}])); 
(Subscript[b, {1,1}])=δ{0,0,1};
(Subscript[B2, 2])[u_]=(u^2); (Subscript[B2, 1])[u_]=2u (1-u); (Subscript[B2, 0])[u_]=(((1-u))^2);
p[u_,v_,a_,b_,γ1_,γ2_,δ_]=Sum[(Subscript[B2, i])[((u+1)/2)](Subscript[B2, j])[((v+1)/2)](Subscript[b, {i,j}]),{i,0,2},{j,0,2}]//Simplify;
f[u_,v_,a_,b_,γ1_,γ2_,δ_]=p[u,v,a,b,γ1,γ2,δ].p[u,v,a,b,γ1,γ2,δ]-1//Simplify;
fs1[u_,a_,γ1_]=f[u,1,a,0,γ1,0,0]//Simplify;
fs2[u_,b_,γ2_]=f[1,u,0,b,0,γ2,0]//Simplify;
fd[u_,a_,b_,γ1_,γ2_,δ_]=f[u,u,a,b,γ1,γ2,δ]//Simplify;
g[u_,v_,a_,b_,γ1_,γ2_,δ_]=(Sqrt[p[u,v,a,b,γ1,γ2,δ].p[u,v,a,b,γ1,γ2,δ]])-1//Simplify;
gs1[u_,a_,γ1_]=g[u,1,a,0,γ1,0,0]//Simplify;
gs2[u_,b_,γ2_]=g[1,u,0,b,0,γ2,0]//Simplify;
gd[u_,a_,b_,γ1_,γ2_,δ_]=g[u,u,a,b,γ1,γ2,δ]//Simplify;

Find the minimal length of the shorter edged such that we have exactly one solution

ShorterEdge[a_,h_]:=Module[{γγ1,γ1,γ2,δ,u,uu,v,b},
{γγ1,uu}={γ1,u}/.FindInstance[
D[h[u,1,a,0,γ1,0,0],u]==0&&
h[u,1,a,0,γ1,0,0]+h[0,1,a,0,γ1,0,0]==0&&0<u<1&&γ1>0,{γ1,u},WorkingPrecision->20][[1]];
b/.FindInstance[
h[0,1,a,0,γγ1,0,0]==h[1,0,0,b,0,γ2,0]&&
(D[h[u,v,a,b,γγ1,γ2,δ],v]/.{u->uu,v->1})==0&&
h[0,1,a,0,γγ1,0,0]==h[0,0,a,b,γγ1,γ2,δ]&&0<b<a&&γ2>0&&δ>0,{b,γ2,δ},WorkingPrecision->20][[1]]]

Find the only solution

OneApproximant[a_,b_,h_,P_]:=Module[{γ1,γγ1,γ2,γγ2,δ,δδ,i,u,uu,v,vv,newvalue},
{γγ1,γγ2,δδ}={γ1,γ2,δ}/.FindInstance[
h[0,1,a,0,γ1,0,0]==h[0,0,a,b,γ1,γ2,δ]&&
h[1,0,0,b,0,γ2,0]==h[0,1,a,0,γ1,0,0]&&
D[h[v,v,a,b,γ1,γ2,δ],v]==0&&
h[v,v,a,b,γ1,γ2,δ]+h[0,1,a,0,γ1,0,0]==0&&0<v<1&&γ1>0&&γ2>0&&δ>0,{γ1,γ2,δ,v},WorkingPrecision->20][[1]];
For[i=1,i<=10,i++,
{uu,vv}={u,v}/.FindInstance[D[h[u,v,a,b,γγ1,γγ2,δδ],u]==0&&D[h[u,v,a,b,γγ1,γγ2,δδ],v]==0&&0.2<=u<=1&&0.2<=v<=1,{u,v},WorkingPrecision->20][[1]];
newvalue=(h[0,1,a,0,γγ1,0,0]-h[uu,vv,a,b,γγ1,γγ2,δδ])/2;
γγ1=γ1/.FindInstance[h[0,1,a,0,γ1,0,0]==newvalue&&γ1>0,γ1,WorkingPrecision->20][[1]];
γγ2=γ2/.FindInstance[h[1,0,0,b,0,γ2,0]==newvalue&&γ2>0,γ2,WorkingPrecision->20][[1]];
δδ=δ/.FindInstance[h[uu,vv,a,b,γγ1,γγ2,δ]+newvalue==0&&δ>0,δ,WorkingPrecision->20][[1]]
];
If[P==1,
Print["Parameters of the best approximant are (Subscript[γ, 1])=",γγ1,", (Subscript[γ, 2])=",γγ2,", and δ=",δδ,"."];
Print["A minimum of the error function is at a point (",uu,",",vv,"). And the error is ",h[0,0,a,b,γγ1,γγ2,δδ],"."];
Plot3D[h[u,v,a,b, γγ1,γγ2,δδ],{u,-1,1},{v,-1,1}],
{γγ1,γγ2,δδ,uu,vv}]]
OneApproximantPlot[a_,b_,h_]:=OneApproximant[a,b,h,1]
OneApproximantParameters[a_,b_,h_]:=OneApproximant[a,b,h,0]

Find more solutions

MoreApproximant[a_,b_,h_]:=Module[{γ1,γγ1,γ2,γγ21,γγ22,γγ23,δ,δδ1,δδ2,δδ3,u,uu,v},
{γγ1,uu}={γ1,u}/.FindInstance[
D[h[u,1,a,0,γ1,0,0],u]==0&&
h[u,1,a,0,γ1,0,0]+h[0,1,a,0,γ1,0,0]==0&&0<u<1&&γ1>0,{γ1,u},WorkingPrecision->20][[1]];
{γγ21,δδ1}={γ2,δ}/.FindInstance[
h[0,1,a,0,γγ1,0,0]==h[1,0,0,b,0,γ2,0]&&
h[0,1,a,0,γγ1,0,0]==h[0,0,a,b,γγ1,γ2,δ]&&γ2>0&&δ>0,{γ2,δ},WorkingPrecision->20][[1]];
{γγ22,δδ2}={γ2,δ}/.FindInstance[
h[0,1,a,0,γγ1,0,0]==h[1,0,0,b,0,γ2,0]&&
(D[h[u,v,a,b,γγ1,γ2,δ],v]==0/.{u->uu,v->1})&&γ2>0&&δ>0,{γ2,δ},WorkingPrecision->20][[1]];
{γγ23,δδ3}={γ2,δ}/.FindInstance[
h[0,1,a,0,γγ1,0,0]==h[0,0,a,b,γγ1,γ2,δ]&&
(D[h[u,v,a,b,γγ1,γ2,δ],v]==0/.{u->uu,v->1})&&γ2>0&&δ>0,{γ2,δ},WorkingPrecision->20][[1]];
Print["The parameter (Subscript[γ, 1])=", γγ1," and for the other two parameters we have more possibilities."]; Print["Every pair ((Subscript[γ, 2]),δ) in the triangle with vertices (",N[γγ21,4],",",N[δδ1,4],"), (",N[γγ22,4],",",N[δδ2,4],"), (",N[γγ23,4],",",N[δδ3,4],") induces a best approximant." ];
GraphicsRow[{
Plot3D[h[u,v,a,b, γγ1,γγ21,δδ1],{u,-1,1},{v,-1,1}],
Plot3D[h[u,v,a,b, γγ1,γγ22,δδ2],{u,-1,1},{v,-1,1}],
Plot3D[h[u,v,a,b, γγ1,γγ23,δδ3],{u,-1,1},{v,-1,1}]}]]

Figure 3 : The domain (light blue) in   for which pairs  imply several optimal approximants of the spherical rectangle . The grey dashed curve is the graph of the function  and together with the -axis determine the admissible domain of pairs .

points={{0,0}};
For[aa=0.02,aa<=0.99,aa+=0.01,AppendTo[points,{aa,ShorterEdge[aa,g]}]]
For[aa=0.991,aa<=0.999,aa+=0.001,AppendTo[points,{aa,ShorterEdge[aa,g]}]]
AppendTo[points,{1,0}]

Show[Plot[Min[x,Sqrt[1-(x^2)]],{x,0,1},AspectRatio->Automatic,PlotStyle->{Gray,Thin}],ListLinePlot[points,Filling->Axis]]

BestApproximant[a_,b_,h_]:=Module[{aa=Max[a,b],bb=Min[a,b],ba},
If[(aa^2)+(bb^2)>=1,Print["No approximant"],
ba=ShorterEdge[aa,h];
If[bb>=ba,
OneApproximantPlot[aa,bb,h],
MoreApproximant[aa,bb,h]]
]]

Figure 4: Graphs of  for three optimal approximants of the spherical rectangle given by parameters  and . Here  and the left graph corresponds to the pair , the middle one to , and the right one to .

BestApproximant[0.8,0.6,g]

BestApproximant[0.5,0.5,g]

MoreApproximant[0.75,0.2,g]

6. Numerical examples

(Subscript[a, max])=(1/(Sqrt[2]));

tabel={{"a","((Subscript[γ, r]^( * )),(Subscript[δ, r]^( * )))","|g|","r"}};
x=OneApproximantParameters[N[(Subscript[a, max])],N[(Subscript[a, max])],g];
d1=Abs[g[0,0,(Subscript[a, max]),(Subscript[a, max]),x[[1]],x[[1]],x[[3]]]];
AppendTo[tabel,
{"(Subscript[a, max])","("<>ToString[N[x[[1]],5]]<>","<>ToString[N[x[[3]],5]]<>")",ScientificForm[d1],"-"}];
For[i=1,i<=6,i++,
x=OneApproximantParameters[N[(1/(2^i))(Subscript[a, max])],N[(1/(2^i))(Subscript[a, max])],g];
d2=Abs[g[0,0,(1/(2^i))(Subscript[a, max]),(1/(2^i))(Subscript[a, max]),x[[1]],x[[1]],x[[3]]]];
AppendTo[tabel,
{1/(2^i)"(Subscript[a, max])","("<>ToString[N[x[[1]],5]]<>","<>ToString[N[x[[3]],5]]<>")",ScientificForm[d2],Log[d1/d2]/Log[2]}];
d1=d2]

Table 1: Optimal parameters  and  according to the radial error , the corresponding radial distance  and the numerical rate of convergence  for several parameters , , with .

Grid[tabel,Frame->All]