求解半線性特征值問題的搜索延拓法

謝資清

（湖南師范大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院 計(jì)算與隨機(jī)數(shù)學(xué)教育部重點(diǎn)實(shí)驗(yàn)室，湖南 長沙410081）

考慮如下半線性特征值問題

其中，L是一致橢圓算子，Ω?Rd是一個(gè)光滑有界區(qū)域，‖·‖為L2（Ω）范數(shù)，f是一個(gè)光滑非線性函數(shù)且滿足f（0）=0.

模型問題（1）廣泛出現(xiàn)在現(xiàn)代科學(xué)的各個(gè)領(lǐng)域，例如：激光傳輸［1］、電子結(jié)構(gòu)計(jì)算［2］、玻色-愛因斯坦凝聚（BEC）［3］等.近年來，人們對模型問題（1）的一個(gè)主要研究興趣是設(shè)計(jì)計(jì)算最小能量特征態(tài)（基態(tài)）的數(shù)值方法和建立相關(guān)理論［4-8］.隨著科學(xué)技術(shù)的發(fā)展，具有更高能量的特征態(tài)（激發(fā)態(tài)）越來越受到科學(xué)家們關(guān)注.Bao 等［9］對一些用于刻畫BEC的半線性薛定諤特征值問題，使用Hermite 函數(shù)構(gòu)造好的初值，從數(shù)值上模擬了問題的基態(tài)和幾種激發(fā)態(tài).Chang 等［10］設(shè)計(jì)了計(jì)算BEC 能級的自適應(yīng)同倫算法，能有效計(jì)算BEC 的不同能級及相應(yīng)的定態(tài)解.此外，Zhou 等［11-12］將他們提出的計(jì)算非線性方程多解的局部極小極大方法推廣到一類非線性特征值問題，計(jì)算了多個(gè)特征對，并建立了相關(guān)理論.

由于SEM簡單、高效的優(yōu)勢，本文對原有SEM進(jìn)行改進(jìn)，使其適用于求解半線性特征值問題（1）.為此，首先在由線性算子L 的少數(shù)幾個(gè)特征函數(shù)張成的子空間中逼近模型問題，得到小規(guī)模非線性代數(shù)特征值問題，求解該問題獲得模型問題特征對的多個(gè)初值；然后，采用L 的更多的特征函數(shù)，進(jìn)一步逼近模型問題的特征對以得到更好的初值；再用合適的數(shù)值方法離散模型問題，得到相對大規(guī)模的非線性代數(shù)方程組（非線性代數(shù)特征值問題）；最后用數(shù)值延拓法求解此非線性代數(shù)方程組，得到每個(gè)初值對應(yīng)的特征對.為提高計(jì)算效率，本為將提出一種插值系數(shù)Legendre -Galerkin 譜方法用于離散模型問題（1）.

本文結(jié)構(gòu)安排如下：第1 節(jié)給出求解模型問題（1）的SEM 的算法步驟和二維情形的插值系數(shù)Legendre-Galerkin譜離散格式；第2 節(jié)給出求解模型問題（1）的SEM 算法步驟和二維情形的插值系數(shù)Legendre-Galerkin 譜離散格式；第3 節(jié)應(yīng)用改進(jìn)的SEM 求解一類立方非線性特征值問題，給出計(jì)算實(shí)現(xiàn)的細(xì)節(jié)，并展示和分析數(shù)值結(jié)果.

1 求解半線性特征值問題的搜索延拓法

不失一般性，以下考慮L = -Δ的情形，則模型問題（1）的弱形式為：求（u，μ）∈（Ω）×R使其中（·，·）表示L2（Ω）內(nèi)積，A（u，v）：=（▽u，▽v）.

1.1 SEM算法步驟基于原始SEM的思想，求解半線性特征值問題（1）的SEM 依次在3 層子空間Sk?Sn?SN上逼近模型問題，分為以下4 個(gè)步驟：

第一步 計(jì)算線性算子-Δ的特征對眾所周知，線性特征值問題

有可數(shù)無限多特征值：

第二步 在較小的子空間Sk中逼近模型問題以獲得特征對的多個(gè)初值.

定義Sk=span｛φnj：1≤j≤k｝，其中nj為互不相同的正整數(shù).求（uk，μk）∈Sk×R使

注意到

在（3）式中，令

并取v=φni，i=1，2，…，k，得到關(guān)于k+1 個(gè)未知數(shù)（a1，…，ak，μk）的非線性代數(shù)方程組

其中a=（a1，a2，…，ak），

因?yàn)閗較小，可通過解析方法或牛頓法等求得非線性方程組（4）的多解.特別地，當(dāng)f 為多項(xiàng)式函數(shù)時(shí)，可利用數(shù)值代數(shù)幾何方法［19］得到（4）式的所有解.顯然，每組解）對應(yīng)模型問題（1）的特征對的一個(gè)初值

第三步 在較大子空間Sn中進(jìn)一步逼近模型問題的特征對以得到更好的初值.

為進(jìn)一步逼近模型問題的特征對，取一個(gè)適當(dāng)大的子空間Sn使Sk?Sn.設(shè)

求（un，μn）∈Sn×R使

類似于第二步的討論，記

得到關(guān)于z=（z1，…，zn，μn）T∈Rn+1的非線性代數(shù)方程組

第四步 在大子空間SN中用合適數(shù)值方法離散模型問題并求解相應(yīng)的非線性代數(shù)方程組得到最終數(shù)值解.

在SN中使用有限元、插值系數(shù)有限元或譜方法等離散模型問題（1），并用第二步或第三步得到的每個(gè)逼近解（u0，μ0）為初值，采用數(shù)值延拓法求解相應(yīng)的非線性代數(shù)方程組，得到最終的數(shù)值解（uN，μN(yùn)）∈SN×R.

由于第三步和第四步都涉及用數(shù)值延拓法求解非線性代數(shù)方程組，以方程組（5）為例介紹一類典型的牛頓型數(shù)值延拓法［15］.記J（z）為F（z）的雅可比矩陣.令G（z）= J（z0）（z - z0）.構(gòu)造凸同倫函數(shù)

顯然，H（z0，0）=0.選取適當(dāng)大的正整數(shù)m 和剖分0 =t0＜t1＜t2＜…＜tm=1.考慮子問題

給定精度0 ＜ε?ε1?1.求解非線性代數(shù)方程組（5）的牛頓型數(shù)值延拓法的步驟如下：

1）依次對i=1，2，…，m -1，以zi，0=zi-1為初值，用如下牛頓迭代求解子問題（Pi）

zi，k+1= zi，k-（Ji，k）-1H（zi，k，ti）， k = 0，1，…，其中Ji，k=Hz′（zi，k，ti）=tiJ（zi，k）+（1 -ti）J（z0），直 到 條 件｜zi，k+1- zi，k｜ ＜ε1成 立，得 近 似 解

zi：=zi，k+1.

當(dāng)｜zk+1-zk｜＜ε，得問題（5）的近似解z*：=zk+1.

1.2 插值系數(shù)Legendre-Galerkin譜離散由于譜方法具有高精度的優(yōu)點(diǎn)，受插值系數(shù)有限元方法［1，16］的啟發(fā)，本文提出一種插值系數(shù)Legendre -Galerkin譜方法用于離散半線性特征值問題（1）.下面給出二維情形的離散格式，三維和更高維情形可類似得到.

為簡單起見，假設(shè)Ω =（-1，1）2.給定適當(dāng)大的正整數(shù)N，考慮多項(xiàng)式逼近空間

其中PN為定義在I =［-1，1］上且次數(shù)不超過N的多項(xiàng)式的全體.令為對應(yīng)于Legendre -Gauss-Lobatto（LGL）點(diǎn)的拉格朗日多項(xiàng)式基（節(jié)點(diǎn)基），則二維張量形式的LGL 點(diǎn)和相應(yīng)的節(jié)點(diǎn)基分別為

提出如下插值系數(shù)Legendre - Galerkin 譜方法：求（uN，μN(yùn)）∈SN×R使

其中IN：C（ˉΩ）→PN?PN是基于二維LGL 點(diǎn)的拉格朗日插值算子，滿足：

注意，（6）式中的INf（uN）是對f（uN）作整體插值.這是（6）式與標(biāo)準(zhǔn)Legendre -Galerkin 譜方法的重要區(qū)別.

注意到f（0）=0，有

因?yàn)閡N∈SN，所以

代入（6）式，并取vN=hlm，l，m =1，2，…，q =N -1，得到如下矩陣形式的非線性代數(shù)方程組

其中“∶”表示矩陣的Frobenius內(nèi)積，定義為

記K=B?A+A?B，M=B?B，其中“?”表示Kronecker乘積算子，定義為則方程組（7）可重寫為

容易得出F的雅可比矩陣J為

其中

由此可見，插值系數(shù)法的運(yùn)用，使得雅可比矩陣的表達(dá)式（8）比標(biāo)準(zhǔn)Galerkin譜方法簡單.事實(shí)上，由于K、M是常數(shù)矩陣，每步牛頓迭代中，更新雅可比矩陣的主要工作量僅是計(jì)算對角矩陣Df（ˉu）.這是插值系數(shù)法的主要優(yōu)點(diǎn).

2 一類立方非線性特征值問題的計(jì)算

取f（u）=u3，Ω=（0，π）2，考慮半線性特征值問題

下面應(yīng)用第2 節(jié)描述的SEM 求（9）式的多特征對（u，μ）.

首先，在SEM的第一步，直接計(jì)算得線性特征問題

的特征值和相應(yīng)的規(guī)范正交特征函數(shù)分別為

在第二步，考慮Sk是由-Δ的對應(yīng)同一特征值λ的特征函數(shù)張成的子空間.設(shè)λ 是k 重特征值，則存在k個(gè)不同的正整數(shù)對（ip，jp）∈N+×N+使得

其中a=（a1，a2，…，ak），

定理2.1非線性代數(shù)方程組（10）至少有3k-1 個(gè)解.

證明設(shè)（ip，jp）∈N+×N+滿足

考慮積分

直接計(jì)算可知［15，20］：

i）若（i1，j1）=（i2，j2）=（i3，j3）=（i4，j4），Iλ=9/（4π2）；

ii）若有且僅有2 對（ip，jp）分別相等，Iλ=1/π2；

iii）其他情形，Iλ=0.于是

將（11）式代入（10）式得

設(shè)r是向量a=（a1，a2，…，ak）的非零分量個(gè)數(shù)，顯然1≤r≤k.不妨假設(shè)，i =1，2，…，r，，可以得到解

下面分別對k =1，2，3，4 情形，考慮Sk由-Δ的k重特征值對應(yīng)的特征函數(shù)張成，給出相應(yīng)的數(shù)值結(jié)果.在數(shù)值計(jì)算中，算法第三步的Sn通常由5～20 個(gè)基函數(shù)張成.在第四步，取N =32 進(jìn)行離散.數(shù)值延拓法子問題個(gè)數(shù)取為10、牛頓迭代的精度控制數(shù)取為以下特征函數(shù)的圖像均在一個(gè)64 ×64 網(wǎng)格上呈現(xiàn).

由此得到2 個(gè)（9）式的特征對數(shù)值解.圖1 給出了i=1，2，3 情形得到的（9）的特征函數(shù)uii的圖像，相應(yīng)的特征值和初值信息列于表1.

圖1 特征函數(shù)u11（左）、u22（中）及u33（右）Fig.1 Eigenfunctions u11（left），u22（middle）and u33（right）

表1 圖1 所示特征函數(shù)u對應(yīng)的特征值μ及初值（u0，μ0）Tab.1 Eigenvalues μ corresponding to eigenfunctions u shown in Fig.1 and their initial guesses（u0，μ0）

2）重特征情形.取

根據(jù)定理3.1，算法第二步得到32-1 =8 個(gè)特征對的初值

由此得到8 個(gè)（9）式的特征對數(shù)值解.對（i，j）=（1，2），（1，3）和（1，5），對應(yīng)的幾種特征函數(shù)的數(shù)值解圖像繪于圖2 ～4，相應(yīng)的特征值和初值信息列于表2.

圖2 特征函數(shù)u12（左）、u12+21（中）及u12-21（右）Fig.2 Eigenfunctions u12（left），u12+21（middle）and u12-21（right）

圖3 特征函數(shù)u13（左）、u13+31（中）、及u13-31（右）Fig.3 Eigenfunctions u13（left），u13+31（middle）and u13-31（right）

圖4 特征函數(shù)u15（左）、u15+51（中）及u15-51（右）Fig.4 Eigenfunctions u15（left），u15+51（middle）and u15-51（right）

由此得到26 個(gè)（9）式的特征對數(shù)值解，其中6 種特征函數(shù)的數(shù)值解圖像繪于圖5 和6，相應(yīng)的特征值和初值信息列于表3.

表2 圖2 ～4 所示特征函數(shù)u對應(yīng)的特征值μ及初值（u0，μ0）Tab.2 Eigenvalues μ corresponding to eigenfunctions u shown in Figs.2-4 and their initial guesses（u0，μ0）

圖5 特征函數(shù)u17（左）、u55（中）及u17+71（右）Fig.5 Eigenfunctions u17（left），u55（middle）and u17+71（right）

圖6 特征函數(shù)u17-71（左）、u17+55+71（中）及u17-55+71（右）Fig.6 Eigenfunctions u17-71（left），u17+55+71（middle）and u17-55+71（right）

一般地，根據(jù)定理3.1，對于-Δ的任何k重特征值λ，可在算法第二步構(gòu)造3k-1 個(gè)（9）式的特征對的初值.因此，理論上可以繼續(xù)使用-Δ的更高重?cái)?shù)的特征值和其對應(yīng)的特征函數(shù)構(gòu)造初值，并計(jì)算相應(yīng)的特征對的數(shù)值解.然而，由于這時(shí)解的高度不穩(wěn)定性，需要使用更大的逼近空間SN才能準(zhǔn)確計(jì)算相應(yīng)的數(shù)值解.例如最小的5 重特征值和最小的6 重特征值分別為

由于它們的值較大，本文暫未計(jì)算它們對應(yīng)的解.

表3 圖5 ～6 所示特征函數(shù)u對應(yīng)的特征值μ及初值（u0，μ0）Tab.3 Eigenvalues μ corresponding to eigenfunctions u shown in Figs.5-6 and their initial guesses（u0，μ0）

圖7 特征函數(shù)u18（左）、u47（中）及u18+47（右）Fig.7 Eigenfunctions u18（left），u47（middle）and u18+47（right）

圖8 特征函數(shù)u47+74（左）、u18+47+74（中）及u18+47+81（右）Fig.8 Eigenfunctions u47+74（left），u18+47+74（middle）and u18+47+81（right）

圖9 特征函數(shù)u18+47+74+81（左）、u18-47+74+81（中）及u18-47-74+81（右）Fig.9 Eigenfunctions u18+47+74+81（left），u18-47+74+81（middle）and u18-47-74+81（right）

表4 圖7 ～9 所示特征函數(shù)u對應(yīng)的特征值μ及初值（u0，μ0）Tab.4 Eigenvalues μ corresponding to eigenfunctions u shown in Figs.7-9 and their initial guesses（u0，μ0）

基于上述數(shù)值結(jié)果，觀測到如下現(xiàn)象：

若λ是-Δ的k 重特征值，相應(yīng)地，半線性特征值問題（9）至少存在3k-1 個(gè)不同的特征對（若不區(qū)分同一特征值對應(yīng)的正負(fù)特征函數(shù)，則問題（9）至少存在對應(yīng)于λ的（3k-1）/2 個(gè)不同的特征對）.也就是說，半線性特征值問題（9）的特征對能按-Δ特征值大小排序，且密度遠(yuǎn)大于-Δ的特征對.它們的個(gè)數(shù)和分布類似樹形結(jié)構(gòu)，且在-Δ的多重特征值的地方派生出更多樹枝.

值得一提的是，這一關(guān)于半線性特征值問題的特征對分布的數(shù)值觀測，與Chen 等［15］提出并被Zhang等［20］證明的關(guān)于立方非線性橢圓方程多解分布的猜想非常類似.對上述數(shù)值觀測的嚴(yán)格數(shù)學(xué)分析或證明還有待進(jìn)一步研究.