第06章MB-方程的数值解

131第六章MB方程的数值解前面我们建立了光学Bloch方程和MB方程，并求得了在特定条件下的解析解，应用这些解析解的结果，很好地解释了超短激光脉冲与物质的近共振相互作用的许多物理现象。但是，通常情况下，无法求得光学Bloch方程和MB方程的解析解。只能借助于微分方程的数值解法来进行求解。这一章，讲述光学Blcoh方程、均匀展宽体系MB方程和非均匀展宽体系MB方程的数值解法。6.1光学Bloch方程的数值解法6.1.1常微分方程的数值解法常微分方程的经典数值解法有Euler法，Heun法和Runge－Kutta法。Euler法是最简单的有限差分法，误差是一阶，Heun法是Euler法的改进，也是最简单的预报–校正法，误差是二阶，Runge–Kutta是最为常用的求解方法，其中又以标准四阶Runge–Kutta法最为常用。Runge–Kutta法由于其格式和效率的限制，而且，标准四阶对于较小的步长需要大量计算，数值模拟时需要的时间较长。人们常常根据具体不同的数学物理方程，在标准格式的基本上提出修正的方法，如修正的四阶校正-预报法，该方法是以Euler法、Heun法为初始猜测，通过预报-校正二步来实现数值求解。对于如下形式常微分方程的初值问题，几种微分方程的基本原理d(,)dyftyt(6.1.1))0()0(yty(6.1.2)Euler法的计算步骤),(1kkkkythfyyhttkk1,1,,1,0Nk(6.1.3)式中，NTh，为步长，欧拉法是一阶格式的。Heun法的计算步骤是以Euler法作为预报，然后用梯形公式进行校正：),()0(1kkkkythfyy---预报格式(6.1.4)),(),(2)(11)1(1pkkkkkpkytfytfhyy---校正格式(6.1.5)式中，p为迭代次数，在步长选取适当的情况下，通常只需迭代两三次即可。标准RK4法由合适的泰勒方法推导而来，其计算步骤如下，6)22(43211ffffhyykk(6.1.6)132式中，),(1kkytff(6.1.7a))2,2(12fhyhtffkk(6.1.7b))2,2(23fhyhtffkk(6.1.7c)),(34fhyhtffkk(6.1.7d)标准RK4法对于较小的步长，需要大量的计算，并且当结果不够好时，必须重新计算。对于标准RK4法，已有了多种修正的方法。其中一种高效的算法是修正的四阶校正-预报法(PC4)。其计算步骤简述如下。第一步，用Heun法和中点规则，得到初始预报所需的初始点。第二步，用下面预报计算公式，计算预报点),(9),(38111)0(2/1kkkkkkytfytfhyy(6.1.8)),(10),(8),(432/12/1111)0(1kkkkkkkkytfytfytfhyy(6.1.9)第三步，用下面校正计算公式，计算得到结果),(),(8),(524)(11)(2/12/1)1(2/1pkkpkkkkkpkytfytfytfhyy(6.1.10)),(),(4),(6)(11)(2/12/1)1(1pkkpkkkkkpkytfytfytfhyy(6.1.11)6.1.2光学Bloch数值求解的结果我们首先将解析解与数值方法的结果进行比较。假设入射到二能级体系时矩形光脉冲，解析解I–IV与数值解计算结果如图6.1.1所示。从图中我们可以看出，除Euler法外，其余四种方法，其线型都与解析解重合，而Euler解也能在较小的误差范围内保持其线型，说明四种方法的计算结果是可靠的。计算有参数如下：光场的归一化面积均为3.2π，脉宽为50fs。步长为0.01fs。我们计算四种数值模拟方法的与解析解的绝对误差，如图6.1.2。从图中，我们可以看出，Euler法是步长一阶的，Heun法是步长二阶的，RK4法和PC4法都是步长四阶的，与这些算法的误差理论分析的结果是一致的。四种算法的误差都显示出一个规则的变化，但都随着积分时间的增长而增大，依然在误差的阶数量级内。RK4以及PC4，步长变成原来的十分之一，误差变成原来的万分之一，变小的幅度远小于Euler法以及Heun法。这直接验证了这四种算法在求解光学Bloch方程时，是可以保证收敛以及结果的可靠性。计算有参数如下：光场的归一化面积均为3.2π，脉宽为50fs。步长为0.01fs。133-1.0-0.50.00.51.0-1.0-0.50.00.51.001020304050-1.0-0.50.00.51.01020304050-1.0-0.50.00.51.0EulerwTime(fs)IIIIIIIV图6.1.1光学Bloch方程解析解I–IV和四种数值方法的计算结果0.02.0x10-24.0x10-26.0x10-28.0x10-21.0x10-10.02.0x10-44.0x10-46.0x10-48.0x10-4010203040500.04.0x10-98.0x10-91.2x10-81.6x10-82.0x10-810203040500.08.0x10-91.6x10-82.4x10-83.2x10-8AbsoluteerrorofwTime(fs)IIIIIIIVEulerHeunRK4PC4图6.1.2w分量四种数值解的相对误差矩形脉冲只是一种理想的近似，而与物质相互作用的超短脉冲都是具有一定形状的，通常是高斯与双曲正割型的。对于具有一定形状的脉冲，只有共振并忽略弛豫时间的情况下，才能得到解析解，表达为，0)(tu，)(sin)(ttv，)(cos)(ttw我们计算了高斯型和双曲正割型脉冲作用下，数值求解与解析解的不同数值结果，如图6.1.3，从图中，我们可以看出，由于脉冲形状的不同，其在脉冲相干区域的线型也不同，但四种数值算法的规律与矩形脉冲134是一样的。图6.1.3的上方是高斯型入射脉冲，下面的是双曲正割入射脉冲。同之前一样，其线型都与解析解重合，Euler也能在较小的误差范围内保持线型。计算有参数如下：光场的归一化面积均为3.2π，脉宽为50fs。步长为0.01fs。-1.0-0.50.00.51.002004006008001000-1.0-0.50.00.51.0wTime(fs)Euler图6.1.3高斯型和双曲正割脉冲作用下光学Bloch方程四种数值方法与解析解的比较我们计算在高斯型脉冲和双曲正割脉冲作用下，四种算法与解析解的绝对误差。其阶规律同矩形算法是一样的，但RK4法和PC4法的不再显示规则的变化，而是随机变化的，如图6.1.4所示。计算有参数如下：光场的归一化面积均为3.2π，脉宽为50fs。步长为0.01fs。0.02.0x10-24.0x10-20.02.0x10-34.0x10-36.0x10-3020040060080010000.02.0x10-720040060080010000.02.0x10-7AbsoluteerrowofwTime(fs)GaussSechIIIIIIIV图6.1.4高斯型和双曲正割脉冲作用下，光学Bloch方程四种数值方法与解析解的绝对误差从图中我们可以看出，Euler和Heun法的误差变换类似矩形脉冲，逐步累积，但是RK4以及PC4的误差不在显示规则变换，而是随机的。从以上计算结果我们得到，PC4法在求解微分方程上是可靠，可以应用135其来研究光学Bloch方程所描述的光学瞬态动力学过程。我们用PC4模拟一个算例，数值求解了没有任何近似的光学Bloch方程，即是失谐量，弛豫时间均不为零，如图6.1.5，计算时参数设置，光脉冲面积分别是π,2π,3.2π，脉冲形状是双曲正割型的，脉宽为50fs，9.0fs-1，1002Tfs，2001Tfs0.00.5-0.2-0.10.00.102004006008001000-1.0-0.8-0.620040060080010000.00.51.0pi2pi3.2piwTime(fs)uvWE图6.1.5在双曲正割脉冲作用下的光学瞬态相干动力学过程6.2均匀展宽体系MB方程的数值求解法6.2.1算法的建立偏微分方程的数值求解有多种方法，经典的标准算法：Euler法，即简单差分法，其误差是一阶的；4阶龙格-库塔法(RK4法，其误差是四阶的。均匀展宽体系MB方程的Euler法的计算公式为：1,11,,11,()/2ksksksksvh(6.2.1a)1,11,21,1,()ksksksksuuuvh(6.2.1b)1,11,21,1,11,()kskskskskksvvvuwh(6.2.1c)1,11,11,1,1,[(1))ksksksksks(6.2.1d)式中，hzzss1，httkk1，1,...,1,0,Nks。NTh/，为步长。RK4法计算公式为：1361,11,,11234[(22)/6]/2ksksksqqqqh(6.2.2a)1,11,1234(22)/6ksksuuffffh(6.2.2b)1,11,1234(22)/6ksksvvggggh(6.2.2c)1,11,1234(22)/6kskswwpppph(6.2.2d)式中，),(,11skkutff(6.2.3a)21,1(,)22kkshhfftuf(6.2.3b)31,2(,)22kkshhfftuf(6.2.3c)41,3(,)kksffthuhf(6.2.3d)),,(,1,11skskkwtgg(6.2.3e)21,11,1(,,)222kkskshhhggtqwp(6.2.3f)31,21,2(,,)222kkskshhhggtqwp(6.2.3i)41,31,3(,,)kksksggthhqwhp(6.2.3j)),,(,1,11skskkvtpp(6.2.3k)21,11,1(,,)222kkskshhhpptqvg(6.2.3l)31,21,2(,,)222kkskshhhpptqvg(6.2.3m)41,31,3(,,)kksksppthhqvhg(6.2.3n)),(,11skkvtqq(6.2.3o)21,1(,)22kkshhqqtvg(6.2.3p)31,2(,)22kkshhqqtvg(6.2.3q)41,3(,)kksqqthvhg(6.2.3r)经研究表明，上述两方法求解MB的可靠性并不好。为了建立可靠、稳定、高效的数值算法。我们建立了“预报校正-龙格-库塔(PCRK4)”算法。137均匀展宽MB方程中的Bloch方程部分，这部分可以看成是时间的一阶常微分方程组，对应于每一个空间格点的u、v、w可用标准RK4法进行求解。其格式如下：1234(,1)(,)(22)6tuuuuzzttzztthffffuu(6.2.4a)1234(,1)(,)(22)6tvvvvzzttzztthffffvv(6.2.4b)1234(,1)(,)(22)6t(6.2.4c)定义2ufuv，2vfvuw，1wfwu分别为Bloch方程的右边，th为时间格点步长。Bloch方程中各项12(,)(,)uzzttzztttfuvh(6.2.5a)22(,)1(,)1()()22ttuzzttuzzttvthhfufvfh(6.2.5b)32(,)2(,)2()()2