Computing Principal Eigenvalue in Interval Arithmetic

Jiří Benedikt, University of West Bohemia, Plzeň, Czech Republic

ABSTRACT

The solvability of boundary value problems for elliptic and parabolic differential equations (e.g., diffusion equation) is closely related to the corresponding eigenvalues. In particular, the principal eigenvalue in which the geometry of the energy functional is substantially changed. In the lecture we present a numerical method of finding guaranteed upper and lower estimates of the principal eigenvalue in Mathematica using interval arithmetic.

Interval arithmetic

The number 3 is a root of the following polynomials:

Expand[(x-3)^3]

%/.x->3.

pol[x_]=Expand[(x-3)^30];
pol/@{1.,2.,3.,4.}

pol/@{1.`20,2.`20,3.`20,4.`20}

Interval[3.]
%//InputForm

Interval[N[Pi]]//InputForm

Interval[2.]*Interval[3.]//InputForm

pol/@{Interval[1.],Interval[2.],Interval[3.],Interval[4.]}

pol/@{Interval[2.`30],Interval[3.`30],Interval[4.`30]}

Eigenvalues of boundary value problems

Eigenfunctions and eigenvalues of a Dirichlet boundary value problem of ordinary second-order differential equation with constant coefficient:

eqncc=-u''[t]==λ u[t];
bc={u[0]==0,u[1]==0};
DSolveValue[{eqncc,bc},u[t],t]

The principal (first) eigenfunction:

pEigf[t_]=Sin[Pi t]
-pEigf''[t]

The principal eigenvalue is (π^2).

N[Pi^2]

Plot[pEigf[t],{t,0,1}]

Example with a non-constant coefficient:

w[t_]=t+Cos[t];
eqnnc=-u''[t]==λ w[t]u[t];
DSolveValue[{eqnnc,bc},u[t],t]

Upper estimate of the principal eigenvalue

Variational characterization of the principal eigenvalue for -u''(t)=λ w(t) u(t):

(Subscript[λ, 1])=min(( (∫_0^1)v'(((t))^2)dt)/((∫_0^1)w(t)((v(t))^2)dt))

where the minimum of the so-called Rayleigh quotient is taken over all non-zero v∈(Subscript[W, 0]^(1,2))(0,1). It is achieved when v is the principal eigenfunction:

{Integrate[pEigf'[t]^2,{t,0,1}],Integrate[pEigf[t]^2,{t,0,1}]}

Using a similar test function:

v[t_]=4t(1-t);
Plot[{v[t],pEigf[t]},{t,0,1}]

Upper estimate for the constant coefficent (exact value is (π^2)≐9.8696):

Integrate[v'[t]^2,{t,0,1}]/Integrate[v[t]^2,{t,0,1}]

Upper estimate for the coefficient w(t)=t+cos(t):

FullSimplify[Integrate[v'[t]^2,{t,0,1}]/Integrate[w[t]v[t]^2,{t,0,1}]]
upper=N[%]

Lower estimate of the principal eigenvalue

The well known Picone identity implies that for any test function v(t), positive and (C^2) on (0,1),

(Subscript[λ, 1])≥(inf_(t∈(0,1)))(( -v''(t))/(w(t)v(t))).

Similarly, equality holds when v is the principal eigenfunction. Using again a similar test function:

Plot[-v''[t]/v[t],{t,0,1}]

the lower estimate for the constant coefficent is

FindMinimum[{-v''[t]/v[t],0<t<1},t]

Upper estimate for the coefficient w(t)=t+cos(t):

Plot[-v''[t]/(w[t]v[t]),{t,0,1}]

FindMinimum[{-v''[t]/(w[t]v[t]),0<t<1},t]
lower=%[[1]];

Hence, the principal eigenvalue for the non-constant coefficient w(t)=t+cos(t) is in

Interval[{lower,upper}]

How to narrow the interval? We need a test function which is more similar to the principal eigenfunction.

Constructing a test function from the numerical solution

vao[λ_?NumericQ]:=NDSolveValue[{-u''[t]==λ w[t]u[t],u[0]==0,u'[0]==1},u[t],{t,0,1}]/.t->1

(Subscript[λ, e])=λ/.FindRoot[vao[λ]==0,{λ,upper}];
{(Subscript[λ, e])-lower,upper-(Subscript[λ, e])}
(Subscript[λ, e])//InputForm

pEigfN[t_]=NDSolveValue[{-u''[t]==(Subscript[λ, e]) w[t]u[t],u[0]==0,u'[0]==1},u[t],{t,0,1}];
Plot[pEigfN[t],{t,0,1}]

To compute both the upper and the lower estimate, we need a (C^2) function, positive in (0,1), zero at 0 and 1. We construct it as a cubic spline interpolation function given by the control points:

cp[n_?IntegerQ]:=Flatten[{0,Table[pEigfN[t/n],{t,1,n-1}],0}]
cp[4]

First we need the (n-1)×(n-1) matrix

A[n_?IntegerQ]:=DiagonalMatrix[Table[4,n-1]]+DiagonalMatrix[Table[1,n-2],1]+DiagonalMatrix[Table[1,n-2],-1]
A[4]//MatrixForm

and it’s inverse

Inverse[A[4]]//MatrixForm

Next we compute the second derivatives at the inner control points and we choose the second derivatives at the endpoints to be zero since the second derivative of the eigenfunction is zero:

s[n_?IntegerQ]:=Flatten[{0,Inverse[A[n]].Block[{cpn=cp[n]},6n^2(Drop[cpn,2]-2Drop[cpn,{1,n+1,n}]+Drop[cpn,-2])],0}]
s[4]

The coefficients of the cubic functions are

coef[n_?IntegerQ]:=Block[{cpn=cp[n],sn=s[n]},{Drop[cpn,-1],(Drop[cpn,1]-Drop[cpn,-1])*n-(2Drop[sn,-1]+Drop[sn,1])/(6n),Drop[sn,-1]/2,(Drop[sn,1]-Drop[sn,-1])*n/6}]
coef[4]

Block[{n=4,coefn=coef[n]},Plot[Sum[coefn[[i,Floor[n*t]+1]]*(t-Floor[n*t]/n)^(i-1),{i,1,4}],{t,0,1}]]

Upper estimate using the cubic spline as a test function

Let us compute the Rayleigh quotient:

Block[{n=4,coefn=coef[n]},ParallelSum[Integrate[Sum[(j-1)coefn[[j,i]]*t^(j-2),{j,2,4}]^2,{t,0,1/n}],{i,1,n}]/ParallelSum[Integrate[w[t+(i-1)/n]*Sum[coefn[[j,i]]*t^(j-1),{j,1,4}]^2,{t,0,1/n}],{i,1,n}]]//InputForm

This upper estimate is not guaranteed due to rounding errors. Hence, we use the interval arithmetic.

cp[n_?IntegerQ]:=Flatten[{0,Interval/@Table[pEigfN[t/n],{t,1,n-1}],0}]
cp[4]

s[4]

coef[4]

(*upperInt[n_?IntegerQ,m_?IntegerQ]:=Block[{coefn=coef[n]},Sum[Integrate[Sum[(j-1)coefn[[j,i]]*x^(j-2),{j,2,4}]^2,{x,0,1/n}],{i,1,n}]/Sum[Sum[Interval[{w[(i-1+(k-1)/m)/n],w[(i-1+k/m)/n]}]*Integrate[Sum[coefn[[j,i]]*x^(j-1),{j,1,4}]^2,{x,(k-1)/(m*n),k/(m*n)}],{k,1,m}],{i,1,n}]]*)
ParallelEvaluate[Off[General::munfl]];
intCS=Integrate[Sum[a[j]*t^(j-1),{j,1,4}]^2,t];
intDCS=Integrate[Sum[(j-1)a[j]*t^(j-2),{j,2,4}]^2,t];
upperInt[n_?IntegerQ,m_?IntegerQ]:=Block[{coefn=coef[n]},ParallelSum[(intDCS/.Table[a[j]->coefn[[j,i]],{j,2,4}])/.t->1/n,{i,1,n}]/ParallelSum[Sum[Interval[{w[(i-1+(k-1)/m)/n],w[(i-1+k/m)/n]}]*((#[[2]]-#[[1]])&@((intCS/.Table[a[j]->coefn[[j,i]],{j,1,4}])/.t->{(k-1)/(m*n),k/(m*n)})),{k,1,m}],{i,1,n}]]
upperInt[4,1000]

This interval must contain the exact value of the Rayleigh quotient for the spline function and, consequently, the right endpoint of this interval, i.e.,

Max[upperInt[4,1000]]
%-(Subscript[λ, e])

is a guaranteed upper estimate of the principal eigenvalue. If we want a more precise estimate, we can use more control points and/or a finer division in the integration:

upperInt[4,2000]
upperInt[8,500]
AbsoluteTiming[upperInt[20,10000]]
upperSpl=Max[%[[2]]];
upperSpl-(Subscript[λ, e])
InputForm[upperSpl]

Lower estimate using the cubic spline as a test function

Let us find a guaranteed lower estimate for((-v''(t))/(w(t)v(t))):

lowerInt[n_?IntegerQ,m_?IntegerQ]:=Block[{coefn=coef[n]},Min[Table[Min[Table[(Sum[-(j-1)(j-2)coefn[[j,i]]*t^(j-3),{j,3,4}]/Sum[coefn[[j,i]]*t^(j-1),{j,1,4}]/.t->Interval[{(k-1)/(m*n),k/(m*n)}])/Interval[{w[(i-1+(k-1)/m)/n],w[(i-1+k/m)/n]}],{k,1,m}]],{i,1,n}]]]
lowerInt[4,1000]

We get -∞ because the denominator is zero at the endpoints. Hence, zero is inside the first and the last corresponding interval and the ratio cannot be bounded from below:

Interval[{1,2}]/Interval[{-1,1}]

We can shift the test function by a small positive constant:

lowerInt[n_?IntegerQ,m_?IntegerQ,d_?NumberQ]:=Block[{coefn=coef[n]},Min[Table[Min[Table[(Sum[-(j-1)(j-2)coefn[[j,i]]*t^(j-3),{j,3,4}]/(Sum[coefn[[j,i]]*t^(j-1),{j,1,4}]+d)/.t->Interval[{(k-1)/(m*n),k/(m*n)}])/Interval[{w[(i-1+(k-1)/m)/n],w[(i-1+k/m)/n]}],{k,1,m}]],{i,1,n}]]]
lowerInt[4,1000,0.001]

We get a negative number near zero because the numerator is zero at the endpoints, too. But we choose the second derivative of the spline at the endpoints to be negative:

s[n_?IntegerQ,end_?NumberQ]:=Flatten[{end,Inverse[A[n]].Block[{cpn=cp[n]},6n^2(Drop[cpn,2]-2Drop[cpn,{1,n+1,n}]+Drop[cpn,-2])],end}]
s[4,-1]

coef[n_?IntegerQ,end_?NumberQ]:=Block[{cpn=cp[n],sn=s[n,end]},{Drop[cpn,-1],(Drop[cpn,1]-Drop[cpn,-1])*n-(2Drop[sn,-1]+Drop[sn,1])/(6n),Drop[sn,-1]/2,(Drop[sn,1]-Drop[sn,-1])*n/6}]
coef[4,-1]

lowerInt[n_?IntegerQ,m_?IntegerQ,d_?NumberQ,end_?NumberQ]:=Block[{coefn=coef[n,end]},Min[ParallelTable[Min[Table[(Sum[-(j-1)(j-2)coefn[[j,i]]*t^(j-3),{j,3,4}]/(Sum[coefn[[j,i]]*t^(j-1),{j,1,4}]+d)/.t->Interval[{(k-1)/(m*n),k/(m*n)}])/Interval[{w[(i-1+(k-1)/m)/n],w[(i-1+k/m)/n]}],{k,1,m}]],{i,1,n}]]]
lowerInt[4,1000,0.001,-1]

Manipulate[Block[{d=10^ld-10^-4,end=-(10^le-10^-4),coefn=coef[n,end]},Show[Plot[{(Subscript[λ, e]),Sum[-(i-1)(i-2)Min[coefn[[i,Floor[n*t]+1]]]*(t-Floor[n*t]/n)^(i-3),{i,3,4}]/(Sum[Min[coefn[[i,Floor[n*t]+1]]]*(t-Floor[n*t]/n)^(i-1),{i,1,4}]+d)/w[t]},{t,0,1},PlotRange->{-2,15}],Graphics[{Opacity[0.3],Green,EdgeForm[Thin],Flatten[Table[Table[Rectangle[{(i-1+(k-1)/m)/n,Max[Min[#],-100]},{(i-1+k/m)/n,Min[Max[#],100]}]&@((Sum[-(j-1)(j-2)coefn[[j,i]]*t^(j-3),{j,3,4}]/(Sum[coefn[[j,i]]*t^(j-1),{j,1,4}]+d)/.t->Interval[{(k-1)/(m*n),k/(m*n)}])/Interval[{w[(i-1+(k-1)/m)/n],w[(i-1+k/m)/n]}]),{k,1,m}],{i,1,n}]]}]]],{{n,4},3,6,1},{{m,3},1,10,1},{ld,-4,0},{le,-4,0}]

AbsoluteTiming[lowerInt[100,10000,10^-9,-1]]
lowerSpl=%[[2]];

Finally, the principal eigenvalue is guaranteed to lie in

Interval[{lowerSpl,upperSpl}]//InputForm

(Subscript[λ, e])//InputForm