引 言
自然科学和工程技术领域内经常会遇到随着参数变化而发生状态突变的非线性问题，例如描述化学相变过程的反应扩散方程、弹性结构的屈曲和颤振问题以及黏性流体的稳定性及流动分离问题等[1]．分岔分析是研究这类非线性问题及其突变现象的强有力手段，因此这类问题有时又称作非线性分岔问题[2]．分岔分析早已从临界状态邻域内的局部分析发展到了非线性分岔问题的大范围分析，即要求在参数大范围变化内计算分岔点、追踪多解及分析解的稳定性．
非线性分岔分析是建立在含参数的非线性代数方程组和常微分方程组基础之上的，而非线性分岔问题往往为时空坐标下的偏微分方程．这意味着非线性分岔问题的数值求解强烈依赖于空间离散方法的精度．非线性分岔分析的数值计算方法一般可分为两类：一是连续追踪平衡解的伪弧长法[3]；一是直接求解临界参数的扩展方程法[4]，两种方法各有优缺点．伪弧长法是有非线性限元增量法的一般化，需要较高的数值技巧来定位分岔点和分岔方向，是一种间接方法；而扩展方程法将原方程的奇异点化为一个新方程的正则解，可求解极值型、音叉型和Hopf型等多种分岔点，其缺点为新方程的维数会成倍增加．这两种方法都需要离散方程组具有良好性质的系数矩阵及Jacobi矩阵．理论研究则表明，小波Galerkin法得到的系数矩阵条件数是有界的，不会造成矩阵的病态[5]．利用Haar小波对非线性反应扩散方程进行分岔分析也表明，小波方法能在较小的分解尺度下捕捉到方程的极值型和音叉型分岔点[6]．同时笔者在研究过程中观察到，基于Coiflet的小波Galerkin法得到的离散方程组可显式地推导Jacobi矩阵甚至更高阶的Hesse矩阵，而不需要采用数值公式来计算这些矩阵[7]．基于Coiflet的小波Galerkin法及其在强非线性力学问题求解方面的研究进一步表明，小波Galerkin法对问题的非线性强弱不敏感，具有良好的计算精度和收敛阶数，尤其适用于梁板壳大挠度和非线性振动等结构非线性分析[8-12]．这些研究均表明，小波Galerkin法不仅适用于求解非线性分岔问题，而且有可能在结构屈曲和动力失稳等复杂分岔问题求解方面具有更好的计算效率和精度．
本文通过一个典型的Bratu分岔问题来研究小波Galerkin法求解非线性分岔问题的具体过程．Bratu方程可用于描述固体推进剂放热反应和宇宙膨胀过程等物理现象，是有限单元法[13]、谱方法[14]和无网格法[15-16]等方法发展过程中的典型算例之一．不同的是，有限单元法多使用增量法来计算临界参数，而谱方法和无网格法则采用延拓法．此外，一些从事计算数学的学者在研究非线性方程组极值解的问题时，也以热爆炸问题为研究对象提出了一种直接求极值型分岔点的方法[17-18]．这些工作大多以限差分法(FDM)作为空间离散方法．由于Bratu问题具有十分典型的分岔行为以及重要的应用背景，经常用作检验非线性分岔问题新的计算方法[19]．因此本文选择该问题进行分岔分析，首先给出了一维和二维问题的小波离散格式；然后推导了单参数和双参数Bratu问题的伪弧长格式和扩展方程；最后讨论了小波Galerkin法求解的数值结果和计算效率．
1 小波离散格式
1.1 Bratu非线性分岔问题
Bratu方程是描述固体推进剂放热反应热平衡问题的基本方程，其无量纲的形式为
(1)
其中u为系统温升， λ和ε为无量纲参数， λ又称为Frank-Kamenetski参数， 这两个参数分别表示系统的反应速率和活化能．为了便于说明， 边界条件采用Dirichlet齐次边界．该问题为经典的稳态热平衡双参数分岔问题．
1.2 分岔问题的小波离散格式
小波Galerkin法采用基于Coiflet的修正基函数作为权函数与试函数．对于有限区间上可积函数空间L2[0,1]内任意函数f(x)，有逼近公式如下[20]：
(2)
其中Pjf(x)为f(x)在小波子空间上的投影；基函数φj,k(x)的表达式为
(3)
修正系数为
(4)
引入β0和β1是为说明边界条件的施加过程，为了满足x=0处的齐次边界，仅仅需要令β0,i=0,而其他元素保持不变即可．根据四点Markov数值微分公式推导了ζ0和ζ1[20]:
(5)
根据小波理论，对于定义在有限区域上的二维函数f(x,y)∈L2[0,1]2,类似地有
(6)
对于一维Bratu分岔问题，对方程中的各项函数采用小波逼近公式进行离散，代入可得
(7)
选择权函数hj,k(x)，表示满足边界条件的修正基函数φj,k(x)．应用Galerkin法，可得
A1U1+λB1 f(U1,ε)+λD1≈0，
(8)
其中系数矩阵和未知向量分别为
D1={b1,l0+b1,l2j},
U1={u1,k=u(xk),xk=k/2j},
f(U1,ε)={f(u1,k,ε)=eu(xk)/(1+εu(xk))}．
类似地，对于二维Bratu分岔问题，离散可得
(9)
应用Galerkin法，可得
A2U2+λB2f(U2,ε)+λD2≈0,
(10)
其中系数矩阵为
A2={a2,rs,lk=crkdsl+drkcsl}，
dsl=
hj,l(y)hj,s(y)dy，
B2={b2,rs,lk=erkesl},
erk=φj,k(x)hj,i(x)dx,
k,l=1,2,…,2j-1，元素为尺度函数与其导乘积的积分，D2={b2,rs,l0+b2,rs,l2m+b2,rs,0k+b2,rs,2jk}，未知向量为f(U2,ε)={f(u2,kl,ε)=eu(xk,yl)/(1+εu(xk,yl))}和U2={u2,kl=u(xk,yl),xk=k/2j,yl=l/2j,kl=(k-1)(2j-1)+l}．本文选用小波函数消失矩数目N=6，尺度函数一阶矩M1=7的Coiflet小波，上述积分系数的计算参见文献[20]．
2 分岔计算方法
将一维和二维Bratu分岔问题写作统一形式：
AiUi+λBif(Ui,ε)+λDi≈0，
(11)
下标1和2表示一维和二维，参数λ是该问题的主要控制参数，固定参数ε利用伪弧长法可在参数λ大范围内追踪其平衡解以及利用扩展方程可以直接求解简单的极值型分岔点．对于双参数问题，利用伪弧长法求解扩充方程则可以追踪双参数平面内的平衡解，利用高阶扩充方程则可以求解尖点型分岔点的值．
2.1 追踪平衡解的伪弧长格式
伪弧长法本质上是一种延拓法，其基本思想是将控制参数λ看作与状态变量一样的变量，并且引入新的弧长参数s．那么对于某个固定参数ε0，式(11)中的分岔问题改写为[3]
F(Ui,λ,ε0)=AiUi+λBif(Ui,ε0)+λDi=0．
(12)
选择适当的约束方程，且令y=[Ui,λ],式(12)可进一步写为
(13)
其中τ是解曲线的单位切向量，定义式为τ=v(y)/‖v(y)‖,切向量为v(y)=[J1,J2,…,Jn+1]T．切向量的元素则为
(14)
其中
表示划去该列．式(13)可化为常微分方程的Cauchy问题:
(15)
为了追踪解曲线，采用预报-修正格式如下[3]:
yk=yk-1+τ(yk-1)·(sk-sk-1)，
(16)
(17)
其中式(16)为Euler预报格式，式(17)为Newton修正格式，迭代序列
收敛于解曲线y(s)上的某一点，DF(yk)则为预报点yk处的Fréchet导数．由式(12)可推得
(18)
一维问题和二维问题中的向量g(Ui,ε0)分别为g(U1,ε0)={g(u1,k,ε0)}和g(U2,ε0)={g(u2,kl,ε0)}, 元素则为g(ui,ε0)=f(ui,ε0)/(1+ε0ui)2,其中ui=u1,k和ui=u2,kl．
2.2 求解临界参数的扩展方程
对于某个固定参数ε0，分岔问题仅存在简单的极值型分岔点，式(12)的扩展方程为[4]
(19)
其中v是系数矩阵的特征向量，l是一个任意向量．为了追踪参数平面内分岔问题的平衡解和稳定边界，可以利用伪弧长法来求解式(19)中的扩充系统．令变量为z=(Ui,v,λ,ε),单位切向量为τ=v(z)/‖v(z)‖,v(z)=[J1,J2,…,Jn+1]T．类似可得，切向量元素为
(20)
以及Fréchet导数表示为
DG(z)=
(21)
其中一维和二维Bratu分岔问题的向量r(Ui,ε)，h(Ui,ε)和w(Ui,ε)分别为r(U1,ε)={r(u1,k,ε)}, h(U1,ε)={h(u1,k,ε)},w(U1,ε)={w(u1,k,ε)},以及r(U2,ε)={r(u2,kl,ε)}，h(U2,ε)={h(u2,kl,ε)}和w(U2,ε)={w(u2,kl,ε)},元素则为
和
对于式(11)中的双参数分岔问题，还存在尖点型分岔点，其扩展方程为[4]
(22)
其中v和ξ为特征向量，l是一个任意向量．这里取l=DλF(Ui,λ,ε)=Bif(Ui,ε)+Di．{Ui, v, λ}和{Ui, v，ξ, λ, ε} 是式(19)和(22)的待求变量向量．二阶Fréchet导数推得为
DxxF(Ui,λ,ε)vv=λBidiag[h(Ui,ε)]v2．
(23)
对于式(19)和(22)中的扩展方程，尽管其维数相对于小波离散方程的维数扩大了一倍，利用换元变换可对其进行降阶求解．本文主要关注小波Galerkin法在求解非线性分岔问题的求解过程和计算精度，因此仍然采用Newton迭代法直接求解上述扩展方程，其中取零向量作为特征向量的初始值．
3 结果与讨论
一维和二维Bratu分岔问题分别描述无限大平板形状和无限长方形杆形状的放热反应系统热爆炸问题，其包含了两个控制参数λ和ε．
3.1 单参数 Bratu分岔问题
对于ε0=0的一维Bratu问题，当参数λ>λ0, λ=λ0和λ<λ0时，无解、有1个解和2个解，其中临界参数λ0满足下列关系式：
(24)
(a) 一维问题 (b) 二维问题
(a) The 1D case (b) The 2D case
图1 单参数Bratu问题的解曲线
Fig. 1 Solution curves of 1-parameter Bratu problems
(a) 一维问题 (b) 二维问题
(a) The 1D case (b) The 2D case
图2 临界参数绝对误差与离散方程维数N的关系
Fig. 2 Relationships between absolute errors of the critical parameter and number of discretized equations N
一维问题临界参数精确值为λ0=3.513 830 719 1，对于ε0=0的二维Bratu问题，尚未找到临界参数的解析表达式[19]，数值结果为λ0=6.808 124 43．单参数Bratu问题的解曲线由Boyd求得，并得到了不同方法的验证[15]．利用式(12)的伪弧长格式，本文也获得了与解析曲线和Boyd曲线完全一致的分岔曲线，如图1所示，其中小波Galerkin法的分解尺度取j=4．临界参数λ0为放热反应系统热爆炸的判据，由式(19)扩展方程中可直接求出．图2给出了小波Galerkin法和间断Galerkin法[19]求得的临界参数绝对误差与离散方程维数N的关系，可见小波Galerkin法不仅具有更高的计算精度，而且具有更快的收敛速度．
(a) 折迭线 (b) u(0.5)与参数ε的关系[17]
(a) Fold lines (b) Relationships between u(0.5) and ε[17]
图3 一维双参数Bratu问题
Fig. 3 The 1D Bratu problem with 2 parameters
(a) 折迭线 (b) u(0.5,0.5)与参数ε的关系[14]
(a) Fold lines (b) Relationships between u(0.5,0.5) and ε[14]
图4 二维双参数Bratu问题
Fig. 4 The 2D Bratu problem with 2 parameters
3.2 双参数Bratu分岔问题
对于双参数问题，一般通过求解参数平面上的折迭线来描述其分岔行为．折迭线的物理意义为非线性分岔问题的稳定边界，其本质上则是分岔点的集合．利用式(19)的伪弧长格式可以连续追踪双参数问题的稳定边界，图3和图4分别给出了一维和二维双参数Bratu问题的折迭线以及中心温升与参数ε的关系曲线．小波Galerkin法的结果与其他方法获得的结果是一致的，显示了当ε=0，0<ε<εcr和ε≥εcr时，双参数Bratu问题分别存在1个和2个以及没有极值型分岔点[17]．当0<ε<εcr时，两个极值型分岔点分别表示升温过程中的点火条件和降温过程中的熄火条件．同时可以看到，随着参数ε靠近临界值εcr，两个分岔点的距离越来越小直至重合．临界值εcr为放热反应系统的转变条件，即系统的热爆炸特性消失，系统处于连续的变化之中．
表1 尖点临界值
Table 1 Critical value of the cusp
method1D Bratu problemεcrλcru(0.5)2D Bratu problemεcrλcru(0.5,0.5)FDM[4]0.245 785.229 594.896 550.241 77910.239 2565.979 54WGM0.245 780 763 25.229 503 360 14.898 866 621 40.241 848 405 110.241 370 044 55.988 446 899 8
临界参数(λcr, εcr)一般称之为双参数分岔问题的尖点，可由式(22)的扩展方程求解．表1给出了小波Galerkin法求得的尖点值，与文献中的结果一致[4]．除了稳定边界之外，利用式(19)的伪弧长格式还可以获得尖点突变的平衡曲面，如图5所示．值得指出的是，可以通过对式(14)和(20)中的矩阵进行LU分解，且乘以矩阵元素的最大值来提高伪弧长格式的稳健性．
(a) 一维问题 (b) 二维问题
(a) The 1D case (b) The 2D case
图5 双参数Bratu问题的解曲面
Fig. 5 Solution surfaces of the Bratu problem with 2 parameters
4 结 语
本文利用基于Coiflet的小波Galerkin法求解了典型的Bratu非线性分岔问题，数值结果表明，小波Galerkin法在求解临界参数时具有更高的计算精度．求解过程也进一步地说明，基于小波Galerkin法的离散格式可直接推导伪弧长法和扩展方程等用于非线性分岔计算的求解格式，从而可直接求解极值型和尖点型等不同类型分岔点以及多参数的尖点突变曲面．由于小波Galerkin法在结构大挠度和非线性振动问题求解方面的良好表现，那么对非线性分岔问题求解中系数矩阵性质以及扩展方程组直接求解的进一步研究，则有望为非线性屈曲和动力失稳等一般结构非线性稳定性问题提供新的求解思路．
[1] 陈予恕, 唐云. 非线性动力学中的现代分析方法[M]. 北京: 科学出版社, 1992.(CHEN Yushu, TANG Yun. Modern Analytical Methods in Nonlinear Dynamics[M]. Beijing: Science Press, 1992.(in Chinese))
[2] SEYDEL R. Practical Bifurcation and Stability Analysis[M]. New York: Springer, 2009.
[3] 季海波, 武际可. 分叉问题及其数值方法[J]. 力学进展, 1993, 23(4): 493-502.(JI Haibo, WU Jike. Bifurcation and its numerical methods[J]. Advances in Mechancis, 1993, 23(4): 493-502.(in Chinese))
[4] ROOSE D, PIESSENS R. Numerical computation of non-simple turning points and cusps[J]. Numerische Mathematik, 1985, 46(2): 189-211.
[5] 陈仲英, 巫斌. 小波分析[M]. 北京: 科学出版社, 2007.(CHEN Zhongying, WU Bin. Wavelet Analysis[M]. Beijing: Science Press, 2007.(in Chinese))
[6] KRISHNAN J, RUNBORG O, KEVREKIDIS I G. Bifurcation analysis of nonlinear reaction-diffusion problems using wavelet: based reduction techniques[J]. Computers & Chemical Engineering, 2004, 28(4): 557-574.
[7] 张磊. 高精度小波数值方法及其在结构非线性分析中的应用[D]. 博士学位论文. 兰州: 兰州大学, 2016.(ZHANG Lei. High-precision wavelet numerical methods and applications to nonlinear structural analysis[D]. PhD Thesis. Lanzhou: Lanzhou University, 2016.(in Chinese))
[8] 刘小靖, 王记增, 周又和. 一种适用于强非线性结构力学问题数值求解的修正小波伽辽金方法[J]. 固体力学学报, 2011, 32(3): 249-257.(LIU Xiaojing, WANG Jizeng, ZHOU Youhe. A modified wavelet Galerkin method for computations in structural mechancis with strong nonlinearity[J]. Chinese Journal of Solid Mechanics, 2011, 32(3): 249-257.(in Chinese))
[9] LIU X J, ZHOU Y H, WANG X M, et al. A wavelet method for solving a class of nonlinear boundary value problems[J]. Communications in Nonlinear Science and Numerical Simulation, 2013, 18(8): 1939-1948.
[10] WANG X M, LIU X J, WANG J Z, et al. A wavelet method for bending of circular plate with large defelction[J]. Acta Mech Solida Sinica, 2015, 28(1): 83-90.
[11] ZHANG L, WANG J Z, ZHOU Y H. Wavelet solution for large deflection bending problems of thin rectangular plates[J]. Archive of Applied Mechanics, 2015, 85(3): 355-365.
[12] MA X L, WU B, ZHANG J H, et al. A new numerical scheme with wavelet-Galerkin followed by spectral deferred correction for solving string vibration problems[J]. Mechanism and Machine Theory, 2019, 142(12): 103623.
[13] ANDERSON C A, ZIENKIEWICZ O C. Spontanenous ignition: finite element solutions for steady state and transient conditions[J]. Journal of Heat Transfer, 1974, 96(3): 398-404.
[14] KAPANIA R K. A pseudo-spectral solution of 2-parameter Bratu’s equation[J]. Computational Mechanics, 1990, 6(1): 55-63.
[15] KARKOWSKI J. Numerical experiments with the Bratu equation in one, two and three dimensions[J]. Computational and Applied Mathematics, 2013, 32(1): 231-244.
[16] 洪文强, 徐绩青, 许锡宾, 等. 求解Bratu型方程的径向基函数逼近法[J]. 应用数学和力学, 2016, 37(6): 617-625.(HONG Wenqiang, XU Jiqing, XU Xibin, et al. The radial basis function approximation method for solving Bratu-type equations[J]. Applied Mathematics and Mechanics, 2016, 37(6): 617-625.(in Chinese))
[17] BODDINGTON T, FENG C G, GRAY P. Thermal explosion and the theory of its initiation by steady intense light[J]. Proceedings of the Royal Society A: Mathematical, 1983, 390(1799): 265-281.
[18] GREENWAY P, SPENCE A. Numerical calculation of critical points for a slab with partial insulation[J]. Combust Flame, 1985, 62(2): 141-156.
[19] CLIFFE K A, HALL E, HOUSTON P, et al. Adaptivity and a posteriori error control for bifurcation problems Ⅰ: the Bratu problem[J]. Communications in Computational Physics, 2010, 8(4): 845-865.
[20] 王记增. 正交小波统一理论与方法及其在压电智能结构等力学研究中的应用[D]. 博士学位论文. 兰州: 兰州大学, 2001.(WANG Jizeng. Generalized theory and arithmetic of orthogonal wavelets and applications to researches of mechanics including piezoeiectric smart structures[D]. PhD Thesis. Lanzhou: Lanzhou University, 2001.(in Chinese))