謝資清

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

考慮如下半線(xiàn)性特征值問(wèn)題

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

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

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

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

1 求解半線(xiàn)性特征值問(wèn)題的搜索延拓法

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

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

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

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

第二步 在較小的子空間Sk中逼近模型問(wèn)題以獲得特征對(duì)的多個(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）的非線(xiàn)性代數(shù)方程組

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

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

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

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

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

類(lèi)似于第二步的討論，記

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

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

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

由于第三步和第四步都涉及用數(shù)值延拓法求解非線(xiàn)性代數(shù)方程組，以方程組（5）為例介紹一類(lèi)典型的牛頓型數(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.考慮子問(wèn)題

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

1）依次對(duì)i=1，2，…，m -1，以zi，0=zi-1為初值，用如下牛頓迭代求解子問(wèn)題（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｜＜ε，得問(wèn)題（5）的近似解z*：=zk+1.

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

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

其中PN為定義在I =［-1，1］上且次數(shù)不超過(guò)N的多項(xiàng)式的全體.令為對(duì)應(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)的拉格朗日插值算子，滿(mǎn)足：

注意，（6）式中的INf（uN）是對(duì)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，得到如下矩陣形式的非線(xiàn)性代數(shù)方程組

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

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

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

其中

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

2 一類(lèi)立方非線(xiàn)性特征值問(wèn)題的計(jì)算

取f（u）=u3，Ω=（0，π）2，考慮半線(xiàn)性特征值問(wèn)題

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

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

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

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

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

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

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

考慮積分

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

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

ii）若有且僅有2 對(duì)（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，，可以得到解

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

由此得到2 個(gè)（9）式的特征對(duì)數(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對(duì)應(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è)特征對(duì)的初值

由此得到8 個(gè)（9）式的特征對(duì)數(shù)值解.對(duì)（i，j）=（1，2），（1，3）和（1，5），對(duì)應(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）式的特征對(duì)數(shù)值解，其中6 種特征函數(shù)的數(shù)值解圖像繪于圖5 和6，相應(yīng)的特征值和初值信息列于表3.

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

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

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

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

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

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