INTERVALOVÁ ARITMETIKA A VÝZKUM DIFERENCIÁLNÍCH ROVNIC

Jiří Benedikt

Katedra matematiky, Západočeská univerzita v Plzni

Author je podporován grantem č. ME 10093.

Aritmetika ve strojové přesnosti (IEEE 754)

Uvažujme následující příklad. Chceme ověřit, že hodnota polynomu v jeho kořenu je nulová.

Expand[(x-4)^3]

%/.x->4.

Stejný výsledek dostaneme pro libovolnou násobnost kořene, např. 30.

Expand[(x-4)^30]/.x->4.

Nyní vyzkoušejme stejný výraz, ale s kořenem 3.

Expand[(x-3)^30]/.x->3.

Proč nevyšla nula? Obecně řečeno je problém v tom, že výsledky aritmetických operací ve strojové přesnosti jsou zaokrouhlovány, protože na počítači lze reprezentovat pouze konečný počet čísel. Neočekávaný výsledek našeho výpočtu vznikl proto, že vysoké mocniny 4 jsou přesně reprezentovatelné v IEEE 754, zatímco vysoké mocniny 3 nejsou.

Intervalová aritmetika

Použitím intervalové aritmetiky v Mathematice dostaneme vždy pravdivou odpověď.

Expand[(x-3)^30]/.x->Interval[3.]

Tento výsledek znamená, že přesný výsledek 0 musí ležet v uvedeném intervalu. Samozřejmě se může stát, že výsledný interval obsahuje 0 i když nedosazujeme kořen:

Expand[(x-3)^30]/.x->Interval[4.]

V takovém případě lze použít intervalovou aritmetiku s vyšší přesností.

Expand[(x-3)^30]/.x->Interval[4.`30]

Tento výsledek zaručuje, že 4 není kořen polynomu (((x-3))^30). Intervalová aritmetika obecně nezaručuje, jaký je přesný výsledek, ale jaký není.

Vlastní funkce p-laplaciánu na kruhu

Zabýváme se vlastními funkcemi Dirichletovy okrajové úlohy

-(Subscript[Δ, p])u=λ|u((|)^(p-2))u              v D
         u=0                                 na ∂D

kde D⊂(ℝ^2) jednotkový otevřený kruh se středem v počátku, (Subscript[Δ, p])u=div(|∇u((|)^(p-2))∇u) je p-laplacián, p>1 a λ∈ℝ je spektrální parametr. Pro p=2 se rovnice v (1) redukuje na lineární rovnici Δu=λu. Říkáme, že λ je vlastní číslo (1), pokud má (1) nenulové řešení, které pak nazýváme příslušnou vlastní funkcí. Je známo, že první vlastní funkce (příslušná prvnímu vlastnímu číslu) nemění znaménko a že je radiálně symetrická.

Plot3D[BesselJ[0,Sqrt[x^2+y^2]*BesselJZero[0,1]],{x,-1,1},{y,-1,1},RegionFunction->Function[{x,y,z},x^2+y^2≤1]]

V článku [Benedikt, J., Drábek, P., Girg, P., The second eigenfunction of the p-Laplacian on the disc is not radial, submitted, 2010] jsme dokázali, že druhá radiálně symetrická vlastní funkce (1)

Plot3D[BesselJ[0,Sqrt[x^2+y^2]*BesselJZero[0,2]],{x,-1,1},{y,-1,1},RegionFunction->Function[{x,y,z},x^2+y^2≤1]]

není druhá vlastní funkce (1), tj. že druhá vlastní funkce (1) není radiálně symetrická.

p-Besselova funkce

Radiálně symetrické vlastní funkce (1) dostaneme rotací řešení počáteční úlohy

-(r|u'((|)^(p-2))u')'=r|u((|)^(p-2))u,   r>0,
u(0)=1,    u'(0)=0.

Řesení značíme (Subscript[J, 0,p]) a nazýváme ho p-Besselovou funkcí, protože pro p=2 je (Subscript[J, 0,2])=(Subscript[J, 0]) klasická Besselova funkce prvního druhu.

Použitím jistých variačních argumentů jsme v [Benedikt, J., Drábek, P., Girg, P., The second eigenfunction of the p-Laplacian on the disc is not radial, submitted, 2010] dokázali, že pokud první dva nulové body (Subscript[ν, 1,])(Subscript[ν, 2]) funkce (Subscript[J, 0,p]) splňují nerovnost (Subscript[ν, 2])>2(Subscript[ν, 1]), pak druhá vlastní funkce (1) nemůže být radiálně symetrická. Pro p=2 jsou první dva kořeny

Table[N[BesselJZero[0,n]],{n,1,2}]

Nerovnost (Subscript[ν, 2])>2(Subscript[ν, 1]) pro p blízko 1 a blízko ∞ jsme dokázali asymptotickou analýzou, zatímco pro p∈[1.01,226] jsme použili intervalovou aritmetiku v Mathematice.

Intervalová numerická integrace počáteční úlohy

Uvažujme obecnou počáteční úlohu pro rovnici (soustavu rovnic) prvního řádu

u'(r)=f(r,u),
u(0)=(Subscript[u, 0]).

Nechť F je tzv. intervalové rozšíření f, tj. pro dané intervaly R a U, F(R,U) je interval splňující r ∈ R, u ∈ U ⟹ f (r, u) ∈ F (R, U). Mějme h>0. Nejdříve určíme interval (co nejmenší), ve kterém musí ležet u(r) pro všechna r∈(Subscript[R, 0])=[0,h]. K tomu potřebujeme a priorní odhad u. Pomocí první integrálu rovnice lze např. snadno dokázat, že hodnoty p-Besselovy funkce leží v intervalu (Subscript[I, 0])=[-1,1] pro všechna r≥0. Vlastnost F pak zaručuje, že u'(r)∈F((Subscript[R, 0]),(Subscript[I, 0])) pro všechna r∈(Subscript[R, 0]). Tudíž u(r)=(Subscript[u, 0])+(∫_0^r)u'(s)ds∈(Subscript[u, 0])+(Subscript[R, 0])·F((Subscript[R, 0]),(Subscript[I, 0])) pro všechna r∈(Subscript[R, 0]). Nyní víme, že u(r)∈(Subscript[I, 1])=((Subscript[u, 0])+(Subscript[R, 0])·F((Subscript[R, 0]),(Subscript[I, 0])))⋂(Subscript[I, 0]) pro r∈(Subscript[R, 0]), kde interval (Subscript[I, 1]) je pravděpodobně menší, než (Subscript[I, 0]) (maximálně stejný). Odtud u'(r)∈F((Subscript[R, 0]),(Subscript[I, 1])) a celý postup můžeme opakovat, až dostaneme u(r)∈(Subscript[I, 2])=((Subscript[u, 0])+(Subscript[R, 0])·F((Subscript[R, 0]),(Subscript[I, 1])))⋂(Subscript[I, 1]). Podle naší zkušenosti jsou tři iterace ((Subscript[I, 3])) v jistém smyslu optimální. Pak pokračujeme analogicky na (Subscript[R, 1])=[h,2h] s počáteční podmínkou u(h)∈(Subscript[u, 0])+h·F((Subscript[R, 0]),(Subscript[I, 2])) a tak dále, až pokryjeme dostatečně velký interval. Pro dostatečně malý krok h dostaneme takové intervaly pro (Subscript[ν, 1,])(Subscript[ν, 2]) (ve kterých (Subscript[ν, 1,])(Subscript[ν, 2]) musí ležet), že levý koncový bod intervalu pro (Subscript[ν, 2]) je ostře větší, než pravý koncový bod intervalu pro (Subscript[ν, 1]). To dokazuje, že (Subscript[ν, 2])>2(Subscript[ν, 1]).

Implementace

Manipulate[Block[{pa,dim,step,IntOb,It,dr,Mezi,s,c,sol,fi,ad},
fi[pp_?NumberQ,s_?NumberQ]=Abs[s]^(pp-1)*Sign[s];
fi[pp_?NumberQ,Interval[{sl_?NumberQ,sh_?NumberQ}]]=Interval[{fi[pp,sl],fi[pp,sh]}];
fi[Interval[{ppl_?NumberQ,pph_?NumberQ}],Interval[{sl_?NumberQ,sh_?NumberQ}]]=Interval[{fi[If[sl>-1&&sl<0||sl>1,ppl,pph],sl],fi[If[sh>-1&&sh<0||sh>1,pph,ppl],sh]}];
ad[pp_?NumberQ]=pp/(pp-1);
ad[Interval[{ppl_?NumberQ,pph_?NumberQ}]]=Interval[{pph/(pph-1),ppl/(ppl-1)}];
pa=ad[p];
dim=2;
step=length/steps;
IntOb[i_]:=Interval[{Min[Union@@i],Max[Union@@i]}];
It[{u_,v_,r_}]:=Nest[{IntervalIntersection[Interval[{-1,1}],u+fi[pa,IntOb[{v,#[[2]]}]/Interval[{r,#[[3]]}]^(dim-1)]*step],IntervalIntersection[Interval[{-1,1}*#[[3]]^dim/dim],v-Interval[{r,#[[3]]}]^(dim-1)*fi[p,IntOb[{u,#[[1]]}]]*step],#[[3]]}&,{Interval[{-1,1}],Interval[{-1,1}*(r+step)^dim/dim],r+step},iterations];
dr=Prepend[NestList[It,{Interval[{1.-step^pa/pa/dim^(1/(p-1)),1.}],Interval[{-step^dim/dim*1.,0.}],step},steps-1],{Interval[1.],Interval[0.],0}];
Mezi[l_]:=IntOb/@Transpose[{Drop[l,-1],Drop[l,1]}];
s=Mezi[Transpose[dr][[1]]];
c=Mezi[Transpose[dr][[2]]];
sol={u,v}/.NDSolve[{v'[r]==-r^(dim-1)*fi[p,u[r]],u'[r]==fi[p/(p-1),v[r]/r^(dim-1)],u[10^-6]==1,v[10^-6]==0},{u,v},{r,10^-6,length}][[1]];
Show[Graphics[{Green,Rectangle[{#[[1]]*step,Min[#[[2]]]},{(#[[1]]+1)*step,Max[#[[2]]]}]&/@Transpose[{Table[i,{i,0,steps-1}],s}]},PlotRange->{{0,length},{-1.1,1.1}},Axes->True],Plot[sol[[1]][r],{r,0,length}]]],{{p,2},1.1,4},{{length,8},1,10},{{steps,2000},1,3000,1},{{iterations,3},0,5,1,Setter}]